/***********************************************************************
**                                                                    **
**          /----------------------------------------------\          **
**         |  Delphes, a framework for the fast simulation  |         **
**         |       of a  generic collider experiment        |         **
**          \-------------  arXiv:0903.2225v1  ------------/          **
**                                                                    **
**                                                                    **
**   This package uses:                                               **
**   ------------------                                               **
**   ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86        **
**   FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210]      **
**   Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]       **
**   FROG: [hep-ex/0901.2718v1]                                       **
**   HepMC: Comput. Phys. Commun.134 (2001) 41                        **
**                                                                    **
** ------------------------------------------------------------------ **
**                                                                    **
**   Main authors:                                                    **
**   -------------                                                    **
**                                                                    **
**                Severine Ovyn                Xavier Rouby           **
**          severine.ovyn@uclouvain.be      xavier.rouby@cern         **
**                                                                    **
**         Center for Particle Physics and Phenomenology (CP3)        **
**                Universite catholique de Louvain (UCL)              **
**                     Louvain-la-Neuve, Belgium                      **
**                                                                    **
**                      Copyright (C) 2008-2010,                      **
**                        All rights reserved.                        **
**                                                                    **
***********************************************************************/

/// \file SmearUtil.cc
/// \brief RESOLution class, and some generic definitions


#include "SmearUtil.h"
#include "TStopwatch.h"
#include "TF2.h"

#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <map>
#include <vector>
#include <cmath>
#include <cstdlib> // for exit()
using namespace std;

//------------------------------------------------------------------------------

RESOLution::RESOLution() {

  // Detector characteristics
  CEN_max_tracker   = 2.5;                // Maximum tracker coverage          
  CEN_max_calo_cen  = 1.7;                // central calorimeter coverage
  CEN_max_calo_ec   = 3.0;                // calorimeter endcap coverage
  CEN_max_calo_fwd  = 5.0;                // forward calorimeter pseudorapidity coverage
  CEN_max_mu        = 2.4;                // muon chambers pseudorapidity coverage
  
  // Energy resolution for electron/photon
  // \sigma/E = C + N/E + S/\sqrt{E}
  ELG_Scen          = 0.05;              // S term for central ECAL
  ELG_Ncen          = 0.25;              // N term 
  ELG_Ccen          = 0.005;             // C term 
  ELG_Sec           = 0.05;              // S term for central ECAL endcap
  ELG_Nec           = 0.25;              // S term 
  ELG_Cec           = 0.005;              // S term 
  ELG_Sfwd          = 2.084;             // S term for FCAL
  ELG_Nfwd          = 0.0;               // N term 
  ELG_Cfwd          = 0.107;             // C term 
  ELG_Szdc          = 0.70;              // S term for ZDC
  ELG_Nzdc          = 0.0;               // N term 
  ELG_Czdc          = 0.08;              // C term

  // Energy resolution for hadrons in ecal/hcal/fwd
  // \sigma/E = C + N/E + S/\sqrt{E}
  HAD_Scen          = 1.5;               // S term for central HCAL
  HAD_Ncen          = 0.;                // N term
  HAD_Ccen          = 0.05;              // C term
  HAD_Sec           = 1.5;               // S term for HCAL endcap
  HAD_Nec           = 0.;                // N term 
  HAD_Cec           = 0.05;              // C term 
  HAD_Sfwd          = 2.7;               // S term for FCAL 
  HAD_Nfwd          = 0.;                // N term
  HAD_Cfwd          = 0.13;              // C term
  HAD_Szdc          = 1.38;              // S term for ZDC
  HAD_Nzdc          = 0.;                // N term
  HAD_Czdc          = 0.13;              // C term

  // Muon smearing
  MU_SmearPt        = 0.01;

  // time resolution
  ZDC_T_resolution   = 0;	// resolution for time measurement [s]
  RP220_T_resolution = 0;
  RP420_T_resolution = 0;

  // Tracking efficiencies
  TRACK_ptmin       = 0.9;               // minimal pt needed to reach the calorimeter in GeV
  TRACK_eff         = 100;               // efficiency associated to the tracking

  // Calorimetric towers
  TOWER_number      = 40;
  const float tower_eta_edges[41] = { 
    0., 0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609, 0.696, 0.783, 0.870, 0.957, 1.044, 1.131, 1.218, 1.305, 1.392, 1.479, 1.566,
    1.653, 1.740, 1.830, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.868, 2.950, 3.125, 3.300, 3.475, 3.650, 3.825, 4.000, 4.175,
    4.350, 4.525,  4.700, 5.000}; // temporary object
  TOWER_eta_edges = new float[TOWER_number+1];
  for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = tower_eta_edges[i];
  
  const float tower_dphi[40] = { 
    5, 5, 5, 5, 5,  5, 5, 5, 5, 5,  5, 5, 5, 5, 5,  5, 5, 5, 5, 10,
    10,10,10,10,10, 10,10,10,10,10, 10,10,10,10,10, 10,10,10,20, 20 }; // temporary object
  TOWER_dphi = new float[TOWER_number];
  for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = tower_dphi[i];


  // Thresholds for reconstructed objetcs (GeV)
  PTCUT_elec      = 10.0;
  PTCUT_muon      = 10.0;
  PTCUT_jet       = 20.0;
  PTCUT_gamma     = 10.0;
  PTCUT_taujet    = 10.0;

  ZDC_gamma_E     = 20; // GeV
  ZDC_n_E	  = 50; // GeV

  // Isolation
  ISOL_trk_PT         = 2.0;      //minimal pt of tracks for isolation criteria
  ISOL_trk_Cone   = 0.5;      //Cone  for isolation criteria
  ISOL_calo_Cone  = 0.5;      //Cone  for isolation criteria
  ISOL_calo_ET    = 1E99;     //minimal tower energy for isolation criteria. Default off = 1E99
  ISOL_calo_Grid  = 3;        //Grid size (N x N) for calorimetric isolation -- should be odd

  // General jet variable
  JET_coneradius   = 0.7;	     // generic jet radius ; not for tau's !!!
  JET_jetalgo      = 2;              // 1 for Cone algorithm, 2 for MidPoint algorithm, 3 for SIScone algorithm, 4 for kt algorithm
  JET_seed         = 1.0;            // minimum seed to start jet reconstruction
  JET_Eflow        = 0;		     // 1 for Energy flow in jets reco ; 0 if not

  // Tagging definition
  BTAG_func_b      = "0.4  + 0*x + 0*y";	// efficiency functions (x=Pt, y=eta)  in ROOT::TF2 format
  BTAG_func_c      = "0.1  + 0*x + 0*y";
  BTAG_func_l      = "0.01 + 0*x + 0*y";
  BTAG_func_g      = "0.01 + 0*x + 0*y";


  // FLAGS
  FLAG_bfield      = 1;                       //1 to run the bfield propagation else 0
  FLAG_vfd         = 1;                       //1 to run the very forward detectors else 0
  FLAG_RP          = 1;                       //1 to run the zero degree calorimeter else 0
  FLAG_trigger     = 1;                       //1 to run the trigger selection else 0
  FLAG_frog        = 1;                       //1 to run the FROG event display
  FLAG_lhco        = 1;

  // In case BField propagation allowed
  TRACK_radius      = 129;                   //radius of the BField coverage
  TRACK_length      = 300;                   //length of the BField coverage
  TRACK_bfield_x    = 0;                     //X composant of the BField
  TRACK_bfield_y    = 0;                     //Y composant of the BField
  TRACK_bfield_z    = 3.8;                   //Z composant of the BField

  // In case Very forward detectors allowed
  VFD_min_calo_vfd  = 5.2;                   // very forward calorimeter (if any) like CASTOR
  VFD_max_calo_vfd  = 6.6;
  VFD_min_zdc       = 8.3;
  VFD_s_zdc         = 140;                   // distance of the Zero Degree Calorimeter, from the Interaction poin, in [m]

  RP_220_s          = 220;                   // distance of the RP to the IP, in meters
  RP_220_x          = 0.002;                 // distance of the RP to the beam, in meters
  RP_420_s          = 420;                   // distance of the RP to the IP, in meters
  RP_420_x          = 0.004;                 // distance of the RP to the beam, in meters
  RP_IP_name        = "IP5"; 
  RP_beam1Card      = "data/LHCB1IR5_v6.500.tfs";
  RP_beam2Card      = "data/LHCB1IR5_v6.500.tfs";

  // In case FROG event display allowed
  NEvents_Frog      = 10;

  // Number of events to be processed 
  NEvents           = -1;

  //********************************************
  //jet stuffs not defined in the input datacard
  //********************************************
  
  JET_overlap     = 0.75;
  // MidPoint algorithm definition
  JET_M_coneareafraction = 0.25;
  JET_M_maxpairsize      = 2;
  JET_M_maxiterations    = 100;
  // Define Cone algorithm.
  JET_C_adjacencycut  =  2;
  JET_C_maxiterations = 100;
  JET_C_iratch        =  1;
  //Define SISCone algorithm.
  JET_S_npass         = 0;
  JET_S_protojet_ptmin= 0.0;
  
  //For Tau-jet definition
  TAU_energy_scone = 0.15;  // radius R of the cone for tau definition, based on energy threshold
  TAU_track_scone  = 0.4;   // radius R of the cone for tau definition, based on track number
  TAU_track_pt     = 2;     // minimal pt [GeV] for tracks to be considered in tau definition
  TAU_energy_frac  = 0.95;  // fraction of energy required in the central part of the cone, for tau jets
  
  PT_QUARKS_MIN =   2.0  ;                   // minimal pt needed by quarks to do b-tag

  //for very forward detectors
  RP_offsetEl_s     = 120;    // distance of beam separation point, from IP
  RP_offsetEl_x     = -0.097; // half distance of separation of beams in horizontal plan, in m
  RP_offsetEl_y     = 0;      // half distance of separation of beams in vertical plan, in m
  RP_cross_x        = -500;   // IP offset in horizontal plane, in micrometers
  RP_cross_y        = 0.0;    // IP offset in vertical plane, in micrometers
  RP_cross_ang_x    = 142.5;  // half-crossing angle in horizontal plane, in microrad
  RP_cross_ang_y    = 0.0;    // half-crossing angle in vertical plane, in microrad


  PdgTableFilename  = "data/particle.tbl";
  inputfilelist     = "";
  detectorcard      = "";
  triggercard       = "";

  grandom = new TRandom3(0);  // a new seed is set everytime Delphes is run
}


