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

/*
  ---- Delphes ----
  A Fast Simulator for general purpose LHC detector
  S. Ovyn ~~~~ severine.ovyn@uclouvain.be

  Center for Particle Physics and Phenomenology (CP3)
  Universite Catholique de Louvain (UCL)
  Louvain-la-Neuve, Belgium
*/

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


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

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


ParticleUtil::ParticleUtil(const TLorentzVector &genMomentum, int pid) {
  _pid=pid;
  _e = genMomentum.E();
  _px = genMomentum.Px();
  _py = genMomentum.Py();
  _pz = genMomentum.Pz();
  _pt = genMomentum.Pt();

  //_e, _px, _py, _pz, _pt;
  //float _eta, _etaCalo, _phi, _phiCalo;
  //int _pid;
}

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

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_Cfwd          = 0.107;             // S term for forward ECAL
  ELG_Sfwd          = 2.084;             // C term for forward ECAL
  ELG_Nfwd          = 0.0;               // N term for central ECAL

  // 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 // hadronic calorimeter
  HAD_Nhcal         = 0.;                // N term for central HCAL
  HAD_Chcal         = 0.05;              // C term for central HCAL
  HAD_Shf           = 2.7;               // S term for HF // forward calorimeter
  HAD_Nhf           = 0.;                // N term for HF
  HAD_Chf           = 0.13;              // C term for HF

  // Muon smearing
  MU_SmearPt        = 0.01;

  // 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
  PTCUT_elec      = 10.0;
  PTCUT_muon      = 10.0;
  PTCUT_jet       = 20.0;
  PTCUT_gamma     = 10.0;
  PTCUT_taujet    = 10.0;

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

  // 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_trigger     = 1;                       //1 to run the trigger selection else 0
  FLAG_frog        = 1;                       //1 to run the FROG event display

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

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

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


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;

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

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

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

  // 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_trigger     = DET.FLAG_trigger;
  FLAG_frog        = DET.FLAG_frog;

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

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

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

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;

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

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

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

  // 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_trigger     = DET.FLAG_trigger;
  FLAG_frog        = DET.FLAG_frog;

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

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

  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;
  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 <<"**        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;
    
    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(),"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(),"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(),"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(),"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(),"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_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(),"NEvents_Frog"))     {curstring >> varname >> ivalue; NEvents_Frog     = ivalue;}
  }
  
  //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<<"#*********************************************************************"<<"\n";
  f_out<<"#                                                                    *"<<"\n";
  f_out<<"#                  ---- DELPHES release 1.0 ----                     *"<<"\n";
  f_out<<"#                                                                    *"<<"\n";
  f_out<<"#        A Fast Simulator for general purpose LHC detector           *"<<"\n";
  f_out<<"#                 Written by S. Ovyn and X. Rouby                    *"<<"\n";
  f_out<<"#                   severine.ovyn@uclouvain.be                       *"<<"\n";
  f_out<<"#                                                                    *"<<"\n";
  f_out<<"#                   http:                                            *"<<"\n";
  f_out<<"#                                                                    *"<<"\n";
  f_out<<"#         Center for Particle Physics and Phenomenology (CP3)        *"<<"\n";
  f_out<<"#              Universite Catholique de Louvain (UCL)                *"<<"\n";
  f_out<<"#                     Louvain-la-Neuve, Belgium                      *"<<"\n";
  f_out<<"#                                                                    *"<<"\n";
  f_out<<"#....................................................................*"<<"\n";
  f_out<<"#                                                                    *"<<"\n";
  f_out<<"#     This package uses:                                             *"<<"\n";
  f_out<<"#     FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210]    *"<<"\n";
  f_out<<"#     Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]     *"<<"\n";
  f_out<<"#     ExRootAnalysis                                                 *"<<"\n";
  f_out<<"#                                                                    *"<<"\n";
  f_out<<"#....................................................................*"<<"\n";
  f_out<<"#                                                                    *"<<"\n";
  f_out<<"# This file contains all the running parameters (detector and cuts)  *"<<"\n";
  f_out<<"# necessary to reproduce the detector simulation                     *"<<"\n";
  f_out<<"#                                                                    *"<<"\n";
  f_out<<"#....................................................................*"<<"\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_vfd==1){
    f_out<<"#**********************************                                  *"<<"\n";
    f_out<<"# Very forward detector switches 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 << left << setw(55) <<"* Distance of the RP to the IP, in meters:   "<<""
          << left << setw(5) <<RP_220_s             <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* Distance of the RP to the beam, in meters: "<<""
          << left << setw(5) <<RP_220_x             <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* Distance of the RP to the IP, in meters:   "<<""
          << left << setw(5) <<RP_420_s             <<""<< right << setw(10)<<"*"<<"\n";
    f_out << left << setw(55) <<"* Distance of the RP to the beam, in meters: "<<""
          << left << setw(5) <<RP_420_x             <<""<< right << setw(10)<<"*"<<"\n";
    f_out<<"*                                                                    *"<<"\n";
  }
  else {
    f_out<<"#***********************************                                 *"<<"\n";
    f_out<<"# Very forward detector switches 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 forward ECAL: "<<""
        << left << setw(30) <<ELG_Sfwd       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* N term for forward ECAL: "<<""
        << left << setw(30) <<ELG_Nfwd       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* C term for forward ECAL: "<<""
        << left << setw(30) <<ELG_Cfwd       <<""<< 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 forward HCAL: "<<""
        << left << setw(30) <<HAD_Shf         <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* N term for forward HCAL: "<<""
        << left << setw(30) <<HAD_Nhf         <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(30) <<"* C term for forward HCAL: "<<""
        << left << setw(30) <<HAD_Chf         <<""<< right << setw(10)<<"*"<<"\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 switches 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 switches 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<<"*                                                                    *"<<"\n";
  f_out<<"#***************                                                     *"<<"\n";
  f_out<<"# Jet definition                                                     *"<<"\n";
  f_out<<"#***************                                                     *"<<"\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";
  
}

