/*
  ---- 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 "interface/SmearUtil.h"
#include "TRandom.h"

#include <iostream>
#include <sstream>
#include <fstream>

#include <cstdio>
#include <cstdio>
#include <TCanvas.h>
#include <TLorentzVector.h>
#include <string>
#include <fstream>
#include <istream>
#include <iostream>
#include <TMath.h>
#include <vector>
#include "TRandom.h"
#include <cstdlib>
#include <math.h>
#include <time.h>

#include <iostream>
#include <iomanip>

#include <iostream>
#include <sstream>

#include <fstream>
#include <vector>

using namespace std;

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

RESOLution::RESOLution() {

MAX_TRACKER  = 	2.5;  // tracker coverage
MAX_CALO_CEN =  3.0;  // central calorimeter coverage
MAX_CALO_FWD =  5.0;  // forward calorimeter pseudorapidity coverage
MAX_MU       =  2.4;  // muon chambers pseudorapidity coverage
MIN_CALO_VFWD=  5.2;  // very forward calorimeter (if any), like CASTOR
MAX_CALO_VFWD=  6.6;  // very forward calorimeter (if any), like CASTOR
MIN_ZDC      =  8.3;  // zero-degree calorimeter, coverage

ZDC_S	     =  140.; // ZDC  distance to IP
RP220_S      =  220;  // distance of the RP to the IP, in meters
RP220_X      =  0.002;// distance of the RP to the beam, in meters
FP420_S      =  420;  // distance of the RP to the IP, in meters
FP420_X      =  0.004;// distance of the RP to the beam, in meters


ELG_Scen =        0.05;                  // S term for central ECAL
ELG_Ncen =        0.25 ;                 // N term for central ECAL
ELG_Ccen =        0.0055 ;                // 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

HAD_Shcal =       1.5 ;                    // S term for central HCAL // hadronic calorimeter
HAD_Nhcal  =      0.0   ;                    // N term for central HCAL
HAD_Chcal  =      0.05 ;                   // C term for central HCAL
HAD_Shf  =        2.7  ;                    // S term for central HF // forward calorimeter
HAD_Nhf  =        0.0   ;                    // N term for central HF
HAD_Chf  =        0.13 ;                    // C term for central HF

MU_SmearPt =      0.01 ;

ELEC_pt =         10.0;
MUON_pt =         10.0;
JET_pt =          20.0;
TAUJET_pt =       10.0;


TAU_CONE_ENERGY = 0.15  ;                  // Delta R = radius of the cone // for "electromagnetic collimation"
TAU_EM_COLLIMATION = 0.95;
TAU_CONE_TRACKS=  0.4  ;                   //Delta R for tracker isolation for tau's
PT_TRACK_TAU =    2.0  ;                   // GeV // 6 GeV ????


PT_TRACKS_MIN =   0.9  ;                   // minimal pt needed to reach the calorimeter, in GeV
PT_QUARKS_MIN =   2.0  ;                   // minimal pt needed by quarks to reach the tracker, in GeV (??????)
TRACKING_EFF  =   90;


TAGGING_B    =    40;
MISTAGGING_C =    10;
MISTAGGING_L =    1;


CONERADIUS   = 0.7; 	// generic jet radius ; not for tau's !!!
JETALGO      =   1;     // 1 for Cone algorithm, 2 for MidPoint algorithm, 3 for SIScone algorithm, 4 for kt algorithm

//General jet parameters
SEEDTHRESHOLD     = 1.0;
OVERLAPTHRESHOLD = 0.75;

// Define Cone algorithm.
C_ADJACENCYCUT     = 2;
C_MAXITERATIONS    = 100;
C_IRATCH           = 1;

//Define MidPoint algorithm.
M_CONEAREAFRACTION  =   0.25;
M_MAXPAIRSIZE       =   2;
M_MAXITERATIONS     =   100;

}

//------------------------------------------------------------------------------
void RESOLution::ReadDataCard(const string datacard) {
  
  string temp_string;
  istringstream curstring;
  
  ifstream fichier_a_lire(datacard.c_str());
  if(!fichier_a_lire.good()) {
	cout << datacard << "Datadard " << datacard << " 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;
    
    if(strstr(temp_string.c_str(),"#")) { }
    else if(strstr(temp_string.c_str(),"MAX_TRACKER")){curstring >> varname >> value; MAX_TRACKER = value;}
    else if(strstr(temp_string.c_str(),"MAX_CALO_CEN")){curstring >> varname >> value; MAX_CALO_CEN = value;}
    else if(strstr(temp_string.c_str(),"MAX_CALO_FWD")){curstring >> varname >> value; MAX_CALO_FWD = value;}
    else if(strstr(temp_string.c_str(),"MAX_MU")){curstring >> varname >> value; MAX_MU = 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(),"TAU_CONE_ENERGY")){curstring >> varname >> value; TAU_CONE_ENERGY = value;}
    else if(strstr(temp_string.c_str(),"TAU_CONE_TRACKS")){curstring >> varname >> value; TAU_CONE_TRACKS = value;}
    else if(strstr(temp_string.c_str(),"PT_TRACK_TAU")){curstring >> varname >> value; PT_TRACK_TAU = value;}
    else if(strstr(temp_string.c_str(),"PT_TRACKS_MIN")){curstring >> varname >> value; PT_TRACKS_MIN = value;}
    else if(strstr(temp_string.c_str(),"TAGGING_B")){curstring >> varname >> value; TAGGING_B = (int)value;}
    else if(strstr(temp_string.c_str(),"MISTAGGING_C")){curstring >> varname >> value; MISTAGGING_C = (int)value;}
    else if(strstr(temp_string.c_str(),"MISTAGGING_L")){curstring >> varname >> value; MISTAGGING_L = (int)value;}
    else if(strstr(temp_string.c_str(),"CONERADIUS")){curstring >> varname >> value; CONERADIUS = value;}
    else if(strstr(temp_string.c_str(),"JETALGO")){curstring >> varname >> value; JETALGO = (int)value;}
    else if(strstr(temp_string.c_str(),"TRACKING_EFF")){curstring >> varname >> value; TRACKING_EFF = (int)value;}
    else if(strstr(temp_string.c_str(),"ELEC_pt")){curstring >> varname >> value; ELEC_pt = value;}
    else if(strstr(temp_string.c_str(),"MUON_pt")){curstring >> varname >> value; MUON_pt = value;}
    else if(strstr(temp_string.c_str(),"JET_pt")){curstring >> varname >> value; JET_pt = value;}
    else if(strstr(temp_string.c_str(),"TAUJET_pt")){curstring >> varname >> value; TAUJET_pt = value;}

  }

// General jet variables
 SEEDTHRESHOLD     = 1.0;
 OVERLAPTHRESHOLD = 0.75;

// Define Cone algorithm.
 C_ADJACENCYCUT     = 2;
 C_MAXITERATIONS    = 100;
 C_IRATCH           = 1;

//Define MidPoint algorithm.
 M_CONEAREAFRACTION  =   0.25;
 M_MAXPAIRSIZE       =   2;
 M_MAXITERATIONS     =   100;

//Define SISCone algorithm.
 NPASS = 0;               
 PROTOJET_PTMIN = 0.0; 


}

void RESOLution::Logfile(string LogName) { 
  
  ofstream f_out(LogName.c_str());
  
  f_out<<"#*********************************************************************"<<"\n";
  f_out<<"#                                                                    *"<<"\n";
  
  string temp_string;
  istringstream curstring;
  ifstream fichier_a_lire("VERSION");
  if(fichier_a_lire.good()) {
  string Version;
  while (getline(fichier_a_lire,temp_string)) {
    curstring.str(temp_string);
    curstring >> Version;}
    f_out<<"#"<< right << setw(39) <<"---- DELPHES release "<<Version<< left <<" ----                     *"<<"\n";
    
  }
  else { f_out<<"#                ---- DELPHES release not defined ----               *"<<"\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) <<MAX_TRACKER   <<""<< right << setw(15)<<"*"<<"\n";
  f_out << left << setw(30) <<"* Maximum central calorimeter:               "<<""
	<< left << setw(10) <<MAX_CALO_CEN  <<""<< right << setw(15)<<"*"<<"\n";
  f_out << left << setw(30) <<"* Maximum forward calorimeter:               "<<""
	<< left << setw(10) <<MAX_CALO_FWD  <<""<< right << setw(15)<<"*"<<"\n";
  f_out << left << setw(30) <<"* Muon chambers coverage:                    "<<""
	<< left << setw(10) <<MAX_MU        <<""<< right << setw(15)<<"*"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out<<"#*************************************                               *"<<"\n";
  f_out<<"# Very forward detector caracteristics                               *"<<"\n";
  f_out<<"#*************************************                               *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out << left << setw(55) <<"* Minimum very forward calorimeter:          "<<""
	<< left << setw(5) <<MIN_CALO_VFWD <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(55) <<"* Maximum very forward calorimeter:          "<<""
	<< left << setw(5) <<MAX_CALO_VFWD <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(55) <<"* Distance of the ZDC to the IP, in meters:  "<<""
	<< left << setw(5) <<ZDC_S         <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(55) <<"* Distance of the RP to the IP, in meters:   "<<""
	<< left << setw(5) <<RP220_S       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(55) <<"* Distance of the RP to the beam, in meters: "<<""
	<< left << setw(5) <<RP220_X       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(55) <<"* Distance of the RP to the IP, in meters:   "<<""
	<< left << setw(5) <<FP420_S       <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(55) <<"* Distance of the RP to the beam, in meters: "<<""
	<< left << setw(5) <<FP420_X       <<""<< right << setw(10)<<"*"<<"\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<<"# Tracking system acceptance                                         *"<<"\n";
  f_out<<"#***************************                                         *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out << left << setw(55) <<"* Minimal pT needed to reach the calorimeter [GeV]: "<<""
	<< left << setw(10) <<PT_TRACKS_MIN         <<""<< right << setw(5)<<"*"<<"\n";
  f_out << left << setw(55) <<"* Efficiency associated to the tracking: "<<""
	<< left << setw(10) <<TRACKING_EFF          <<""<< right << setw(5)<<"*"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out<<"#*************************                                           *"<<"\n";
  f_out<<"# Muon smearing parameters                                           *"<<"\n";
  f_out<<"#*************************                                           *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  //MU_SmearPt       0.01
  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_CONE_ENERGY                           <<""<< right << setw(20)<<"*"<<"\n";
  f_out << left << setw(45) <<"* Fraction of energy in the small cone: "<<""
        << left << setw(5) <<TAU_EM_COLLIMATION*100                    <<""<< right << setw(20)<<"! not in datacard  *"<<"\n";
  f_out << left << setw(45) <<"* Cone radius for tracking tagging:     "<<""
	<< left << setw(5) <<TAU_CONE_TRACKS                           <<""<< right << setw(20)<<"*"<<"\n";
  f_out << left << setw(45) <<"* Minimum track pT [GeV]:               "<<""
	<< left << setw(5) <<PT_TRACK_TAU                              <<""<< right << setw(20)<<"*"<<"\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) <<ELEC_pt         <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(40) <<"* Minimum pT for muons: "<<""
	<< left << setw(20) <<MUON_pt         <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(40) <<"* Minimum pT for jets: "<<""
	<< left << setw(20) <<JET_pt         <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(40) <<"* Minimum pT for Tau-jets: "<<""
	<< left << setw(20) <<TAUJET_pt         <<""<< right << setw(10)<<"*"<<"\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) <<TAGGING_B         <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(50) <<"* Efficiency to mistag a c-jet as a b-jet: "<<""
	<< left << setw(10) <<MISTAGGING_C         <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(50) <<"* Efficiency to mistag a light jet as a b-jet: "<<""
	<< left << setw(10) <<MISTAGGING_L         <<""<< 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(JETALGO) {
  default:
  case 1: {
  f_out<<"* CDF JetClu jet algorithm with parameters:                          *"<<"\n";
  f_out << left << setw(40) <<"*        - Seed threshold:    "<<""
        << left << setw(10) <<SEEDTHRESHOLD                   <<""<< right << setw(20)<<"! not in datacard  *"<<"\n";
  f_out << left << setw(40) <<"*        - Cone radius:       "<<""
        << left << setw(10) <<CONERADIUS                      <<""<< right << setw(20)<<"*"<<"\n";
  f_out << left << setw(40) <<"*        - Adjacency cut:     "<<""
        << left << setw(10) <<C_ADJACENCYCUT                  <<""<< right << setw(20)<<"! not in datacard  *"<<"\n";
  f_out << left << setw(40) <<"*        - Max iterations:    "<<""
        << left << setw(10) <<C_MAXITERATIONS                 <<""<< right << setw(20)<<"! not in datacard  *"<<"\n";
  f_out << left << setw(40) <<"*        - Iratch:            "<<""
        << left << setw(10) <<C_IRATCH                        <<""<< right << setw(20)<<"! not in datacard  *"<<"\n";
  f_out << left << setw(40) <<"*        - Overlap threshold: "<<""
        << left << setw(10) <<OVERLAPTHRESHOLD                <<""<< 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) <<SEEDTHRESHOLD                   <<""<< right << setw(10)<<"! not in datacard  *"<<"\n";
  f_out << left << setw(40) <<"*        - Cone radius:       "<<""
        << left << setw(20) <<CONERADIUS                      <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(40) <<"*        - Cone area fraction:"<<""
        << left << setw(20) <<M_CONEAREAFRACTION              <<""<< right << setw(10)<<"! not in datacard  *"<<"\n";
  f_out << left << setw(40) <<"*        - Maximum pair size: "<<""
        << left << setw(20) <<M_MAXPAIRSIZE                   <<""<< right << setw(10)<<"! not in datacard  *"<<"\n";
  f_out << left << setw(40) <<"*        - Max iterations:    "<<""
        << left << setw(20) <<M_MAXITERATIONS                 <<""<< right << setw(10)<<"! not in datacard  *"<<"\n";
  f_out << left << setw(40) <<"*        - Overlap threshold: "<<""
        << left << setw(20) <<OVERLAPTHRESHOLD                <<""<< 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) <<CONERADIUS                            <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(40) <<"*        - Overlap threshold:       "<<""
        << left << setw(20) <<OVERLAPTHRESHOLD                      <<""<< right << setw(10)<<"! not in datacard  *"<<"\n";
  f_out << left << setw(40) <<"*        - Number pass max:         "<<""
	<< left << setw(20) <<NPASS                                 <<""<< right << setw(10)<<"! not in datacard  *"<<"\n";
  f_out << left << setw(40) <<"*        - Minimum pT for protojet: "<<""
	<< left << setw(20) <<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) <<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) <<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) <<CONERADIUS                <<""<< right << setw(10)<<"*"<<"\n";
  }
    break;


  }
  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()) < MAX_TRACKER) { // if the electron is inside the tracker
    energyS = gRandom->Gaus(energy, sqrt( 
					 pow(ELG_Ncen,2) + 
					 pow(ELG_Ccen*energy,2) + 
					 pow(ELG_Scen*sqrt(energy),2) ));
  } else { // outside the tracker
    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 = 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()) < 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);
  } else { 
    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);
}

// **********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() < SEEDTHRESHOLD) continue;
      if((DeltaR(phi,eta,towers[i].fourVector.phi(),towers[i].fourVector.eta()) < TAU_CONE_ENERGY)) {
          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_CONE_TRACKS)
        )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 < 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 < (TAGGING_B+1)    && Bjets(subarray,JET.Eta(),JET.Phi())==pB ) return true; // b-tag of b-jets is 40%
  else if( rand()%100 < (MISTAGGING_C+1) && Bjets(subarray,JET.Eta(),JET.Phi())==pC ) return true; // b-tag of c-jets is 10%
  else if( rand()%100 < (MISTAGGING_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(Float_t phi,Float_t eta,const vector<TLorentzVector> &tracks,float PT_TRACK2)
{
   bool isolated = false;
   Float_t 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_TRACK2)continue;
       Float_t genDeltaR = DeltaR(phi,eta,tracks[i].Phi(),tracks[i].Eta());    // slower to evaluate
         if(
             (genDeltaR > deltar) ||
             (genDeltaR==0)
           ) continue ;
           deltar=genDeltaR;
      }
   if(deltar > 0.5)isolated = true; // returns the closest distance
   return isolated;
}


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


float Charge(const long int pid) {
   // source: RPP chap 34 Monte Carlo Particle Numbering Scheme
/*   switch (abs(pid)) {
	case 1: case 3: case 5: case 7: return (float) sign(pid)*(-1/3); break;   // d, s, b, b'
	case 2: case 4: case 6: case 8: return (float) sign(pid)*2/3; break;      // u, c, t, t'

	case 11: case 13: case 15:      return (float) sign(pid)*(-1); break;     // e, mu, tau
	case 12: case 14: case 16:      return (float) 0; break;     // nu_e, nu_mu, nu_tau

	case 9: case 21: case 22: case 23: case 25: 
                case 32: case 33: case 35: case 36: return (float) 0; break; // neutral gauge/higgs bosons
	case 24: case 34: case 37: return (float) sign(pid); break; // charged gauge/higgs bosons
   }
*/
	return 0;
}