RESOLution::RESOLution(const RESOLution & DET) {
  // Detector characteristics
  CEN_max_tracker   = DET.CEN_max_tracker; 
  CEN_max_calo_cen  = DET.CEN_max_calo_cen;
  CEN_max_calo_ec   = DET.CEN_max_calo_ec;
  CEN_max_calo_fwd  = DET.CEN_max_calo_fwd;
  CEN_max_mu        = DET.CEN_max_mu;
  
  // Energy resolution for electron/photon
  ELG_Scen          = DET.ELG_Scen;
  ELG_Ncen          = DET.ELG_Ncen;
  ELG_Ccen          = DET.ELG_Ccen;
  ELG_Sec           = DET.ELG_Sec;
  ELG_Nec           = DET.ELG_Nec;
  ELG_Cec           = DET.ELG_Cec;
  ELG_Cfwd          = DET.ELG_Cfwd;
  ELG_Sfwd          = DET.ELG_Sfwd;
  ELG_Nfwd          = DET.ELG_Nfwd;
  ELG_Czdc          = DET.ELG_Czdc;
  ELG_Szdc          = DET.ELG_Szdc;
  ELG_Nzdc          = DET.ELG_Nzdc;

  // Energy resolution for hadrons in ecal/hcal/fwd/zdc
  HAD_Scen          = DET.HAD_Scen;
  HAD_Ncen          = DET.HAD_Ncen;
  HAD_Ccen          = DET.HAD_Ccen;
  HAD_Sec           = DET.HAD_Sec;
  HAD_Nec           = DET.HAD_Nec;
  HAD_Cec           = DET.HAD_Cec;
  HAD_Sfwd           = DET.HAD_Sfwd;
  HAD_Nfwd           = DET.HAD_Nfwd;
  HAD_Cfwd           = DET.HAD_Cfwd;
  HAD_Szdc          = DET.HAD_Szdc;
  HAD_Nzdc          = DET.HAD_Nzdc;
  HAD_Czdc          = DET.HAD_Czdc;

  // time resolution
  ZDC_T_resolution   = DET.ZDC_T_resolution;       // resolution for time measurement [s]
  RP220_T_resolution = DET.RP220_T_resolution;
  RP420_T_resolution = DET.RP420_T_resolution;

  // Muon smearing
  MU_SmearPt        = DET.MU_SmearPt;

  // Tracking efficiencies
  TRACK_ptmin       = DET.TRACK_ptmin;
  TRACK_eff         = DET.TRACK_eff;

  // Calorimetric towers
  TOWER_number      = DET.TOWER_number;
  TOWER_eta_edges = new float[TOWER_number+1];
  for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = DET.TOWER_eta_edges[i];

  TOWER_dphi = new float[TOWER_number];
  for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = DET.TOWER_dphi[i];

  // Thresholds for reconstructed objetcs
  PTCUT_elec      = DET.PTCUT_elec;
  PTCUT_muon      = DET.PTCUT_muon;
  PTCUT_jet       = DET.PTCUT_jet;
  PTCUT_gamma     = DET.PTCUT_gamma;
  PTCUT_taujet    = DET.PTCUT_taujet;

  ZDC_gamma_E     = DET.ZDC_gamma_E;
  ZDC_n_E         = DET.ZDC_n_E;

  // Isolation
  ISOL_trk_PT     = DET.ISOL_trk_PT;     // tracking isolation
  ISOL_trk_Cone   = DET.ISOL_trk_Cone;  
  ISOL_calo_Cone  = DET.ISOL_calo_Cone;  
  ISOL_calo_ET    = DET.ISOL_calo_ET;  // calorimeter isolation, defaut off
  ISOL_calo_Grid  = DET.ISOL_calo_Grid; 


  // General jet variable
  JET_coneradius   = DET.JET_coneradius;
  JET_jetalgo      = DET.JET_jetalgo;
  JET_seed         = DET.JET_seed;
  JET_Eflow        = DET.JET_Eflow;

  // Tagging definition
  BTAG_func_b      = DET.BTAG_func_b;
  BTAG_func_c      = DET.BTAG_func_c;
  BTAG_func_l      = DET.BTAG_func_l;
  BTAG_func_g      = DET.BTAG_func_g;


  // FLAGS
  FLAG_bfield      = DET.FLAG_bfield;
  FLAG_vfd         = DET.FLAG_vfd;
  FLAG_RP          = DET.FLAG_RP;
  FLAG_trigger     = DET.FLAG_trigger;
  FLAG_frog        = DET.FLAG_frog;
  FLAG_lhco        = DET.FLAG_lhco;

  // In case BField propagation allowed
  TRACK_radius      = DET.TRACK_radius;
  TRACK_length      = DET.TRACK_length;
  TRACK_bfield_x    = DET.TRACK_bfield_x;
  TRACK_bfield_y    = DET.TRACK_bfield_y;
  TRACK_bfield_z    = DET.TRACK_bfield_z;

  // In case Very forward detectors allowed
  VFD_min_calo_vfd  = DET.VFD_min_calo_vfd;
  VFD_max_calo_vfd  = DET.VFD_max_calo_vfd;
  VFD_min_zdc       = DET.VFD_min_zdc;
  VFD_s_zdc         = DET.VFD_s_zdc;

  RP_220_s          = DET.RP_220_s;
  RP_220_x          = DET.RP_220_x;
  RP_420_s          = DET.RP_420_s;
  RP_420_x          = DET.RP_420_x;
  RP_beam1Card      = DET.RP_beam1Card;
  RP_beam2Card      = DET.RP_beam2Card;
  RP_offsetEl_s     = DET.RP_offsetEl_s;
  RP_offsetEl_x     = DET.RP_offsetEl_x;
  RP_offsetEl_y     = DET.RP_offsetEl_y;
  RP_cross_x        = DET.RP_cross_x;
  RP_cross_y        = DET.RP_cross_y;
  RP_cross_ang_x    = DET.RP_cross_ang_x;
  RP_cross_ang_y    = DET.RP_cross_ang_y;
  RP_IP_name        = DET.RP_IP_name;

  // In case FROG event display allowed
  NEvents_Frog      = DET.NEvents_Frog;

  // Number of events to be processed
  NEvents           = DET.NEvents;

  JET_overlap      = DET.JET_overlap;
  // MidPoint algorithm definition
  JET_M_coneareafraction = DET.JET_M_coneareafraction;
  JET_M_maxpairsize      = DET.JET_M_maxpairsize;
  JET_M_maxiterations    = DET.JET_M_maxiterations;
  // Define Cone algorithm.
  JET_C_adjacencycut     = DET.JET_C_adjacencycut;
  JET_C_maxiterations    = DET.JET_C_maxiterations;
  JET_C_iratch           = DET.JET_C_iratch;
  //Define SISCone algorithm.
  JET_S_npass            = DET.JET_S_npass;
  JET_S_protojet_ptmin   = DET.JET_S_protojet_ptmin;

  //For Tau-jet definition
  TAU_energy_scone = DET.TAU_energy_scone;
  TAU_track_scone  = DET.TAU_track_scone;
  TAU_track_pt     = DET.TAU_track_pt;
  TAU_energy_frac  = DET.TAU_energy_frac;

  PT_QUARKS_MIN    = DET.PT_QUARKS_MIN;
  PdgTableFilename = DET.PdgTableFilename;
  PDGtable         = DET.PDGtable;
  inputfilelist    = DET.inputfilelist;
  detectorcard     = DET.detectorcard;
  triggercard      = DET.triggercard;

  grandom          = new TRandom3(*(DET.grandom));
}

RESOLution& RESOLution::operator=(const RESOLution& DET) {
  if(this==&DET) return *this;
  // Detector characteristics
  CEN_max_tracker   = DET.CEN_max_tracker; 
  CEN_max_calo_cen  = DET.CEN_max_calo_cen;
  CEN_max_calo_ec   = DET.CEN_max_calo_ec;
  CEN_max_calo_fwd  = DET.CEN_max_calo_fwd;
  CEN_max_mu        = DET.CEN_max_mu;
  
  // Energy resolution for electron/photon
  ELG_Scen          = DET.ELG_Scen;
  ELG_Ncen          = DET.ELG_Ncen;
  ELG_Ccen          = DET.ELG_Ccen;
  ELG_Sec           = DET.ELG_Sec;
  ELG_Nec           = DET.ELG_Nec;
  ELG_Cec           = DET.ELG_Cec;
  ELG_Cfwd          = DET.ELG_Cfwd;
  ELG_Sfwd          = DET.ELG_Sfwd;
  ELG_Nfwd          = DET.ELG_Nfwd;
  ELG_Czdc          = DET.ELG_Czdc;
  ELG_Szdc          = DET.ELG_Szdc;
  ELG_Nzdc          = DET.ELG_Nzdc;

  // Energy resolution for hadrons in ecal/hcal/fwd/zdc
  HAD_Scen          = DET.HAD_Scen ;
  HAD_Ncen          = DET.HAD_Ncen;
  HAD_Ccen          = DET.HAD_Ccen;
  HAD_Sec           = DET.HAD_Sec;
  HAD_Nec           = DET.HAD_Nec;
  HAD_Cec           = DET.HAD_Cec;
  HAD_Sfwd          = DET.HAD_Sfwd;
  HAD_Nfwd          = DET.HAD_Nfwd;
  HAD_Cfwd          = DET.HAD_Cfwd;
  HAD_Szdc          = DET.HAD_Szdc;
  HAD_Nzdc          = DET.HAD_Nzdc;
  HAD_Czdc          = DET.HAD_Czdc;

  // time resolution
  ZDC_T_resolution   = DET.ZDC_T_resolution;       // resolution for time measurement [s]
  RP220_T_resolution = DET.RP220_T_resolution;
  RP420_T_resolution = DET.RP420_T_resolution;

  // Muon smearing
  MU_SmearPt        = DET.MU_SmearPt;

  // Tracking efficiencies
  TRACK_ptmin       = DET.TRACK_ptmin;
  TRACK_eff         = DET.TRACK_eff;

  // Calorimetric towers
  TOWER_number      = DET.TOWER_number;
  TOWER_eta_edges = new float[TOWER_number+1];
  for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = DET.TOWER_eta_edges[i];

  TOWER_dphi = new float[TOWER_number];
  for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = DET.TOWER_dphi[i];

  // Thresholds for reconstructed objetcs
  PTCUT_elec      = DET.PTCUT_elec;
  PTCUT_muon      = DET.PTCUT_muon;
  PTCUT_jet       = DET.PTCUT_jet;
  PTCUT_gamma     = DET.PTCUT_gamma;
  PTCUT_taujet    = DET.PTCUT_taujet;

  ZDC_gamma_E     = DET.ZDC_gamma_E;
  ZDC_n_E         = DET.ZDC_n_E;

  // Isolation
  ISOL_trk_PT     = DET.ISOL_trk_PT;       // tracking isolation
  ISOL_trk_Cone   = DET.ISOL_trk_Cone;    
  ISOL_calo_Cone  = DET.ISOL_calo_Cone;    
  ISOL_calo_ET    = DET.ISOL_calo_ET;  // calorimeter isolation, defaut off
  ISOL_calo_Grid  = DET.ISOL_calo_Grid;

  // General jet variable
  JET_coneradius   = DET.JET_coneradius;
  JET_jetalgo      = DET.JET_jetalgo;
  JET_seed         = DET.JET_seed;
  JET_Eflow        = DET.JET_Eflow;

  // Tagging definition
  BTAG_func_b      = DET.BTAG_func_b;
  BTAG_func_c      = DET.BTAG_func_c;
  BTAG_func_l      = DET.BTAG_func_l;
  BTAG_func_g      = DET.BTAG_func_g;

  // FLAGS
  FLAG_bfield      = DET.FLAG_bfield;
  FLAG_vfd         = DET.FLAG_vfd;
  FLAG_RP          = DET.FLAG_RP;
  FLAG_trigger     = DET.FLAG_trigger;
  FLAG_frog        = DET.FLAG_frog;
  FLAG_lhco        = DET.FLAG_lhco;

  // In case BField propagation allowed
  TRACK_radius      = DET.TRACK_radius;
  TRACK_length      = DET.TRACK_length;
  TRACK_bfield_x    = DET.TRACK_bfield_x;
  TRACK_bfield_y    = DET.TRACK_bfield_y;
  TRACK_bfield_z    = DET.TRACK_bfield_z;

  // In case Very forward detectors allowed
  VFD_min_calo_vfd  = DET.VFD_min_calo_vfd;
  VFD_max_calo_vfd  = DET.VFD_max_calo_vfd;
  VFD_min_zdc       = DET.VFD_min_zdc;
  VFD_s_zdc         = DET.VFD_s_zdc;

  RP_220_s          = DET.RP_220_s;
  RP_220_x          = DET.RP_220_x;
  RP_420_s          = DET.RP_420_s;
  RP_420_x          = DET.RP_420_x;
  RP_offsetEl_s     = DET.RP_offsetEl_s;
  RP_offsetEl_x     = DET.RP_offsetEl_x;
  RP_offsetEl_y     = DET.RP_offsetEl_y;
  RP_beam1Card      = DET.RP_beam1Card;
  RP_beam2Card      = DET.RP_beam2Card;
  RP_cross_x        = DET.RP_cross_x;
  RP_cross_y        = DET.RP_cross_y;
  RP_cross_ang_x    = DET.RP_cross_ang_x;
  RP_cross_ang_y    = DET.RP_cross_ang_y;
  RP_IP_name        = DET.RP_IP_name;


  // In case FROG event display allowed
  NEvents_Frog      = DET.NEvents_Frog;

  // Number of events to be processed
  NEvents           = DET.NEvents;

  JET_overlap      = DET.JET_overlap;
  // MidPoint algorithm definition
  JET_M_coneareafraction = DET.JET_M_coneareafraction;
  JET_M_maxpairsize      = DET.JET_M_maxpairsize;
  JET_M_maxiterations    = DET.JET_M_maxiterations;
  // Define Cone algorithm.
  JET_C_adjacencycut     = DET.JET_C_adjacencycut;
  JET_C_maxiterations    = DET.JET_C_maxiterations;
  JET_C_iratch           = DET.JET_C_iratch;
  //Define SISCone algorithm.
  JET_S_npass            = DET.JET_S_npass;
  JET_S_protojet_ptmin   = DET.JET_S_protojet_ptmin;

  //For Tau-jet definition
  TAU_energy_scone = DET.TAU_energy_scone;
  TAU_track_scone  = DET.TAU_track_scone;
  TAU_track_pt     = DET.TAU_track_pt;
  TAU_energy_frac  = DET.TAU_energy_frac;

  PT_QUARKS_MIN    = DET.PT_QUARKS_MIN;

  PdgTableFilename = DET.PdgTableFilename;
  PDGtable	   = DET.PDGtable;

  inputfilelist    = DET.inputfilelist;
  detectorcard     = DET.detectorcard;
  triggercard      = DET.triggercard;

  grandom          = new TRandom3(*(DET.grandom));

  return *this;
}

