source: svn/trunk/src/SmearUtil.cc@ 452

/***********************************************************************
**                                                                    **
**          /----------------------------------------------\          **
**         |  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-2009,                      **
**                        All rights reserved.                        **
**                                                                    **
***********************************************************************/

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


#include "SmearUtil.h"
#include "TRandom.h"
#include "TStopwatch.h"

#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <map>
using namespace std;

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

RESOLution::RESOLution() {

  // Detector characteristics
  CEN_max_tracker   = 2.5;                // Maximum tracker coverage          
  CEN_max_calo_cen  = 3.0;                // central calorimeter 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 for central ECAL
  ELG_Ccen          = 0.005;             // C term for central ECAL
  ELG_Sfwd          = 2.084;             // S term for FCAL
  ELG_Nfwd          = 0.0;               // N term for FCAL
  ELG_Cfwd          = 0.107;             // C term for FCAL
  ELG_Szdc          = 0.70;              // S term for ZDC
  ELG_Nzdc          = 0.0;               // N term for ZDC
  ELG_Czdc          = 0.08;              // C term for ZDC 

  // Energy resolution for hadrons in ecal/hcal/hf
  // \sigma/E = C + N/E + S/\sqrt{E}
  HAD_Shcal         = 1.5;               // S term for central HCAL
  HAD_Nhcal         = 0.;                // N term for central HCAL
  HAD_Chcal         = 0.05;              // C term for central HCAL
  HAD_Shf           = 2.7;               // S term for FCAL 
  HAD_Nhf           = 0.;                // N term for FCAL
  HAD_Chf           = 0.13;              // C term for FCAL
  HAD_Szdc          = 1.38;              // S term for ZDC
  HAD_Nzdc          = 0.;                // N term for ZDC 
  HAD_Czdc          = 0.13;              // C term for ZDC 

  // 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_PT         = 2.0;      //minimal pt of tracks for isolation criteria
  ISOL_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      = 1;              // 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        = 1;              // 1 for Energy flow in jets reco ; 0 if not

  // Tagging definition
  BTAG_b           = 40;
  BTAG_mistag_c    = 10;
  BTAG_mistag_l    = 1;

  // 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";
  
}


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_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_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/hf/zdc
  HAD_Shcal         = DET.HAD_Shcal;
  HAD_Nhcal         = DET.HAD_Nhcal;
  HAD_Chcal         = DET.HAD_Chcal;
  HAD_Shf           = DET.HAD_Shf;
  HAD_Nhf           = DET.HAD_Nhf;
  HAD_Chf           = DET.HAD_Chf;
  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_PT         = DET.ISOL_PT;     // tracking isolation
  ISOL_Cone       = DET.ISOL_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_b           = DET.BTAG_b;
  BTAG_mistag_c    = DET.BTAG_mistag_c;
  BTAG_mistag_l    = DET.BTAG_mistag_l;

  // 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;
}

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_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_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/hf
  HAD_Shcal         = DET.HAD_Shcal;
  HAD_Nhcal         = DET.HAD_Nhcal;
  HAD_Chcal         = DET.HAD_Chcal;
  HAD_Shf           = DET.HAD_Shf;
  HAD_Nhf           = DET.HAD_Nhf;
  HAD_Chf           = DET.HAD_Chf;
  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_PT         = DET.ISOL_PT;       // tracking isolation
  ISOL_Cone       = DET.ISOL_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_b           = DET.BTAG_b;
  BTAG_mistag_c    = DET.BTAG_mistag_c;
  BTAG_mistag_l    = DET.BTAG_mistag_l;

  // 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;
  return *this;
}