// **********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)
{
  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];
               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 !!!!!
// previously, the argument 'phi' was before the argument 'eta'
// this has been changed for consistency with the other functions
// double check your running code that uses NumTracks !
unsigned int RESOLution::NumTracks(const vector<TLorentzVector> &tracks, const float pt_track, const float eta, const float phi) {
  unsigned int numtrack=0;
  for(unsigned int i=0; i < tracks.size(); i++) { 
      if((tracks[i].Pt() < pt_track )||
         (DeltaR(phi,eta,tracks[i].Phi(),tracks[i].Eta()) > TAU_track_scone)
        )continue;
      numtrack++;
  }
  return numtrack;
}


//*** 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<TRootGenParticle> &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<TRootGenParticle> &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 float phi, const float eta,const vector<TLorentzVector> &tracks, const float pt_second_track)
{
   bool isolated = false;
   float deltar=5000.; // Initial value; should be high; no further repercussion
   // loop on all final charged particles, with p_t >2, close enough from the electron
   for(unsigned int i=0; i < tracks.size(); i++) 
      {
       if(tracks[i].Pt() < pt_second_track)continue;
       float genDeltaR = DeltaR(phi,eta,tracks[i].Phi(),tracks[i].Eta()); 
         if(
             (genDeltaR > deltar) ||
             (genDeltaR==0)
           ) continue ;
           deltar=genDeltaR;
      }
   if(deltar > 0.5) isolated = true;
   return isolated;
}


  //********** 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 = -100; 
   int index=-100;
   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];
                index = i-1;
                //cout << setw(15) << left << eta << "\t" << iEta << endl;
                break;
        }
   }
   if(index==-100) return; 
   iPhi = -100;
   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; 
                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 Charge(const int pid)
{
  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;
            
}