void RESOLution::setNames(const string& list, const string& det, const string& trig) {
  inputfilelist    = list;
  detectorcard     = det;
  triggercard      = trig;
}

//------------------------------------------------------------------------------
void RESOLution::ReadDataCard(const string datacard) {
  
  string temp_string;
  istringstream curstring;
  
  ifstream fichier_a_lire(datacard.c_str());
  if(!fichier_a_lire.good()) {
    cout <<"**        WARNING: Datadard  not found, use default values         **" << endl;
    return;
  }
  bool CEN_max_calo_ec_flag = false;
  
  while (getline(fichier_a_lire,temp_string)) {
    curstring.clear(); // needed when using several times istringstream::str(string)
    curstring.str(temp_string);
    string varname;
    float value; int ivalue; string svalue;

    if(strstr(temp_string.c_str(),"#")) { }
    else if(strstr(temp_string.c_str(),"CEN_max_tracker"))  {curstring >> varname >> value; CEN_max_tracker   = value;}
    else if(strstr(temp_string.c_str(),"CEN_max_calo_cen")) {curstring >> varname >> value; CEN_max_calo_cen  = value;}
    else if(strstr(temp_string.c_str(),"CEN_max_calo_ec"))  {CEN_max_calo_ec_flag=true; curstring >> varname >> value; CEN_max_calo_ec   = value;}
    else if(strstr(temp_string.c_str(),"CEN_max_calo_fwd")) {curstring >> varname >> value; CEN_max_calo_fwd  = value;}
    else if(strstr(temp_string.c_str(),"CEN_max_mu"))       {curstring >> varname >> value; CEN_max_mu        = value;}

    else if(strstr(temp_string.c_str(),"VFD_min_calo_vfd")) {curstring >> varname >> value; VFD_min_calo_vfd  = value;}
    else if(strstr(temp_string.c_str(),"VFD_max_calo_vfd")) {curstring >> varname >> value; VFD_max_calo_vfd  = value;}
    else if(strstr(temp_string.c_str(),"VFD_min_zdc"))      {curstring >> varname >> value; VFD_min_zdc       = value;}
    else if(strstr(temp_string.c_str(),"VFD_s_zdc"))        {curstring >> varname >> value; VFD_s_zdc         = value;}
    
    else if(strstr(temp_string.c_str(),"RP_220_s"))         {curstring >> varname >> value; RP_220_s          = value;}
    else if(strstr(temp_string.c_str(),"RP_220_x"))         {curstring >> varname >> value; RP_220_x          = value;}
    else if(strstr(temp_string.c_str(),"RP_420_s"))         {curstring >> varname >> value; RP_420_s          = value;}
    else if(strstr(temp_string.c_str(),"RP_420_x"))         {curstring >> varname >> value; RP_420_x          = value;}
    else if(strstr(temp_string.c_str(),"RP_beam1Card"))     {curstring >> varname >> svalue;RP_beam1Card      = svalue;}
    else if(strstr(temp_string.c_str(),"RP_beam2Card"))     {curstring >> varname >> svalue;RP_beam2Card      = svalue;}
    else if(strstr(temp_string.c_str(),"RP_IP_name"))       {curstring >> varname >> svalue;RP_IP_name        = svalue;}

    else if(strstr(temp_string.c_str(),"RP_offsetEl_s"))    {curstring >> varname >> value; RP_offsetEl_s     = value;}
    else if(strstr(temp_string.c_str(),"RP_offsetEl_x"))    {curstring >> varname >> value; RP_offsetEl_x     = value;}
    else if(strstr(temp_string.c_str(),"RP_offsetEl_y"))    {curstring >> varname >> value; RP_offsetEl_y     = value;}
    else if(strstr(temp_string.c_str(),"RP_cross_x"))       {curstring >> varname >> value; RP_cross_x        = value;}
    else if(strstr(temp_string.c_str(),"RP_cross_y"))       {curstring >> varname >> value; RP_cross_y        = value;}
    else if(strstr(temp_string.c_str(),"RP_cross_ang_x"))   {curstring >> varname >> value; RP_cross_ang_x    = value;}
    else if(strstr(temp_string.c_str(),"RP_cross_ang_y"))   {curstring >> varname >> value; RP_cross_ang_y    = value;}
    
    else if(strstr(temp_string.c_str(),"ELG_Scen"))         {curstring >> varname >> value; ELG_Scen          = value;}
    else if(strstr(temp_string.c_str(),"ELG_Ncen"))         {curstring >> varname >> value; ELG_Ncen          = value;}
    else if(strstr(temp_string.c_str(),"ELG_Ccen"))         {curstring >> varname >> value; ELG_Ccen          = value;}
    else if(strstr(temp_string.c_str(),"ELG_Sec"))          {curstring >> varname >> value; ELG_Sec           = value;}
    else if(strstr(temp_string.c_str(),"ELG_Nec"))          {curstring >> varname >> value; ELG_Nec           = value;}
    else if(strstr(temp_string.c_str(),"ELG_Cec"))          {curstring >> varname >> value; ELG_Cec           = value;}
    else if(strstr(temp_string.c_str(),"ELG_Sfwd"))         {curstring >> varname >> value; ELG_Sfwd          = value;}
    else if(strstr(temp_string.c_str(),"ELG_Cfwd"))         {curstring >> varname >> value; ELG_Cfwd          = value;}
    else if(strstr(temp_string.c_str(),"ELG_Nfwd"))         {curstring >> varname >> value; ELG_Nfwd          = value;}
    else if(strstr(temp_string.c_str(),"ELG_Szdc"))         {curstring >> varname >> value; ELG_Szdc          = value;}
    else if(strstr(temp_string.c_str(),"ELG_Czdc"))         {curstring >> varname >> value; ELG_Czdc          = value;}
    else if(strstr(temp_string.c_str(),"ELG_Nzdc"))         {curstring >> varname >> value; ELG_Nzdc          = value;}

    else if(strstr(temp_string.c_str(),"HAD_Shcal"))        {warning("HAD_Shcal","HAD_Scen"); curstring >> varname >> value; HAD_Scen          = value;}
    else if(strstr(temp_string.c_str(),"HAD_Nhcal"))        {warning("HAD_Nhcal","HAD_Ncen"); curstring >> varname >> value; HAD_Ncen          = value;}
    else if(strstr(temp_string.c_str(),"HAD_Chcal"))        {warning("HAD_Chcal","HAD_Ccen"); curstring >> varname >> value; HAD_Ccen          = value;}
    else if(strstr(temp_string.c_str(),"HAD_Shf"))          {warning("HAD_Shf","HAD_Sfwd"); curstring >> varname >> value; HAD_Sfwd          = value;}
    else if(strstr(temp_string.c_str(),"HAD_Nhf"))          {warning("HAD_Nhf","HAD_Nfwd"); curstring >> varname >> value; HAD_Nfwd          = value;}
    else if(strstr(temp_string.c_str(),"HAD_Chf"))          {warning("HAD_Chf","HAD_Cfwd"); curstring >> varname >> value; HAD_Cfwd          = value;}



    else if(strstr(temp_string.c_str(),"HAD_Scen"))         {curstring >> varname >> value; HAD_Scen          = value;}
    else if(strstr(temp_string.c_str(),"HAD_Ncen"))         {curstring >> varname >> value; HAD_Ncen          = value;}
    else if(strstr(temp_string.c_str(),"HAD_Ccen"))         {curstring >> varname >> value; HAD_Ccen          = value;}
    else if(strstr(temp_string.c_str(),"HAD_Sec"))          {curstring >> varname >> value; HAD_Sec           = value;}
    else if(strstr(temp_string.c_str(),"HAD_Nec"))          {curstring >> varname >> value; HAD_Nec           = value;}
    else if(strstr(temp_string.c_str(),"HAD_Cec"))          {curstring >> varname >> value; HAD_Cec           = value;}
    else if(strstr(temp_string.c_str(),"HAD_Sfwd"))         {curstring >> varname >> value; HAD_Sfwd          = value;}
    else if(strstr(temp_string.c_str(),"HAD_Nfwd"))         {curstring >> varname >> value; HAD_Nfwd          = value;}
    else if(strstr(temp_string.c_str(),"HAD_Cfwd"))         {curstring >> varname >> value; HAD_Cfwd          = value;}
    else if(strstr(temp_string.c_str(),"HAD_Szdc"))         {curstring >> varname >> value; HAD_Szdc          = value;}
    else if(strstr(temp_string.c_str(),"HAD_Nzdc"))         {curstring >> varname >> value; HAD_Nzdc          = value;}
    else if(strstr(temp_string.c_str(),"HAD_Czdc"))         {curstring >> varname >> value; HAD_Czdc          = value;}

    else if(strstr(temp_string.c_str(),"ZDC_T_resolution")) {curstring >> varname >> value; ZDC_T_resolution  = value;}
    else if(strstr(temp_string.c_str(),"RP220_T_resolution")) {curstring >> varname >> value; RP220_T_resolution  = value;}
    else if(strstr(temp_string.c_str(),"RP420_T_resolution")) {curstring >> varname >> value; RP420_T_resolution  = value;}
    else if(strstr(temp_string.c_str(),"MU_SmearPt"))       {curstring >> varname >> value; MU_SmearPt        = value;}
    
    else if(strstr(temp_string.c_str(),"TRACK_radius"))     {curstring >> varname >> ivalue;TRACK_radius      = ivalue;}
    else if(strstr(temp_string.c_str(),"TRACK_length"))     {curstring >> varname >> ivalue;TRACK_length      = ivalue;}
    else if(strstr(temp_string.c_str(),"TRACK_bfield_x"))   {curstring >> varname >> value; TRACK_bfield_x    = value;}
    else if(strstr(temp_string.c_str(),"TRACK_bfield_y"))   {curstring >> varname >> value; TRACK_bfield_y    = value;}
    else if(strstr(temp_string.c_str(),"TRACK_bfield_z"))   {curstring >> varname >> value; TRACK_bfield_z    = value;}
    else if(strstr(temp_string.c_str(),"FLAG_bfield"))      {curstring >> varname >> ivalue; FLAG_bfield      = ivalue;}
    else if(strstr(temp_string.c_str(),"TRACK_ptmin"))      {curstring >> varname >> value; TRACK_ptmin       = value;}
    else if(strstr(temp_string.c_str(),"TRACK_eff"))        {curstring >> varname >> value; TRACK_eff         = value;}

    else if(strstr(temp_string.c_str(),"TOWER_number"))     {curstring >> varname >> ivalue;TOWER_number      = ivalue;}
    else if(strstr(temp_string.c_str(),"TOWER_eta_edges")){
	curstring >> varname;	for(unsigned int i=0; i<TOWER_number+1; i++) {curstring >> value; TOWER_eta_edges[i] = value;} }
    else if(strstr(temp_string.c_str(),"TOWER_dphi")){
	curstring >> varname;	for(unsigned int i=0; i<TOWER_number; i++)   {curstring >> value; TOWER_dphi[i] = value;} }

    else if(strstr(temp_string.c_str(),"PTCUT_elec"))       {curstring >> varname >> value; PTCUT_elec        = value;}
    else if(strstr(temp_string.c_str(),"PTCUT_muon"))       {curstring >> varname >> value; PTCUT_muon        = value;}
    else if(strstr(temp_string.c_str(),"PTCUT_jet"))        {curstring >> varname >> value; PTCUT_jet         = value;}
    else if(strstr(temp_string.c_str(),"PTCUT_gamma"))      {curstring >> varname >> value; PTCUT_gamma       = value;}
    else if(strstr(temp_string.c_str(),"PTCUT_taujet"))     {curstring >> varname >> value; PTCUT_taujet      = value;}
    else if(strstr(temp_string.c_str(),"ZDC_gamma_E"))      {curstring >> varname >> value; ZDC_gamma_E      = value;}
    else if(strstr(temp_string.c_str(),"ZDC_n_E"))          {curstring >> varname >> value; ZDC_n_E      = value;}

    else if(strstr(temp_string.c_str(),"ISOL_trk_PT"))      {curstring >> varname >> value; ISOL_trk_PT      = value;}
    else if(strstr(temp_string.c_str(),"ISOL_trk_Cone"))    {curstring >> varname >> value; ISOL_trk_Cone     = value;}
    else if(strstr(temp_string.c_str(),"ISOL_calo_Cone"))    {curstring >> varname >> value; ISOL_calo_Cone     = value;}
    else if(strstr(temp_string.c_str(),"ISOL_calo_ET"))     {curstring >> varname >> value; ISOL_calo_ET      = value;}
    else if(strstr(temp_string.c_str(),"ISOL_calo_Grid"))   {curstring >> varname >> ivalue; ISOL_calo_Grid    = ivalue;}

    else if(strstr(temp_string.c_str(),"JET_coneradius"))   {curstring >> varname >> value; JET_coneradius    = value;}
    else if(strstr(temp_string.c_str(),"JET_jetalgo"))      {curstring >> varname >> ivalue;JET_jetalgo       = ivalue;}
    else if(strstr(temp_string.c_str(),"JET_seed"))         {curstring >> varname >> value; JET_seed          = value;}
    else if(strstr(temp_string.c_str(),"JET_Eflow"))        {curstring >> varname >> ivalue; JET_Eflow        = ivalue;}
 
    else if(strstr(temp_string.c_str(),"BTAG_func_b"))      {curstring >> varname >> svalue; BTAG_func_b     = svalue;}
    else if(strstr(temp_string.c_str(),"BTAG_func_c"))	    {curstring >> varname >> svalue; BTAG_func_c     = svalue;}
    else if(strstr(temp_string.c_str(),"BTAG_func_l"))      {curstring >> varname >> svalue; BTAG_func_l     = svalue;}
    else if(strstr(temp_string.c_str(),"BTAG_func_g"))      {curstring >> varname >> svalue; BTAG_func_g     = svalue;}

    else if(strstr(temp_string.c_str(),"FLAG_vfd"))         {curstring >> varname >> ivalue; FLAG_vfd         = ivalue;}
    else if(strstr(temp_string.c_str(),"FLAG_RP"))          {curstring >> varname >> ivalue; FLAG_RP          = ivalue;}
    else if(strstr(temp_string.c_str(),"FLAG_trigger"))     {curstring >> varname >> ivalue; FLAG_trigger     = ivalue;}
    else if(strstr(temp_string.c_str(),"FLAG_frog"))        {curstring >> varname >> ivalue; FLAG_frog        = ivalue;}
    else if(strstr(temp_string.c_str(),"FLAG_lhco"))        {curstring >> varname >> ivalue; FLAG_lhco        = ivalue;}
    else if(strstr(temp_string.c_str(),"NEvents_Frog"))     {curstring >> varname >> ivalue; NEvents_Frog     = ivalue;}
    else if(strstr(temp_string.c_str(),"NEvents"))          {curstring >> varname >> ivalue; NEvents          = ivalue;}

    else if(strstr(temp_string.c_str(),"PdgTableFilename"))     {curstring >> varname >> svalue; PdgTableFilename = svalue;}
  }

  // for compatibility with old data cards
  if(!CEN_max_calo_ec_flag) {
	cout << "**   Warning \'CEN_max_calo_ec\' not found in datacard.              **"<< endl;
	cout << "**     Same values will be applied for calorimeter endcaps         **"<< endl;
        cout << "**     as for central calorimeters                                 **"<< endl;
        string text = "**     Please update your card ("+ datacard +")";
	cout << left << setw(67) << text << right << setw(2) << "**" << endl;
	cout << "** This change is 100\% backward compatibly with older DetectorCard. **" << endl;
	cout << "** However, please update your DetectorCard. **" << endl;
	CEN_max_calo_ec	 = CEN_max_calo_cen;
	CEN_max_calo_cen = CEN_max_calo_cen/2;	
        ELG_Sec = ELG_Scen;
        ELG_Nec = ELG_Ncen;
        ELG_Cec = ELG_Ccen;
	HAD_Sec = HAD_Scen;
        HAD_Nec = HAD_Ncen;
        HAD_Cec = HAD_Ccen;
  }
 
  if(ISOL_calo_Grid%2 ==0) {
    ISOL_calo_Grid++;
    cout <<"**        WARNING: ISOL_Calo_Grid is not odd. Set it to  "<< ISOL_calo_Grid << "        **" << endl;
  }
 
  //jet stuffs not defined in the input datacard
  JET_overlap     = 0.75;
  // MidPoint algorithm definition
  JET_M_coneareafraction = 0.25;
  JET_M_maxpairsize      = 2;
  JET_M_maxiterations    = 100;
  // Define Cone algorithm.
  JET_C_adjacencycut  =  2;
  JET_C_maxiterations = 100;
  JET_C_iratch        =  1;
  //Define SISCone algorithm.
  JET_S_npass         = 0;
  JET_S_protojet_ptmin= 0.0;
  
  //For Tau-jet definition
  TAU_energy_scone = 0.15;  // radius R of the cone for tau definition, based on energy threshold
  TAU_track_scone  = 0.4;   // radius R of the cone for tau definition, based on track number
  TAU_track_pt     = 2;     // minimal pt [GeV] for tracks to be considered in tau definition
  TAU_energy_frac  = 0.95;  // fraction of energy required in the central part of the cone, for tau jets
  
}