//------------------------------------------------------------------------------
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;
  }
  
  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_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_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"))        {curstring >> varname >> value; HAD_Shcal         = value;}
    else if(strstr(temp_string.c_str(),"HAD_Nhcal"))        {curstring >> varname >> value; HAD_Nhcal         = value;}
    else if(strstr(temp_string.c_str(),"HAD_Chcal"))        {curstring >> varname >> value; HAD_Chcal         = value;}
    else if(strstr(temp_string.c_str(),"HAD_Shf"))          {curstring >> varname >> value; HAD_Shf           = value;}
    else if(strstr(temp_string.c_str(),"HAD_Nhf"))          {curstring >> varname >> value; HAD_Nhf           = value;}
    else if(strstr(temp_string.c_str(),"HAD_Chf"))          {curstring >> varname >> value; HAD_Chf           = 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 >> ivalue;TRACK_eff         = ivalue;}

    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_PT"))          {curstring >> varname >> value; ISOL_PT           = value;}
    else if(strstr(temp_string.c_str(),"ISOL_Cone"))        {curstring >> varname >> value; ISOL_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_b"))           {curstring >> varname >> ivalue;BTAG_b            = ivalue;}
    else if(strstr(temp_string.c_str(),"BTAG_mistag_c"))    {curstring >> varname >> ivalue;BTAG_mistag_c     = ivalue;}
    else if(strstr(temp_string.c_str(),"BTAG_mistag_l"))    {curstring >> varname >> ivalue;BTAG_mistag_l     = ivalue;}

    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;}
  }
 
  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) { 
  //void RESOLution::Logfile(string outputfilename) { 
  
  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.7 of Delphes ---                    **"<< endl;
  f_out <<"**               Last date of change: 7 May 2009                    **"<< 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 is a beta version of Delphes and      **"<< endl;
  f_out <<"** therefore comes without guarantees. Beware of errors and please  **"<< endl;
  f_out <<"**        give us your feedbacks about potential bugs               **"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out <<"**********************************************************************"<< endl;
  f_out <<"**                                                                  **"<< endl;
  f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n"; 
  f_out<<"#  Input PDG table : " << PdgTableFilename << "                               *"<<"\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 forward 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(44) <<"* Beam separation, in meters(hor):"<<""
          << left << setw(6) << RP_offsetEl_x          <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(44) <<"* Beam separation, in meters(ver):"<<""
          << left << setw(6) << RP_offsetEl_y          <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(44) <<"* Distance from IP for Beam separation (m):"<<""
          << left << setw(6) <<RP_offsetEl_s           <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(44) <<"* X offset of beam crossing in micrometers:"<<""
          << left << setw(6) <<RP_cross_x           <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(44) <<"* Y offset of beam crossing in micrometers:"<<""
          << left << setw(6) <<RP_cross_y           <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(44) <<"* X Angle of beam crossing:"<<""
          << left << setw(6) <<RP_cross_ang_x          <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(44) <<"* 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 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_Shcal       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* N term for central HCAL: "<<""
        << left << setw(30) <<HAD_Nhcal       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* C term for central HCAL: "<<""
        << left << setw(30) <<HAD_Chcal       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* S term for FCAL: "<<""
        << left << setw(30) <<HAD_Shf         <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* N term for FCAL: "<<""
        << left << setw(30) <<HAD_Nhf         <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* C term for FCAL: "<<""
        << left << setw(30) <<HAD_Chf         <<""<< 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 tracks [GeV]: "<<""
        << left << setw(20) <<ISOL_PT               <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(40) <<"* Cone for isolation criteria: "<<""
        << left << setw(20) <<ISOL_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 efficiencies [%]                                         *"<<"\n";
  f_out<<"#***************************                                         *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out << left << setw(50) <<"* Efficiency to tag a \"b\" as a b-jet: "<<""
        << left << setw(10) <<BTAG_b               <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(50) <<"* Efficiency to mistag a c-jet as a b-jet: "<<""
        << left << setw(10) <<BTAG_mistag_c            <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(50) <<"* Efficiency to mistag a light jet as a b-jet: "<<""
        << left << setw(10) <<BTAG_mistag_l            <<""<< right << setw(10)<<"*"<<"\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) ));
  } 
  if(fabs(electron.Eta()) > CEN_max_tracker && 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_Nhcal,2) +
                                              pow(HAD_Chcal*energy_hcal,2) + 
                                              pow(HAD_Shcal*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);
  } 
  if(fabs(hadron.Eta()) > CEN_max_calo_cen && fabs(hadron.Eta()) < CEN_max_calo_fwd){ 
    energyS = gRandom->Gaus(energy, sqrt(
                            pow(HAD_Nhf,2) +
                            pow(HAD_Chf*energy,2) +
                            pow(HAD_Shf*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) {
  if(      rand()%100 < (BTAG_b+1)    && Bjets(subarray,JET.Eta(),JET.Phi())==pB ) return true; // b-tag of b-jets is 40%
  else if( rand()%100 < (BTAG_mistag_c+1) && Bjets(subarray,JET.Eta(),JET.Phi())==pC ) return true; // b-tag of c-jets is 10%
  else if( rand()%100 < (BTAG_mistag_l+1) && Bjets(subarray,JET.Eta(),JET.Phi())!=0)   return true; // b-tag of light jets is 1%
  return false;
}

//***********************Isolation criteria***********************
//****************************************************************
bool RESOLution::Isolation(const D_Particle& part, const vector<TRootTracks> &tracks, const float& pt_second_track, const float& isolCone, 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 > deltar) ||
             (genDeltaR==0)  // rejets the track of the particle itself
           ) continue ;
           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_Cone) ptiso += tracks[i].PT; // dR cut on tracks
      }
   if(deltar > isolCone) isolated = true;
   return isolated;
}

