/*
  ---- 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>
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.028;                  // S term for central ECAL
ELG_Ncen =        0.124 ;                 // N term for central ECAL
ELG_Ccen =        0.0026 ;                // 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_Secal =       0.05 ;                    // S term for central ECAL // electromagnetic calorimeter
HAD_Necal =       0.25 ;                    // N term for central ECAL
HAD_Cecal =       0.0055 ;                  // C term for central ECAL
HAD_Shcal =       0.91 ;                    // S term for central HCAL // hadronic calorimeter
HAD_Nhcal  =      0.   ;                    // N term for central HCAL
HAD_Chcal  =      0.038 ;                   // C term for central HCAL
HAD_Shf  =        2.7  ;                    // S term for central HF // forward calorimeter
HAD_Nhf  =        0.   ;                    // N term for central HF
HAD_Chf  =        0.13 ;                    // C term for central HF

MU_SmearPt =      0.01 ;

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
// Define Cone algorithm.
C_SEEDTHRESHOLD     = 1.0;
C_ADJACENCYCUT     = 2;
C_MAXITERATIONS    = 100;
C_IRATCH           = 1;
C_OVERLAPTHRESHOLD = 0.75;

//Define MidPoint algorithm.
M_SEEDTHRESHOLD     =   1.0;
M_CONEAREAFRACTION  =   0.25;
M_MAXPAIRSIZE       =   2;
M_MAXITERATIONS     =   100;
M_OVERLAPTHRESHOLD  =  0.75;

}

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

// Define Cone algorithm.
 C_SEEDTHRESHOLD     = 1.0;
 C_ADJACENCYCUT     = 2;
 C_MAXITERATIONS    = 100;
 C_IRATCH           = 1;
 C_OVERLAPTHRESHOLD = 0.75;

//Define MidPoint algorithm.
 M_SEEDTHRESHOLD     =   1.0;
 M_CONEAREAFRACTION  =   0.25;
 M_MAXPAIRSIZE       =   2;
 M_MAXITERATIONS     =   100;
 M_OVERLAPTHRESHOLD  =  0.75;

}


// **********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) ) ); 
  }
float theta=electron.Theta();
float phi=electron.Phi();
  float px = energyS * sin(theta) * cos(phi);
  float py = energyS * sin(theta) * sin(phi);
  float pz = energyS * cos(theta);
  //cout<<" px "<<px<<" py "<<py<<" pz "<<pz<<endl;
  electron.SetPtEtaPhiE(energyS/cosh(electron.Eta()), electron.Eta(), electron.Phi(), energyS);
  //cout<<"nous "<<electron.Px()<<"  "<<electron.Py()<<"  "<<electron.Pz()<<endl;
  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() < M_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;
}