void RESOLution::Logfile(const string& LogName) { 

  // creates the list of good input files
  // this list is vector<string> inputfiles.
  ifstream infile(inputfilelist.c_str());
  vector<string> inputfiles;
  string filename;
  while(1) {
      infile >> filename;                           // reads the first line of the list
      if(!infile.good()) break;                     // quits when at the end of the list
      ifstream checking_the_file(filename.c_str()); // try to open the file
      if(!checking_the_file.good()) continue;       // skips bad/unknown files
      else checking_the_file.close();               // close file if found
      inputfiles.push_back(filename);               // append the name to the vector
  }
  infile.close();

  ofstream f_out(LogName.c_str());

  f_out <<"**********************************************************************"<< endl;
  f_out <<"**********************************************************************"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"**                            Welcome to                            **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"**     .ddddddd-            lL             hH                       **"<< endl;
  f_out <<"**     -Dd` `dD:           Ll             hH`                       **"<< endl;
  f_out <<"**     dDd   dDd    eeee.  lL  .pp+pp    Hh+hhh`   -eeee- `sssss    **"<< endl;
  f_out <<"**    -Dd   `DD   ee. ee  Ll   .Pp. PP   Hh. HH.  ee. ee  sSs       **"<< endl;
  f_out <<"**    dD`   dDd  eEeee:   lL.  pP.  pP  hH   hH` eEeee:` -sSSSs.    **"<< endl;
  f_out <<"**   .Dd  :dd    eE.     LlL   PpppPP   Hh  Hh   eE         sSS     **"<< endl;
  f_out <<"**   dddddd:.     eee+:  lL.  pp.      hh.  hh    eee+  sssssS      **"<< endl;
  f_out <<"**                            Pp                                    **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"**           Delphes, a framework for the fast simulation           **"<< endl;
  f_out <<"**                 of a  generic collider experiment                **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"**                --- Version 1.9 of Delphes ---                    **"<< endl;
  f_out <<"**               Last date of change: 7 May 2010                    **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"**     This package uses:                                           **"<< endl;
  f_out <<"**     ------------------                                           **"<< endl;
  f_out <<"**     FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210]  **"<< endl;
  f_out <<"**     Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]   **"<< endl;
  f_out <<"**     FROG: L. Quertenmont, V. Roberfroid [hep-ex/0901.2718v1]     **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"** ---------------------------------------------------------------- **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"**   Main authors:                                                  **"<< endl;
  f_out <<"**   -------------                                                  **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"**              Séverine Ovyn               Xavier Rouby            **"<< endl;
  f_out <<"**       severine.ovyn@uclouvain.be      xavier.rouby@cern          **"<< endl;
  f_out <<"**       Center for Particle Physics and Phenomenology (CP3)        **"<< endl;
  f_out <<"**       Universite Catholique de Louvain (UCL)                     **"<< endl;
  f_out <<"**       Louvain-la-Neuve, Belgium                                  **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"** ---------------------------------------------------------------- **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"**   Former Delphes versions and documentation can be found on :    **"<< endl;
  f_out <<"**               http://www.fynu.ucl.ac.be/delphes.html             **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"**   Disclaimer: this program comes without guarantees.             **"<< endl;
  f_out <<"**               Beware of errors and please give us                **"<< endl;
  f_out <<"**               your feedbacks about potential bugs.               **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"**********************************************************************"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n"; 
  f_out <<"*                                                                    *"<< endl;
  f_out <<"#********************************                                    *"<<"\n";
  f_out <<"# Input files                                                        *"<<"\n";
  f_out <<"#********************************                                    *"<<"\n";
  f_out << left << setw(22) <<"*  Input list  "<<""
        << left << setw(39) << inputfilelist << "" << right << setw(9) << "*"<<"\n";
  for (unsigned int i =0; i<inputfiles.size(); i++) {
     f_out << left << setw(22) <<"*  - file  "<<""
           << left << setw(43) << inputfiles[i] << "" << right << setw(5) << "*"<<"\n";
  }
  if(detectorcard != "")
     f_out << left << setw(22) <<"*  Detector card   "<<""
           << left << setw(39) << detectorcard << "" << right << setw(9) << "*"<<"\n";
  if(triggercard != "")
     f_out << left << setw(22) <<"*  Trigger card    "<<""
        << left << setw(39) << triggercard << "" << right << setw(9) << "*"<<"\n";
  f_out<<"*  Beam optics :                                                     *"<<"\n";
  f_out << left << setw(22) <<"*  - beam 1   "<<""
        << left << setw(33) << RP_beam1Card << "" << right << setw(15) << "*"<<"\n";
  f_out << left << setw(22) <<"*  - beam 2   "<<""
        << left << setw(33) << RP_beam2Card << "" << right << setw(15) << "*"<<"\n";
  f_out << left << setw(22) <<"*  Input PDG table   " << ""
        << left << setw(39) << PdgTableFilename << "" << right << setw(9) << "*"<<"\n";

  f_out<<"*                                                                    *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out<<"#********************************                                    *"<<"\n";
  f_out<<"# Central detector caracteristics                                    *"<<"\n";
  f_out<<"#********************************                                    *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out << left << setw(30) <<"* Maximum tracking system:                   "<<""
	<< left << setw(10) <<CEN_max_tracker   <<""<< right << setw(15)<<"*"<<"\n";
  f_out << left << setw(30) <<"* Maximum central calorimeter:               "<<""
	<< left << setw(10) <<CEN_max_calo_cen  <<""<< right << setw(15)<<"*"<<"\n";
  f_out << left << setw(30) <<"* Maximum endcap calorimeter:                "<<""
	<< left << setw(10) <<CEN_max_calo_ec  <<""<< right << setw(15)<<"*"<<"\n";
  f_out << left << setw(30) <<"* Maximum central calorimeter:               "<<""
	<< left << setw(10) <<CEN_max_calo_fwd  <<""<< right << setw(15)<<"*"<<"\n";
  f_out << left << setw(30) <<"* Muon chambers coverage:                    "<<""
	<< left << setw(10) <<CEN_max_mu        <<""<< right << setw(15)<<"*"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  if(FLAG_RP==1){
    f_out<<"#************************************                                *"<<"\n";
    f_out<<"# Very forward Roman Pots switched on                                *"<<"\n";
    f_out<<"#************************************                                *"<<"\n";
    f_out<<"*                                                                    *"<<"\n";
    f_out << left << setw(55) <<"* Distance of the 220 RP to the IP in meters:"<<""
	  << left << setw(5) <<RP_220_s             <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* Distance of the 220 RP to the beam in meters:"<<""
	  << left << setw(5) <<RP_220_x             <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* Distance of the 420 RP to the IP in meters:"<<""
	  << left << setw(5) <<RP_420_s             <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* Distance of the 420 RP to the beam in meters:"<<""
	  << left << setw(5) <<RP_420_x             <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* Interaction point at the LHC named: "<<""
          << left << setw(5) <<RP_IP_name             <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(35) <<"* Datacard for beam 1: "<<""
          << left << setw(25) <<RP_beam1Card            <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(35) <<"* Datacard for beam 2: "<<""
          << left << setw(25) <<RP_beam2Card            <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(54) <<"* Beam separation, in meters(hor):"<<""
          << left << setw(6) << RP_offsetEl_x          <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(54) <<"* Beam separation, in meters(ver):"<<""
          << left << setw(6) << RP_offsetEl_y          <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(54) <<"* Distance from IP for Beam separation (m):"<<""
          << left << setw(6) <<RP_offsetEl_s           <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(54) <<"* X offset of beam crossing in micrometers:"<<""
          << left << setw(6) <<RP_cross_x           <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(54) <<"* Y offset of beam crossing in micrometers:"<<""
          << left << setw(6) <<RP_cross_y           <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(54) <<"* X Angle of beam crossing:"<<""
          << left << setw(6) <<RP_cross_ang_x          <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(54) <<"* Y Angle of beam crossing:"<<""
          << left << setw(6) <<RP_cross_ang_y          <<""<< right << setw(10)<<"*"<<"\n";
    f_out<<"*                                                                    *"<<"\n";
  }
  else {
    f_out<<"#*************************************                               *"<<"\n";
    f_out<<"# Very forward Roman Pots switched off                               *"<<"\n";
    f_out<<"#*************************************                               *"<<"\n";
    f_out<<"*                                                                    *"<<"\n";
  }
  if(FLAG_vfd==1){
    f_out<<"#**************************************                              *"<<"\n";
    f_out<<"# Very forward calorimeters switched on                              *"<<"\n";
    f_out<<"#**************************************                              *"<<"\n";
    f_out<<"*                                                                    *"<<"\n";
    f_out << left << setw(55) <<"* Minimum very forward calorimeter:          "<<""
          << left << setw(5) <<VFD_min_calo_vfd <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* Maximum very forward calorimeter:          "<<""
          << left << setw(5) <<VFD_max_calo_vfd <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* Minimum coverage zero_degree calorimeter "<<""
          << left << setw(5) <<VFD_min_zdc          <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* Distance of the ZDC to the IP, in meters:  "<<""
          << left << setw(5) <<VFD_s_zdc            <<""<< right << setw(10)<<"*"<<"\n";
    f_out<<"*                                                                    *"<<"\n";
  }
  else {
    f_out<<"#***************************************                             *"<<"\n";
    f_out<<"# Very forward calorimeters switched off                             *"<<"\n";
    f_out<<"#***************************************                             *"<<"\n";
    f_out<<"*                                                                    *"<<"\n";
  }

  f_out<<"#************************************                                *"<<"\n";
  f_out<<"# Electromagnetic smearing parameters                                *"<<"\n";
  f_out<<"#************************************                                *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  //# \sigma/E = C + N/E + S/\sqrt{E}
  f_out << left << setw(30) <<"* S term for central ECAL: "<<""
	<< left << setw(30) <<ELG_Scen       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* N term for central ECAL: "<<""
	<< left << setw(30) <<ELG_Ncen       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* C term for central ECAL: "<<""
	<< left << setw(30) <<ELG_Ccen       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* S term for ECAL end-cap: "<<""
        << left << setw(30) <<ELG_Sec       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* N term for ECAL end-cap: "<<""
        << left << setw(30) <<ELG_Nec       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* C term for ECAL end-cap: "<<""
        << left << setw(30) <<ELG_Cec       <<""<< right << setw(10)<<"*"<<"\n";


  f_out << left << setw(30) <<"* S term for FCAL: "<<""
	<< left << setw(30) <<ELG_Sfwd       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* N term for FCAL: "<<""
	<< left << setw(30) <<ELG_Nfwd       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* C term for FCAL: "<<""
	<< left << setw(30) <<ELG_Cfwd       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* S term for ZDC: "<<""
        << left << setw(30) <<ELG_Szdc       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* N term for ZDC: "<<""
        << left << setw(30) <<ELG_Nzdc       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* C term for ZDC: "<<""
        << left << setw(30) <<ELG_Czdc       <<""<< right << setw(10)<<"*"<<"\n";

  f_out<<"*                                                                    *"<<"\n";
  f_out<<"#*****************************                                       *"<<"\n";
  f_out<<"# Hadronic smearing parameters                                       *"<<"\n";
  f_out<<"#*****************************                                       *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out << left << setw(30) <<"* S term for central HCAL: "<<""
	<< left << setw(30) <<HAD_Scen       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* N term for central HCAL: "<<""
	<< left << setw(30) <<HAD_Ncen       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* C term for central HCAL: "<<""
	<< left << setw(30) <<HAD_Ccen       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* S term for HCAL endcap:  "<<""
        << left << setw(30) <<HAD_Sec       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* N term for HCAL endcap:  "<<""
        << left << setw(30) <<HAD_Nec       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* C term for HCAL endcap:  "<<""
        << left << setw(30) <<HAD_Cec       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* S term for FCAL: "<<""
	<< left << setw(30) <<HAD_Sfwd        <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* N term for FCAL: "<<""
	<< left << setw(30) <<HAD_Nfwd        <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* C term for FCAL: "<<""
	<< left << setw(30) <<HAD_Cfwd        <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* S term for ZDC: "<<""
        << left << setw(30) <<HAD_Szdc        <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* N term for ZDC: "<<""
        << left << setw(30) <<HAD_Nzdc        <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* C term for ZDC: "<<""
        << left << setw(30) <<HAD_Czdc        <<""<< right << setw(10)<<"*"<<"\n";

  f_out<<"*                                                                    *"<<"\n";
  f_out<<"#*************************                                           *"<<"\n";
  f_out<<"# Time smearing parameters                                           *"<<"\n";
  f_out<<"#*************************                                           *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out << left << setw(55) <<"* Time resolution for ZDC        :              "<<""
        << left << setw(5) <<ZDC_T_resolution       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(55) <<"* Time resolution for RP220      :              "<<""
        << left << setw(5) <<RP220_T_resolution     <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(55) <<"* Time resolution for RP420      :              "<<""
        << left << setw(5) <<RP420_T_resolution     <<""<< right << setw(10)<<"*"<<"\n";
  f_out<<"*                                                                    *"<<"\n";

  f_out<<"*                                                                    *"<<"\n";
  f_out<<"#*************************                                           *"<<"\n";
  f_out<<"# Muon smearing parameters                                           *"<<"\n";
  f_out<<"#*************************                                           *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out << left << setw(55) <<"* PT resolution for muons        :              "<<""
        << left << setw(5) <<MU_SmearPt            <<""<< right << setw(10)<<"*"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  if(FLAG_bfield==1){
    f_out<<"#***************************                                         *"<<"\n";
    f_out<<"# Magnetic field switched on                                         *"<<"\n";
    f_out<<"#***************************                                         *"<<"\n";
    f_out<<"*                                                                    *"<<"\n";
    f_out << left << setw(55) <<"* Radius of the BField coverage:              "<<""
	  << left << setw(5) <<TRACK_radius    <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* Length of the BField coverage:              "<<""
	  << left << setw(5) <<TRACK_length    <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* BField X component:                         "<<""
	  << left << setw(5) <<TRACK_bfield_x  <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* BField Y component:                         "<<""
	  << left << setw(5) <<TRACK_bfield_y  <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* BField Z component:                         "<<""
	  << left << setw(5) <<TRACK_bfield_z  <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* Minimal pT needed to reach the calorimeter [GeV]: "<<""
	  << left << setw(10) <<TRACK_ptmin         <<""<< right << setw(5)<<"*"<<"\n";
    f_out << left << setw(55) <<"* Efficiency associated to the tracking: "<<""
	  << left << setw(10) <<TRACK_eff          <<""<< right << setw(5)<<"*"<<"\n";
    f_out<<"*                                                                    *"<<"\n";
  }
  else {
    f_out<<"#****************************                                        *"<<"\n";
    f_out<<"# Magnetic field switched off                                        *"<<"\n";
    f_out<<"#****************************                                        *"<<"\n";
    f_out << left << setw(55) <<"* Minimal pT needed to reach the calorimeter [GeV]: "<<""
	  << left << setw(10) <<TRACK_ptmin         <<""<< right << setw(5)<<"*"<<"\n";
    f_out << left << setw(55) <<"* Efficiency associated to the tracking: "<<""
	  << left << setw(10) <<TRACK_eff          <<""<< right << setw(5)<<"*"<<"\n";
    f_out<<"*                                                                    *"<<"\n";
  }
  f_out<<"#********************                                                *"<<"\n";
  f_out<<"# Calorimetric Towers                                                *"<<"\n";
  f_out<<"#********************                                                *"<<"\n";
  f_out << left << setw(55) <<"* Number of calorimetric towers in eta, for eta>0: "<<""
        << left << setw(5)  << TOWER_number <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(55) <<"* Tower edges in eta, for eta>0: "<<"" << right << setw(15)<<"*"<<"\n";
  f_out << "*   ";
  for (unsigned int i=0; i<TOWER_number+1; i++) {
    f_out << left << setw(7) << TOWER_eta_edges[i];
    if(!( (i+1) %9 )) f_out << right << setw(3) << "*" << "\n" << "*   ";
  }
  for (unsigned int i=(TOWER_number+1)%9; i<9; i++) f_out << left << setw(7) << "";
  f_out << right << setw(3)<<"*"<<"\n";
  f_out << left << setw(55) <<"* Tower sizes in phi, for eta>0 [degree]:"<<"" << right << setw(15)<<"*"<<"\n";
  f_out << "*   ";
  for (unsigned int i=0; i<TOWER_number; i++) {
    f_out << left << setw(7) << TOWER_dphi[i];
    if(!( (i+1) %9 )) f_out << right << setw(3) << "*" << "\n" << "*   ";
  }
  for (unsigned int i=(TOWER_number)%9; i<9; i++) f_out << left << setw(7) << "";
  f_out << right << setw(3)<<"*"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out<<"#*******************                                                 *"<<"\n";
  f_out<<"# Minimum pT's [GeV]                                                 *"<<"\n";
  f_out<<"#*******************                                                 *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out << left << setw(40) <<"* Minimum pT for electrons: "<<""
	<< left << setw(20) <<PTCUT_elec         <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(40) <<"* Minimum pT for muons: "<<""
	<< left << setw(20) <<PTCUT_muon         <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(40) <<"* Minimum pT for jets: "<<""
	<< left << setw(20) <<PTCUT_jet          <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(40) <<"* Minimum pT for Tau-jets: "<<""
	<< left << setw(20) <<PTCUT_taujet       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(40) <<"* Minimum pT for photons: "<<""
        << left << setw(20) <<PTCUT_gamma        <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(40) <<"* Minimum E for photons in ZDC: "<<""
        << left << setw(20) <<ZDC_gamma_E        <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(40) <<"* Minimum E for neutrons in ZDC: "<<""
        << left << setw(20) <<ZDC_n_E        <<""<< right << setw(10)<<"*"<<"\n";

  f_out<<"*                                                                    *"<<"\n";
  f_out<<"#*******************                                                 *"<<"\n";
  f_out<<"# Isolation criteria                                                 *"<<"\n";
  f_out<<"#*******************                                                 *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out << left << setw(40) <<"* Minimum pT for tracking isolation [GeV]: "<<""
        << left << setw(20) <<ISOL_trk_PT              <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(40) <<"* Cone for tracking isolation criteria: "<<""
        << left << setw(20) <<ISOL_trk_Cone         <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(40) <<"* Cone for calorimetric isolation criteria: "<<""
        << left << setw(20) <<ISOL_calo_Cone         <<""<< right << setw(10)<<"*"<<"\n";


  if(ISOL_calo_ET > 1E98) f_out<<"# No Calorimetric isolation applied                                  *"<<"\n";
  else {
   f_out << left << setw(40) <<"* Minimum ET for towers [GeV]: "<<""
         << left << setw(20) <<ISOL_calo_ET               <<""<< right << setw(10)<<"*"<<"\n";
   f_out << left << setw(40) <<"* Grid size (NxN) for calorimetric isolation: "<<""
         << left << setw(20) <<ISOL_calo_Grid             <<""<< right << setw(4)<<"*"<<"\n";
  }


  f_out<<"*                                                                    *"<<"\n";
  f_out<<"#***************                                                     *"<<"\n";
  f_out<<"# Jet definition                                                     *"<<"\n";
  f_out<<"#***************                                                     *"<<"\n";
  if(JET_Eflow)
    {
      f_out<<"#***************                                                     *"<<"\n";
      f_out<<"#*  Running considering perfect energy flow on the tracker coverage  *"<<"\n";
    }
  else
    {
      f_out<<"#*     Running considering no energy flow on the tracker coverage    *"<<"\n";
      f_out<<"#*         --> jet algo applied on the calorimetric towers           *"<<"\n";
    }
  f_out<<"*                                                                    *"<<"\n";
  f_out<<"*  Six algorithms are currently available:                           *"<<"\n";
  f_out<<"*         - 1) CDF cone algorithm,                                   *"<<"\n";
  f_out<<"*         - 2) CDF MidPoint algorithm,                               *"<<"\n";
  f_out<<"*         - 3) SIScone algorithm,                                    *"<<"\n";
  f_out<<"*         - 4) kt algorithm,                                         *"<<"\n";
  f_out<<"*         - 5) Cambrigde/Aachen algorithm,                           *"<<"\n";
  f_out<<"*         - 6) Anti-kt algorithm.                                    *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out<<"* You have chosen                                                    *"<<"\n";
  switch(JET_jetalgo) {
  default:
  case 1: {
    f_out<<"* CDF JetClu jet algorithm with parameters:                          *"<<"\n";
    f_out << left << setw(40) <<"*        - Seed threshold:    "<<""
	  << left << setw(10) <<JET_seed                        <<""<< right << setw(20)<<"! not in datacard  *"<<"\n";
    f_out << left << setw(40) <<"*        - Cone radius:       "<<""
	  << left << setw(10) <<JET_coneradius                  <<""<< right << setw(20)<<"*"<<"\n";
    f_out << left << setw(40) <<"*        - Adjacency cut:     "<<""
	  << left << setw(10) <<JET_C_adjacencycut              <<""<< right << setw(20)<<"! not in datacard  *"<<"\n";
    f_out << left << setw(40) <<"*        - Max iterations:    "<<""
	  << left << setw(10) <<JET_C_maxiterations             <<""<< right << setw(20)<<"! not in datacard  *"<<"\n";
    f_out << left << setw(40) <<"*        - Iratch:            "<<""
	  << left << setw(10) <<JET_C_iratch                    <<""<< right << setw(20)<<"! not in datacard  *"<<"\n";
    f_out << left << setw(40) <<"*        - Overlap threshold: "<<""
	  << left << setw(10) <<JET_overlap                     <<""<< right << setw(20)<<"! not in datacard  *"<<"\n";
  }
    break;
  case 2: {
    f_out<<"* CDF midpoint jet algorithm with parameters:                        *"<<"\n";
    f_out << left << setw(40) <<"*        - Seed threshold:    "<<""
	  << left << setw(20) <<JET_seed                        <<""<< right << setw(10)<<"! not in datacard  *"<<"\n";
    f_out << left << setw(40) <<"*        - Cone radius:       "<<""
	  << left << setw(20) <<JET_coneradius                  <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(40) <<"*        - Cone area fraction:"<<""
	  << left << setw(20) <<JET_M_coneareafraction          <<""<< right << setw(10)<<"! not in datacard  *"<<"\n";
    f_out << left << setw(40) <<"*        - Maximum pair size: "<<""
	  << left << setw(20) <<JET_M_maxpairsize               <<""<< right << setw(10)<<"! not in datacard  *"<<"\n";
    f_out << left << setw(40) <<"*        - Max iterations:    "<<""
	  << left << setw(20) <<JET_M_maxiterations             <<""<< right << setw(10)<<"! not in datacard  *"<<"\n";
    f_out << left << setw(40) <<"*        - Overlap threshold: "<<""
	  << left << setw(20) <<JET_overlap                     <<""<< right << setw(10)<<"! not in datacard  *"<<"\n";
  }
    break;
  case 3: {
    f_out <<"* SISCone jet algorithm with parameters:                            *"<<"\n";
    f_out << left << setw(40) <<"*        - Cone radius:             "<<""
	  << left << setw(20) <<JET_coneradius                        <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(40) <<"*        - Overlap threshold:       "<<""
	  << left << setw(20) <<JET_overlap                           <<""<< right << setw(10)<<"! not in datacard  *"<<"\n";
    f_out << left << setw(40) <<"*        - Number pass max:         "<<""
	  << left << setw(20) <<JET_S_npass                           <<""<< right << setw(10)<<"! not in datacard  *"<<"\n";
    f_out << left << setw(40) <<"*        - Minimum pT for protojet: "<<""
	  << left << setw(20) <<JET_S_protojet_ptmin                  <<""<< right << setw(10)<<"! not in datacard  *"<<"\n";
  }
    break;
  case 4: {
    f_out <<"* KT jet algorithm with parameters:                                 *"<<"\n";
    f_out << left << setw(40) <<"*        - Cone radius: "<<""
	  << left << setw(20) <<JET_coneradius                <<""<< right << setw(10)<<"*"<<"\n";
  }
    break;
  case 5: {
    f_out <<"* Cambridge/Aachen jet algorithm with parameters:                   *"<<"\n";
    f_out << left << setw(40) <<"*        - Cone radius: "<<""
	  << left << setw(20) <<JET_coneradius                <<""<< right << setw(10)<<"*"<<"\n";
  }
    break;
  case 6: {
    f_out <<"* Anti-kt jet algorithm with parameters:                            *"<<"\n";
    f_out << left << setw(40) <<"*        - Cone radius: "<<""
	  << left << setw(20) <<JET_coneradius                <<""<< right << setw(10)<<"*"<<"\n";
  }
    break;
  }
  f_out<<"*                                                                    *"<<"\n";
  f_out<<"#******************************                                      *"<<"\n";
  f_out<<"# Tau-jet definition parameters                                      *"<<"\n";
  f_out<<"#******************************                                      *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out << left << setw(45) <<"* Cone radius for calorimeter tagging:  "<<""
	<< left << setw(5) <<TAU_energy_scone                           <<""<< right << setw(20)<<"*"<<"\n";
  f_out << left << setw(45) <<"* Fraction of energy in the small cone: "<<""
	<< left << setw(5) <<TAU_energy_frac*100                        <<""<< right << setw(20)<<"! not in datacard  *"<<"\n";
  f_out << left << setw(45) <<"* Cone radius for tracking tagging:     "<<""
	<< left << setw(5) <<TAU_track_scone                            <<""<< right << setw(20)<<"*"<<"\n";
  f_out << left << setw(45) <<"* Minimum track pT [GeV]:               "<<""
	<< left << setw(5) <<TAU_track_pt                                <<""<< right << setw(20)<<"*"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out<<"#*******************************                                     *"<<"\n";
  f_out<<"# B-tagging efficiency functions                                     *"<<"\n";
  f_out<<"#*******************************                                     *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out << left << setw(50) <<"* Tagging a \"b\" as a b-jet: "<<""
        << left << setw(18) <<BTAG_func_b            <<""<< right << setw(2)<<"*"<<"\n";
  f_out << left << setw(50) <<"* Mistagging a c-jet as a b-jet: "<<""
        << left << setw(18) <<BTAG_func_c            <<""<< right << setw(2)<<"*"<<"\n";
  f_out << left << setw(50) <<"* Mistagging a gluon-jet as a b-jet: "<<""
        << left << setw(18) <<BTAG_func_g            <<""<< right << setw(2)<<"*"<<"\n";
  f_out << left << setw(50) <<"* Mistagging a light jet as a b-jet: "<<""
        << left << setw(18) <<BTAG_func_l            <<""<< right << setw(2)<<"*"<<"\n";


  f_out<<"*                                                                    *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out<<"#....................................................................*"<<"\n";
  f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
 
  f_out.close(); 
}