// ******* Calorimetric isolation
float RESOLution::CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers, const float iPhi, const float iEta) {
 // etrat, which is a percentage between 00 and 99. It is the ratio of the transverse energy 
 // in a 3×3 grid surrounding the muon to the pT of the muon. For well-isolated muons, both ptiso and etrat will be small.
   if(ISOL_Calo_ET>1E10) return UNDEFINED; // avoid doing anything unreasonable...
   float et_sum=0;
   // available parameters: ISOL_Calo_ET ,  ISOL_Calo_Grid
   // Get the EtaCalo/PhiCalo of the muon ; 
   // transform it into iEta/iPhi to get the towers, and their neighbourh (i-1, i-2, etc)

   unsigned int N = ISOL_Calo_Grid;
   int index= iUNDEFINED; // index of the central tower of the grid in TOWER_eta_edges[.]; 
                          // !! TOWER_eta_edges is only with eta>0
   // finds the index of the central tower of the NxN grid
   for (unsigned int i=1; i< TOWER_number+1; i++) {
       if(fabs(iEta) >= TOWER_eta_edges[i-1] && fabs(iEta) < TOWER_eta_edges[i]) {
           index = i-1;
           break;
       }
   }
   if(index != iUNDEFINED) {
           // finds the size in phi of the cells for this eta
           float dphi = TOWER_dphi[index]*pi/180.; // in rad

           //cout << "Grid " << " ----------\n";
           for (unsigned int i_eta=0; i_eta<N; i_eta++) {
              unsigned int real_index = (iEta>0) ? index+i_eta-(N-1)/2 : index+1+i_eta-(N-1)/2 ;
              float eta_ith_tower = TOWER_eta_edges[real_index];
              if(iEta<0) eta_ith_tower *= -1;
                                        
              for (unsigned int i_phi=0; i_phi<N; i_phi++) {
                 float phi_ith_tower = iPhi + (float)(i_phi - (N-1.)/2.)*dphi;
                 D_CaloTower calMuon(towers.getElement(eta_ith_tower,phi_ith_tower));
                 if(calMuon.getEta() != UNDEFINED && calMuon.getE() > ISOL_Calo_ET) { 
                    et_sum += calMuon.getE();
                    //cout << "eta/phi = " << eta_ith_tower << "\t" << phi_ith_tower << "\t" << calMuon.getE() << " GeV" << endl;
                 } 
                 //else cout << "eta/phi = " << eta_ith_tower << "\t" << phi_ith_tower << "\tnot active\n";
              }
    
           } // NxN grid
        //cout << "-----------" << etrat << endl;
   } 
   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
  //cout << "etrat = " << et_sum << "\t Pt=" << part.Pt() << endl;
  
  // should return a number between 0 and 99 (due to LHCO definitions)
  // which is Pt(muon) / sum(ET)
  float etrat = 0.;
  if(et_sum==0) etrat = 99.;
  else if(et_sum>0) etrat = 100*part.Pt()/et_sum;
  if(etrat<0) cout << "Error: negative etrat in CaloIsolation (" << etrat <<")\n";
  //else if(etrat>99) cout << "Error: etrat should be in [0;99] in CaloIsolation (" << etrat <<")\n";
  if(etrat>99) etrat = 99;
  return etrat; 
}


  //********** returns a segmented value for eta and phi, for calo towers *****
void 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;
        }
   }
   if(index==UNDEFINED) return; 
   iPhi = UNDEFINED;
   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;
}

//**************************** 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);
    int ID; 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(22) <<"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<<"*                                                                    *"<<"\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.7 of Delphes ---                 **"<< endl;
  cout <<"**                Last date of change:    7 May 2009               **"<< 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 is a beta version of Delphes and     **"<< endl;
  cout <<"** therefore comes without guarantees. Beware of errors and please **"<< endl;
  cout <<"**        give us your feedbacks about potential bugs              **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"*********************************************************************"<< endl;
  cout <<"*********************************************************************"<< endl;

}