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



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

TRACKING_RADIUS  = 129;                   //radius of the BField coverage
TRACKING_LENGTH  = 300;                   //length of the BField coverage
BFIELD_X         = 0.0;
BFIELD_Y         = 0.0;
BFIELD_Z         = 3.8;

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;
GAMMA_pt =        10.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;

//Trigger flag
DOTRIGGER    =     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;

// Define Calorimeter Towers
NTOWERS	= 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[NTOWERS+1];
for(unsigned int i=0; i<NTOWERS+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[NTOWERS];
for(unsigned int i=0; i<NTOWERS; i++) TOWER_DPHI[i] = tower_dphi[i];

}

//------------------------------------------------------------------------------
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; int ivalue;
    
    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(),"TRACKING_RADIUS")){curstring >> varname >> ivalue; TRACKING_RADIUS = ivalue;}
    else if(strstr(temp_string.c_str(),"TRACKING_LENGTH")){curstring >> varname >> ivalue; TRACKING_LENGTH = ivalue;}
    else if(strstr(temp_string.c_str(),"BFIELD_X")){curstring >> varname >> value; BFIELD_X = value;}
    else if(strstr(temp_string.c_str(),"BFIELD_Y")){curstring >> varname >> value; BFIELD_Y = value;}
    else if(strstr(temp_string.c_str(),"BFIELD_Z")){curstring >> varname >> value; BFIELD_Z = 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 >> ivalue; TAGGING_B = ivalue;}
    else if(strstr(temp_string.c_str(),"MISTAGGING_C")){curstring >> varname >> ivalue; MISTAGGING_C = ivalue;}
    else if(strstr(temp_string.c_str(),"MISTAGGING_L")){curstring >> varname >> ivalue; MISTAGGING_L = ivalue;}
    else if(strstr(temp_string.c_str(),"CONERADIUS")){curstring >> varname >> value; CONERADIUS = value;}
    else if(strstr(temp_string.c_str(),"JETALGO")){curstring >> varname >> ivalue; JETALGO = ivalue;}
    else if(strstr(temp_string.c_str(),"TRACKING_EFF")){curstring >> varname >> ivalue; TRACKING_EFF = ivalue;}
    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(),"GAMMA_pt")){curstring >> varname >> value; GAMMA_pt = value;}
    else if(strstr(temp_string.c_str(),"TAUJET_pt")){curstring >> varname >> value; TAUJET_pt = value;}
    else if(strstr(temp_string.c_str(),"NTOWERS")){curstring >> varname >> ivalue; NTOWERS = ivalue;}
    else if(strstr(temp_string.c_str(),"DOTRIGGER")){curstring >> varname >> ivalue; DOTRIGGER = ivalue;}
    else if(strstr(temp_string.c_str(),"TOWER_ETA_EDGES")){
	curstring >> varname;	for(unsigned int i=0; i<NTOWERS+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<NTOWERS; i++) {curstring >> value; TOWER_DPHI[i] = 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) { 
//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) <<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<<"# Magnetic field needed informations                                 *"<<"\n";
  f_out<<"#***********************************                                 *"<<"\n";
  f_out<<"*                                                                    *"<<"\n";
  f_out << left << setw(55) <<"* Radius of the BField coverage:              "<<""
        << left << setw(5) <<TRACKING_RADIUS <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(55) <<"* Length of the BField coverage:              "<<""
        << left << setw(5) <<TRACKING_LENGTH <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(55) <<"* BField X component:                         "<<""
        << left << setw(5) <<BFIELD_X <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(55) <<"* BField Y component:                         "<<""
        << left << setw(5) <<BFIELD_Y <<""<< right << setw(10)<<"*"<<"\n";
  f_out << left << setw(55) <<"* BField Z component:                         "<<""
        << left << setw(5) <<BFIELD_Z <<""<< right << setw(10)<<"*"<<"\n";
  f_out<<"*                                                                    *"<<"\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) << NTOWERS <<""<< 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<NTOWERS+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=(NTOWERS+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<NTOWERS; i++) {
		f_out << left << setw(7) << TOWER_DPHI[i];
		if(!( (i+1) %9 )) f_out << right << setw(3) << "*" << "\n" << "*   ";
	  }
	  for (unsigned int i=(NTOWERS)%9; i<9; i++) f_out << left << setw(7) << "";
	  f_out << right << setw(3)<<"*"<<"\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 << left << setw(40) <<"* Minimum pT for photons: "<<""
        << left << setw(20) <<GAMMA_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) ));
  } 
  if(fabs(electron.Eta()) > MAX_TRACKER && fabs(electron.Eta()) < 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()) < 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()) < 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(abs(hadron.Eta()) > MAX_CALO_CEN && fabs(hadron.Eta()) < 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() < 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 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< NTOWERS+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(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;
            
}