// **********Provides the smeared TLorentzVector for the electrons********
// Smears the electron energy, and changes the 4-momentum accordingly
// different smearing if the electron is central (eta < 2.5) or forward
void RESOLution::SmearElectron(TLorentzVector &electron) {
  // the 'electron' variable will be changed by the function
  float energy = electron.E();  // before smearing
  float energyS = 0.0;		// after smearing  // \sigma/E = C + N/E + S/\sqrt{E}
  
  if(fabs(electron.Eta()) < CEN_max_tracker) { // if the electron is inside the tracker
    energyS = grandom->Gaus(energy, sqrt( 
					 pow(ELG_Ncen,2) + 
					 pow(ELG_Ccen*energy,2) + 
					 pow(ELG_Scen*sqrt(energy),2) ));
  } 
  else if(fabs(electron.Eta()) >= CEN_max_tracker && fabs(electron.Eta()) < CEN_max_calo_ec){
    energyS = grandom->Gaus(energy, sqrt(
                                         pow(ELG_Nec,2) +
                                         pow(ELG_Cec*energy,2) +
                                         pow(ELG_Sec*sqrt(energy),2) ) );
  }
  else if(fabs(electron.Eta()) >= CEN_max_calo_ec && fabs(electron.Eta()) < CEN_max_calo_fwd){
    energyS = grandom->Gaus(energy, sqrt(
					 pow(ELG_Nfwd,2) +
					 pow(ELG_Cfwd*energy,2) +
					 pow(ELG_Sfwd*sqrt(energy),2) ) ); 
  }
  electron.SetPtEtaPhiE(energyS/cosh(electron.Eta()), electron.Eta(), electron.Phi(), energyS);
  if(electron.E() < 0)electron.SetPxPyPzE(0,0,0,0); // no negative values after smearing !
}


// **********Provides the smeared TLorentzVector for the muons********
//  Smears the muon pT and changes the 4-momentum accordingly
void RESOLution::SmearMu(TLorentzVector &muon) {
  // the 'muon' variable will be changed by the function
  float pt = muon.Pt(); // before smearing
  float ptS=pt;

  if(fabs(muon.Eta()) < CEN_max_mu )
   {
     ptS = grandom->Gaus(pt, MU_SmearPt*pt ); // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
   }
     muon.SetPtEtaPhiE(ptS, muon.Eta(), muon.Phi(), ptS*cosh(muon.Eta()));
  
  if(muon.E() < 0)muon.SetPxPyPzE(0,0,0,0); // no negative values after smearing !
}


// **********Provides the smeared TLorentzVector for the hadrons********
//  Smears the hadron 4-momentum
void RESOLution::SmearHadron(TLorentzVector &hadron, const float frac)
  // the 'hadron' variable will be changed by the function
  // the 'frac' variable describes the long-living particles. Should be 0.7 for K0S and Lambda, 1. otherwise
{
  float energy = hadron.E();  // before smearing
  float energyS = 0.0;        // after smearing  // \sigma/E = C + N/E + S/\sqrt{E}
  float energy_ecal = (1.0 - frac)*energy; // electromagnetic calorimeter
  float energy_hcal = frac*energy;	   // hadronic calorimeter
  // frac takes into account the decay of long-living particles, that decay in the calorimeters
  // some of the particles decay mostly in the ecal, some mostly in the hcal
  
  float energyS1,energyS2;
  if(fabs(hadron.Eta()) < CEN_max_calo_cen) {
   energyS1 = grandom->Gaus(energy_hcal, sqrt(
					      pow(HAD_Ncen,2) +
					      pow(HAD_Ccen*energy_hcal,2) + 
					      pow(HAD_Scen*sqrt(energy_hcal),2) )) ;
      
   energyS2 =  grandom->Gaus(energy_ecal, sqrt(
				      pow(ELG_Ncen,2) +
				      pow(ELG_Ccen*energy_ecal,2) + 
				      pow(ELG_Scen*sqrt(energy_ecal),2) ) );     

  energyS = ((energyS1>0)?energyS1:0) + ((energyS2>0)?energyS2:0);
  } 
  else if(fabs(hadron.Eta()) >= CEN_max_calo_cen && fabs(hadron.Eta()) < CEN_max_calo_ec) {
   energyS1 = grandom->Gaus(energy_hcal, sqrt(
                                              pow(HAD_Nec,2) +
                                              pow(HAD_Cec*energy_hcal,2) +
                                              pow(HAD_Sec*sqrt(energy_hcal),2) )) ;

   energyS2 =  grandom->Gaus(energy_ecal, sqrt(
                                      pow(ELG_Nec,2) +
                                      pow(ELG_Cec*energy_ecal,2) +
                                      pow(ELG_Sec*sqrt(energy_ecal),2) ) );

  energyS = ((energyS1>0)?energyS1:0) + ((energyS2>0)?energyS2:0);
  }
  else if(fabs(hadron.Eta()) >= CEN_max_calo_ec && fabs(hadron.Eta()) < CEN_max_calo_fwd){ 
    energyS = grandom->Gaus(energy, sqrt(
			    pow(HAD_Nfwd,2) +
			    pow(HAD_Cfwd*energy,2) +
			    pow(HAD_Sfwd*sqrt(energy),2) ));
}
  


  hadron.SetPtEtaPhiE(energyS/cosh(hadron.Eta()),hadron.Eta(), hadron.Phi(), energyS);
  
  if(hadron.E() < 0)hadron.SetPxPyPzE(0,0,0,0);
}

//******************************************************************************************

//void RESOLution::SortedVector(vector<ParticleUtil> &vect)
void RESOLution::SortedVector(vector<D_Particle> &vect)
{
  int i,j = 0;
  TLorentzVector tmp;
  bool en_desordre = true;
  int entries=vect.size();
  for(i = 0 ; (i < entries) && en_desordre; i++)
    {
      en_desordre = false;
      for(j = 1 ; j < entries - i ; j++)
        {
          if ( vect[j].Pt() > vect[j-1].Pt() )
             {
               //ParticleUtil tmp = vect[j-1];
               D_Particle tmp = vect[j-1];
               vect[j-1] = vect[j];
               vect[j] = tmp;
               en_desordre = true;
              }
        }
    }
}

// **********Provides the energy in the cone of radius TAU_CONE_ENERGY for the tau identification********
// to be taken into account, a calo tower should
// 	1) have a transverse energy \f$ E_T = \sqrt{E_X^2 + E_Y^2} \f$ above a given threshold
// 	2) be inside a cone with a radius R and the axis defined by (eta,phi)
double RESOLution::EnergySmallCone(const vector<PhysicsTower> &towers, const float eta, const float phi) {
  double Energie=0;
  for(unsigned int i=0; i < towers.size(); i++) {
      if(towers[i].fourVector.pt() < JET_seed) continue;
      if((DeltaR(phi,eta,towers[i].fourVector.phi(),towers[i].fourVector.eta()) < TAU_energy_scone)) {
          Energie += towers[i].fourVector.E;
      }
  }
  return Energie;
}


// **********Provides the number of tracks in the cone of radius TAU_CONE_TRACKS for the tau identification********
// to be taken into account, a track should
// 	1) avec a transverse momentum \$f p_T \$ above a given threshold
// 	2) be inside a cone with a radius R and the axis defined by (eta,phi)
// IMPORTANT REMARK !!!!!
// NEW : "charge" will contain the sum of all charged tracks in the cone TAU_track_scone
unsigned int RESOLution::NumTracks(float& charge, const vector<TRootTracks> &tracks, const float pt_track, const float eta, const float phi) {
  unsigned int numbtrack=0; // number of track in the tau-jet cone, which is smaller than R;
  charge=0;
  for(unsigned int i=0; i < tracks.size(); i++) { 
      if(tracks[i].PT < pt_track ) continue;
      //float dr = DeltaR(phi,eta,tracks[i].PhiOuter,tracks[i].EtaOuter);
      float dr = DeltaR(phi,eta,tracks[i].Phi,tracks[i].Eta);
      if (dr > TAU_track_scone) continue; 
	numbtrack++;
      	charge += tracks[i].Charge; // total charge in the cone for Tau-jet
  }
  return numbtrack;
}

//*** Returns the PID of the particle with the highest energy, in a cone with a radius CONERADIUS and an axis (eta,phi)  *********
//used by Btaggedjet
///// Attention : bug removed => CONERADIUS/2 -> CONERADIUS !!
int RESOLution::Bjets(const TSimpleArray<TRootC::GenParticle> &subarray, const float& eta, const float& phi) {
  float emax=0;
  int Ppid=0;
  if(subarray.GetEntries()>0) {
      for(int i=0; i < subarray.GetEntries();i++) { // should have pt>PT_JETMIN and a small cone radius (r<CONE_JET)
	  float genDeltaR = DeltaR(subarray[i]->Phi,subarray[i]->Eta,phi,eta);
          if(genDeltaR < JET_coneradius && subarray[i]->E > emax) {
              emax=subarray[i]->E;
              Ppid=abs(subarray[i]->PID);
          }
      }
  }
  return Ppid; 
}


//******************* Simulates the b-tagging efficiency for real bjet, or the misendentification for other jets****************

bool RESOLution::Btaggedjet(const TLorentzVector &JET, const TSimpleArray<TRootC::GenParticle> &subarray) {

  const float ptmax = 9E6; // GeV/c
  bool tag = false;
  string efficiency_formula="";

  // selects the correct formula
  int jet_id = Bjets(subarray,JET.Eta(),JET.Phi()); // gets the particle id of the most energetic parton in the jet
  switch (jet_id) {
    case pB: efficiency_formula = BTAG_func_b; break;
    case pC: efficiency_formula = BTAG_func_c; break;
    case pGLUON: efficiency_formula = BTAG_func_g; break;
    default: efficiency_formula = BTAG_func_l; break;
  }

  TF2 efficiency("efficiency",efficiency_formula.c_str(),0,ptmax,-CEN_max_tracker,CEN_max_tracker); // efficiency function wrt Pt
  float x = grandom->Uniform(); // randomisation
  tag = (x < efficiency.Eval(JET.Pt(),JET.Eta())) ? true : false;
  //cout << "formula = " << efficiency_formula << "\t Pt = " << JET.Pt() 
  //     << "\t value = " << efficiency.Eval(JET.Pt())  
  //     << "\t x = " << x << "\t tag= " << tag << endl;

  return tag;
}


//***********************Isolation criteria***********************
//****************************************************************
bool RESOLution::Isolation(const D_Particle& part, const vector<TRootTracks> &tracks, const float& pt_second_track, float& ptiso )
{
   bool isolated = false;
   ptiso = 0; // sum of all track pt in isolation cone
   float deltar=1E99; // Initial value; should be high; no further repercussion

   // loop on all tracks, with p_t above threshold, close enough from the charged lepton
   for(unsigned int i=0; i < tracks.size(); i++) {
       if(tracks[i].PT < pt_second_track) continue; // ptcut on tracks
       float genDeltaR = DeltaR(part.Phi(),part.Eta(),tracks[i].Phi,tracks[i].Eta); 
         if(
             (genDeltaR==0)  // rejets the track of the particle itself
           ) continue ;
           if (genDeltaR < deltar) deltar=genDeltaR;	// finds the closest track

         // as long as (genDeltaR==0) is put above, the particle itself is not taken into account
	 if( genDeltaR < ISOL_trk_Cone && tracks[i].PT>ISOL_trk_PT) ptiso += tracks[i].PT; // dR cut on tracks
      }
   if(deltar > ISOL_trk_Cone) isolated = true;
   return isolated;
}

float RESOLution::IsolationSumEt(const float iPhi, const float iEta, D_CaloTowerList & towers)
{
   float etiso = 0; // sum of all track pt in isolation cone
   for (unsigned int i=0; i< towers.size(); i++) {
        if( (DeltaR(iPhi,iEta,towers[i].getPhi(),towers[i].getEta()) < ISOL_calo_Cone) && (towers[i].getET() > ISOL_calo_ET) ) {
           etiso += towers[i].getET(); // dR cut on tracks
          }
  }
   return etiso;
}



// ******* Calorimetric isolation -- for LHCO only
float RESOLution::CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers, const float iPhi, const float iEta, const int iTower) {
  // iPhi and iEta are the coordinates of the center of the tower crossed by the particle
  // iPhi is in [-pi ; pi]    and   iEta is in [-etamax ; etamax]
  // iTower is the index of the tower, in [0, n_tower]. iTower points only towers in positive range

   float et_sum=0;

   unsigned int N = ISOL_calo_Grid;
   if(N%2 ==0 || N==0) { cout << "Warning: ISOL_Calo_Grid must be odd. ISOL_Calo_Grid = 3 will be used\n"; N=3;}
   int index= iTower; // index of the central tower of the grid in TOWER_eta_edges[.]; 
                      // !! TOWER_eta_edges is only with eta>0

   //cout << "iEta= " << iEta << "\tiPhi= " << iPhi << endl;
   if(index != iUNDEFINED) {
      int range = (int) (N-1)/2;           // the (int) conversion is needed, as N is "unsigned int"
      for(int i= -range; i<= range; i++) { // shift of index in eta: e.g.: i = -1, 0, 1 if N=3;
      // goal : to identify the center of each cell in the NxN grid.
      // the eta/phi coordinates of each center will be in (eta_ith_tower,phi_ith_tower)
         unsigned int i_eta = (index+i>=0) ? index + i : (int) -1*(index + i +1); // i_eta is the index in [0; index_eta_max]
         if(i_eta>=TOWER_number) continue; // avoid going too far: TOWER_dphi[TOWER_number], TOWER_eta_edges[TOWER_number+1]
         //cout << "i_eta" << i_eta << "\t TOWER_number = " << TOWER_number << endl;

         float dphi = TOWER_dphi[i_eta]*pi/180.; // in rad // size in phi of the cells for this eta
         float eta_ith_tower = (TOWER_eta_edges[i_eta]+TOWER_eta_edges[i_eta+1])/2.0;
         if(iEta<0) eta_ith_tower *= -1;     // needed if the central tower is in negative region
         if(index+i<0) eta_ith_tower *= -1;  // needed if the "eta=0"-axis is crossed by the grid
         //cout << "pour eta_ith_tower=" << eta_ith_tower << ", on va dphi = " << dphi << endl;
         // !!! at this point, eta_ith_tower is the value in eta of the center of the current calo cell

         for(int j= -range; j<= range; j++) { // shift of index in phi. e.g.: j = -1, 0, 1 if N=3;
            float phi_ith_tower = iPhi + j*dphi; // in radian
            //cout << "eta_ith_tower= "   << eta_ith_tower     << "\tphi_ith_tower= " << phi_ith_tower << "\tj= " << j << "\tdphi=" << dphi << endl;
            if(phi_ith_tower > pi) phi_ith_tower -= 2.*pi;
            else if (phi_ith_tower < -pi) phi_ith_tower += 2.*pi;
            
            D_CaloTower calMuon(towers.getElement(eta_ith_tower,phi_ith_tower));
            //cout << "eta_ith_tower= "   << eta_ith_tower     << "\tphi_ith_tower= " << phi_ith_tower <<"\t"
            //     << "calMuon.getEta= " << calMuon.getEta() << "\tcalMuon.getPhi()= " << calMuon.getPhi() <<"\t";

            //if(calMuon.getEta() != UNDEFINED && calMuon.getET() > ISOL_calo_ET) {
            if(calMuon.getEta() != UNDEFINED && calMuon.getET() > 0.0) {
               et_sum += calMuon.getET();
               //cout << calMuon.getET() << " GeV";
            }
            //else cout << " not active";
         } // j-loop (phi indices)
      } // i-loop (eta indices)
   }  // undefined index
   else {
     if (CEN_max_mu < CEN_max_calo_fwd) 
        cout << "**    ERROR in RESOLution::CaloIsolation: 'muon'-tower not found!  **" << endl; 
   }      // should never happen ! this would be a bug
  return (part.Pt()==0)? UNDEFINED: et_sum/part.Pt();
}


  //********** returns a segmented value for eta and phi, for calo towers *****
int RESOLution::BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta){
   iEta = UNDEFINED;
   int index= iUNDEFINED;
   for (unsigned int i=1; i< TOWER_number+1; i++) {
        if(fabs(eta)>TOWER_eta_edges[i-1] && fabs(eta)<=TOWER_eta_edges[i]) {
                iEta = (eta>0) ? ((TOWER_eta_edges[i-1]+TOWER_eta_edges[i])/2.0) : -((TOWER_eta_edges[i-1]+TOWER_eta_edges[i])/2.0);
                index = i-1;
                break;
        }
   }
   iPhi = UNDEFINED;
   if(index==iUNDEFINED) return index; 

   float dphi = TOWER_dphi[index]*pi/180.; 
   for (unsigned int i=1; i < 360/TOWER_dphi[index]; i++ ) {
        float low = -pi+(i-1)*dphi;
        float high= low+dphi;
        if(phi > low && phi <= high ){
                iPhi = (low+high)/2.0; 
                break;
        } 
   }
   if (phi > pi-dphi) iPhi = pi-dphi;
   return index;
}

//**************************** Returns the delta Phi ****************************
float DeltaPhi(const float phi1, const float phi2) {
  float deltaphi=phi1-phi2; // in here, -pi < phi < pi
  if(fabs(deltaphi) > pi) {
	deltaphi=2.*pi -fabs(deltaphi);// put deltaphi between 0 and pi
	}
  else deltaphi=fabs(deltaphi);
  
  return deltaphi;
}

//**************************** Returns the delta R****************************
float DeltaR(const float phi1, const float eta1, const float phi2, const float eta2) {
  return sqrt(pow(DeltaPhi(phi1,phi2),2) + pow(eta1-eta2,2));
}

int sign(const int myint) { 
	if (myint >0) return 1; 
	else if (myint <0) return -1; 
	else return 0;
}

int sign(const float myfloat) {
        if (myfloat >0) return 1;
        else if (myfloat <0) return -1;
        else return 0;
}

int ChargeVal(const int pid)
{
  cout << "ChargeVal :: deprecated function, do not use it anymore" << endl;
  int charge;
  if(
     (pid == pGAMMA)    ||
     (pid == pPI0)      ||
     (pid == pK0L)      ||
     (pid == pN)        ||
     (pid == pSIGMA0)   ||
     (pid == pDELTA0)   ||
     (pid == pK0S)   // not charged particles : invisible by tracker
     )
    charge = 0;
  else charge = sign(pid);
  return charge;
            
}

//------------------------------------------------------------------------------
void RESOLution::ReadParticleDataGroupTable() {

  string temp_string;
  istringstream curstring;

  ifstream fichier_a_lire(PdgTableFilename.c_str());
  if(!fichier_a_lire.good()) {
    cout <<"**        ERROR: PDG Table ("<< PdgTableFilename 
         <<  ") not found! exit.                        **" << endl;
    exit(1);
    return;
  }
  // first three lines of the file are useless 
  getline(fichier_a_lire,temp_string);
  getline(fichier_a_lire,temp_string);
  getline(fichier_a_lire,temp_string);


  while (getline(fichier_a_lire,temp_string)) {
    curstring.clear(); // needed when using several times istringstream::str(string)
    curstring.str(temp_string);
    long int ID; std::string name; int charge; float mass; float width; float lifetime;
    // ID name   chg       mass    total width   lifetime
    //  1 d      -1      0.33000     0.00000   0.00000E+00
    //  in the table, the charge is in units of e+/3
    //  the total width is in GeV
    //  the lifetime is ctau in mm
    curstring >> ID >> name >> charge >> mass >> width >> lifetime;
    PdgParticle particle(ID,name,mass,charge/3.,width,lifetime/1000.);
    PDGtable.insert(ID,particle);
    //PdgTable.insert(pair<int,PdgParticle>(ID,particle));
    //cout << PDGtable[ID].name() <<  "\t" << PDGtable[ID].mass() << "\t" << PDGtable[ID].charge() << endl;
  }
  
} // ReadParticleDataGroupTable


// to be improved in order to avoid code repetition
// sorry, no time to do it right now (XR, 19/05/2009)
void time_report(const TStopwatch& global,const TStopwatch& loop,const TStopwatch& trigger,const TStopwatch& frog,const TStopwatch& lhco, const int flag_frog, const int flag_trigger, const int flag_lhco, const string& LogName, const Long64_t allEntries) {

TStopwatch globalwatch(global), loopwatch(loop), triggerwatch(trigger), frogwatch(frog), lhcowatch(lhco);

  cout <<"**                                                                 **"<< endl;
  cout <<"**        ################## Time report #################         **"<< endl;
  cout << left  << setw(32) <<"**              Time report for "<<""
       << left  << setw(15) << allEntries <<""
       << right << setw(22) <<"events         **"<<endl;
  cout <<"**                                                                 **"<< endl;
  cout << left  << setw(10) <<"**"<<""
       << left  << setw(15) <<"Time (s):"<<""
       << right << setw(15) <<"CPU"<<""
       << right << setw(15) <<"Real"<<""
       << right << setw(14)  <<"**"<<endl;
  cout << left  << setw(10) <<"**"<<""
       << left  << setw(15) <<" +  Global:"<<""
       << right << setw(15) <<globalwatch.CpuTime()<<""
       << right << setw(15) <<globalwatch.RealTime()<<""
       << right << setw(14)  <<"**"<<endl;
  cout << left  << setw(10) <<"**"<<""
       << left  << setw(15) <<" +  Events:"<<""
       << right << setw(15) <<loopwatch.CpuTime()<<""
       << right << setw(15) <<loopwatch.RealTime()<<""
       << right << setw(14)  <<"**"<<endl;
  if(flag_trigger == 1)
    {
      cout << left  << setw(10) <<"**"<<""
	   << left  << setw(15) <<" +  Trigger:"<<""
	   << right << setw(15) <<triggerwatch.CpuTime()<<""
	   << right << setw(15) <<triggerwatch.RealTime()<<""
	   << right << setw(14)  <<"**"<<endl;
    }
  if(flag_frog == 1)
    {
      cout << left  << setw(10) <<"**"<<""
	   << left  << setw(15) <<" +  Frog:"<<""
	   << right << setw(15) <<frogwatch.CpuTime()<<""
	   << right << setw(15) <<frogwatch.RealTime()<<""
	   << right << setw(14)  <<"**"<<endl;
    }
  if(flag_lhco == 1)
    {
      cout << left  << setw(10) <<"**"<<""
           << left  << setw(15) <<" +  LHCO:"<<""
           << right << setw(15) <<lhcowatch.CpuTime()<<""
           << right << setw(15) <<lhcowatch.RealTime()<<""
           << right << setw(14)  <<"**"<<endl;
    }


  ofstream f_out(LogName.c_str(),ios_base::app);

  f_out <<"**                                                                   *"<< endl;
  f_out <<"**        ################## Time report #################           *"<< endl;
  f_out << left  << setw(32) <<"**              Time report for "<<""
       << left  << setw(15) << allEntries <<""
       << right << setw(23) <<"events             *"<<endl;
  f_out <<"**                                                                   *"<< endl;
  f_out << left  << setw(10) <<"**"<<""
       << left  << setw(15) <<"Time (s):"<<""
       << right << setw(15) <<"CPU"<<""
       << right << setw(15) <<"Real"<<""
       << right << setw(15)  <<" *"<<endl;
  f_out << left  << setw(10) <<"**"<<""
       << left  << setw(15) <<" +  Global:"<<""
       << right << setw(15) <<globalwatch.CpuTime()<<""
       << right << setw(15) <<globalwatch.RealTime()<<""
       << right << setw(15)  <<" *"<<endl;
  f_out << left  << setw(10) <<"**"<<""
       << left  << setw(15) <<" +  Events:"<<""
       << right << setw(15) <<loopwatch.CpuTime()<<""
       << right << setw(15) <<loopwatch.RealTime()<<""
       << right << setw(15)  <<" *"<<endl;
  if(flag_trigger == 1)
    {
      f_out << left  << setw(10) <<"**"<<""
           << left  << setw(15) <<" +  Trigger:"<<""
           << right << setw(15) <<triggerwatch.CpuTime()<<""
           << right << setw(15) <<triggerwatch.RealTime()<<""
           << right << setw(15)  <<" *"<<endl;
    }
  if(flag_frog == 1)
    {
      f_out << left  << setw(10) <<"**"<<""
           << left  << setw(15) <<" +  Frog:"<<""
           << right << setw(15) <<frogwatch.CpuTime()<<""
           << right << setw(15) <<frogwatch.RealTime()<<""
           << right << setw(15)  <<" *"<<endl;
    }
  if(flag_lhco == 1)
    {
      f_out << left  << setw(10) <<"**"<<""
           << left  << setw(15) <<" +  LHCO:"<<""
           << right << setw(15) <<lhcowatch.CpuTime()<<""
           << right << setw(15) <<lhcowatch.RealTime()<<""
           << right << setw(15)  <<" *"<<endl;
    }
  f_out << left << setw(16) << "**        " << ""
       << left << setw(52) <<  get_time_date() << "**" << endl;


  f_out<<"*                                                                    *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out<<"#....................................................................*"<<"\n";
  f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";

  f_out.close(); 

}

void print_header() {
  cout << endl << endl;
  
  cout <<"*********************************************************************"<< endl;
  cout <<"*********************************************************************"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**                            Welcome to                           **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**     .ddddddd-            lL             hH                      **"<< endl;
  cout <<"**     -Dd` `dD:           Ll             hH`                      **"<< endl;
  cout <<"**     dDd   dDd    eeee.  lL  .pp+pp    Hh+hhh`   -eeee- `sssss   **"<< endl;
  cout <<"**    -Dd   `DD   ee. ee  Ll   .Pp. PP   Hh. HH.  ee. ee  sSs      **"<< endl;
  cout <<"**    dD`   dDd  eEeee:   lL.  pP.  pP  hH   hH` eEeee:` -sSSSs.   **"<< endl;
  cout <<"**   .Dd  :dd    eE.     LlL   PpppPP   Hh  Hh   eE         sSS    **"<< endl;
  cout <<"**   dddddd:.     eee+:  lL.  pp.      hh.  hh    eee+  sssssS     **"<< endl;
  cout <<"**                            Pp                                   **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**           Delphes, a framework for the fast simulation          **"<< endl;
  cout <<"**                 of a  generic collider experiment               **"<< endl;
  cout <<"**                    arXiv:0903.2225v1 [hep-ph]                   **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**                  --- Version 1.9 of Delphes ---                 **"<< endl;
  cout <<"**                Last date of change:   7 May 2010                **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**     This package uses:                                          **"<< endl;
  cout <<"**     ------------------                                          **"<< endl;
  cout <<"**     FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **"<< endl;
  cout <<"**     Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]  **"<< endl;
  cout <<"**     FROG:                               [hep-ex/0901.2718v1]    **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**-----------------------------------------------------------------**"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**   Main authors:                                                 **"<< endl;
  cout <<"**   -------------                                                 **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**              Séverine Ovyn               Xavier Rouby           **"<< endl;
  cout <<"**       severine.ovyn@uclouvain.be      xavier.rouby@cern         **"<< endl;
  cout <<"**       Center for Particle Physics and Phenomenology (CP3)       **"<< endl;
  cout <<"**       Universite Catholique de Louvain (UCL)                    **"<< endl;
  cout <<"**       Louvain-la-Neuve, Belgium                                 **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**-----------------------------------------------------------------**"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**   Former Delphes versions and documentation can be found on :   **"<< endl;
  cout <<"**               http://www.fynu.ucl.ac.be/delphes.html            **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**   Disclaimer: this program comes without guarantees.            **"<< endl;
  cout <<"**               Beware of errors and please give us               **"<< endl;
  cout <<"**               your feedbacks about potential bugs.              **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"*********************************************************************"<< endl;
  cout <<"*********************************************************************"<< endl;

}

string get_time_date() {
  time_t rawtime;
  struct tm * timeinfo;

  time ( &rawtime );
  timeinfo = localtime ( &rawtime );
  
  char temp[100];
  //sprintf(temp,"%i/%i/%i %i:%i:%i",timeinfo->tm_mday,timeinfo->tm_mon+1,timeinfo->tm_year+1900,timeinfo->tm_hour,timeinfo->tm_min,timeinfo->tm_sec);
  sprintf(temp,"%i",timeinfo->tm_mday);
  if(timeinfo->tm_mon+1<10) sprintf(temp,"%s/0%d",temp,timeinfo->tm_mon+1);else sprintf(temp,"%s/%d",temp,timeinfo->tm_mon+1);
  sprintf(temp,"%s/%d %d",temp,timeinfo->tm_year+1900,timeinfo->tm_hour);
  if(timeinfo->tm_min<10)   sprintf(temp,"%s:0%d",temp,timeinfo->tm_min);  else sprintf(temp,"%s:%d",temp,timeinfo->tm_min);
  if(timeinfo->tm_sec<10)   sprintf(temp,"%s:0%d",temp,timeinfo->tm_sec);  else sprintf(temp,"%s:%d",temp,timeinfo->tm_sec);
  string tempstring(temp);
  return tempstring;
}

void warning(const string oldname, const string newname) {
  string text = "** Warning in datacard: " + oldname + " deprecated. Use " + newname + " instead.";
  cout << left  << setw(67) << text 
       << right << setw(2)  <<"**" << endl;
}