Utils/src/PhotonTools.cc

// $Id: PhotonTools.cc,v 1.6 2011/06/29 18:28:05 fabstoec Exp $

#include "MitPhysics/Utils/interface/PhotonTools.h"
#include "MitPhysics/Utils/interface/ElectronTools.h"
#include "MitPhysics/Utils/interface/IsolationTools.h"
#include "MitAna/DataTree/interface/StableData.h"
#include <TFile.h>

ClassImp(mithep::PhotonTools)

using namespace mithep;

//--------------------------------------------------------------------------------------------------
PhotonTools::PhotonTools()  
{
  // Constructor.
}

//--------------------------------------------------------------------------------------------------
Bool_t PhotonTools::PassConversionId(const Photon *p, const DecayParticle *c) {

  if (!c) return kTRUE;
  
  ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
  Double_t deta = c->Eta()-dirconvsc.Eta();
  Double_t dphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
  Double_t eoverp = p->SCluster()->Energy()/c->P();
  
  if (p->IsEB() && eoverp>2.0) return kFALSE;
  if (p->IsEE() && eoverp>3.0) return kFALSE;
  
  if (p->IsEE() && TMath::Abs(deta)>0.01) return kFALSE;
  if (p->IsEE() && TMath::Abs(dphi)>0.01) return kFALSE;

  return kTRUE;
    
}

//--------------------------------------------------------------------------------------------------
Bool_t PhotonTools::PassElectronVeto(const Photon *p, const ElectronCol *els) {

  Bool_t pass = kTRUE;
  for (UInt_t i=0; i<els->GetEntries(); ++i) {
    const Electron *e = els->At(i);
    if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0) {
      pass = kFALSE;
    }
  }
  
  return pass;
}

//--------------------------------------------------------------------------------------------------
Double_t PhotonTools::ElectronVetoCiC(const Photon *p, const ElectronCol *els) {
  
  for (UInt_t i=0; i<els->GetEntries(); ++i) {
    const Electron *e = els->At(i);
    if (e->SCluster()==p->SCluster() ) {
      //if( e->GsfTrk()->NExpectedHitsInner()==0 && e->GsfTrk()->Pt() > 2.5 ) {
      if( e->GsfTrk()->NExpectedHitsInner()==0 ) {
        double dEta = e->DeltaEtaSuperClusterTrackAtVtx();
        double dPhi = e->DeltaPhiSuperClusterTrackAtVtx();
        double dR = TMath::Sqrt(dEta*dEta+dPhi*dPhi);
        return dR;
      }
    }
  }  
  return 99.;
}

//--------------------------------------------------------------------------------------------------
Bool_t PhotonTools::PassElectronVetoConvRecovery(const Photon *p, const ElectronCol *els, const DecayParticleCol *conversions, const BaseVertex *v) {

  Bool_t pass = kTRUE;
  for (UInt_t i=0; i<els->GetEntries(); ++i) {
    const Electron *e = els->At(i);
    if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0 && ElectronTools::PassConversionFilter(e, conversions, 
                                                         v, 0, 1e-6, 2.0, kFALSE, kFALSE) ) {
      pass = kFALSE;
    }
  }
  
  return pass;
}

//--------------------------------------------------------------------------------------------------
Bool_t PhotonTools::PassTriggerMatching(const Photon *p, const TriggerObjectCol *trigobjs)
{
  
  for (UInt_t i=0; i<trigobjs->GetEntries(); ++i) {
    const TriggerObject *trigobj = trigobjs->At(i);
    if (trigobj->TriggerType()==TriggerObject::TriggerCluster || trigobj->TriggerType()==TriggerObject::TriggerElectron || trigobj->TriggerType()==TriggerObject::TriggerPhoton) {
      if (MathUtils::DeltaR(p->SCluster(),trigobj)<0.3) {
        return kTRUE;
      }
    }
  }
  
  return kFALSE;
  
  
}

//--------------------------------------------------------------------------------------------------
const DecayParticle *PhotonTools::MatchedConversion(const Photon *p, const DecayParticleCol *conversions, 
                                               const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
                                               Double_t lxyMin, Double_t dRMin) {
  
  const DecayParticle *match = 0;
  Double_t rhosmallest = 999.;
  for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
    const DecayParticle *c = conversions->At(i);
    ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
    Double_t dr = MathUtils::DeltaR(*c,dirconvsc);
    Double_t rho = c->Position().Rho();
    if (dr<dRMin && rho<rhosmallest && c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
      Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
      Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
      if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
        rhosmallest = rho;
        match = c;
      }
    }
    
  }
  
  return match;
  
}

//--------------------------------------------------------------------------------------------------
const DecayParticle *PhotonTools::MatchedConversion(const Track *t, const DecayParticleCol *conversions, 
                                               const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
                                               Double_t lxyMin) {
  
  for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
    const DecayParticle *c = conversions->At(i);
    if (c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
      Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
      Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
      if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
        const Track *ct1 = dynamic_cast<const ChargedParticle*>(c->Daughter(0))->Trk();
        const Track *ct2 = dynamic_cast<const ChargedParticle*>(c->Daughter(1))->Trk();
        if (t==ct1 || t==ct2) return c;
      }
    }
    
  }
  
  return 0;
  
}

PhotonTools::DiphotonR9EtaCats PhotonTools::DiphotonR9EtaCat(const Photon *p1, const Photon *p2) {
  
  if (p1->IsEB() && p2->IsEB()) {
    if (p1->R9()>0.93 && p2->R9()>0.93) return kCat1;
    else return kCat2;
    
  }
  else {
    if (p1->R9()>0.93 && p2->R9()>0.93) return kCat3;
    else return kCat4; 
  }
  
}

PhotonTools::DiphotonR9EtaConversionCats PhotonTools::DiphotonR9EtaConversionCat(const Photon *p1, const Photon *p2, const DecayParticleCol *conversions, const BaseVertex *v) {
  
  const DecayParticle *conv1 = MatchedConversion(p1, conversions, v);
  const DecayParticle *conv2 = MatchedConversion(p2, conversions, v);
    
  if (p1->IsEB() && p2->IsEB()) {
    if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat1;
    else if (conv1||conv2) return kNewCat2;
    else return kNewCat3;
    
  }
  else {
    if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat4;
    else if (conv1||conv2) return kNewCat5;
    else return kNewCat6;
  }
  
}

PhotonTools::CiCBaseLineCats PhotonTools::CiCBaseLineCat(const Photon *p) {
  if( p->IsEB() ) {
    if ( p->R9() > 0.94 ) return kCiCCat1;
    else return kCiCCat2;
  } else {
    if ( p->R9() > 0.94 ) return kCiCCat3;
    else return kCiCCat4;
  }
}

//--------------------------------------------------------------------------------------------------
const DecayParticle *PhotonTools::MatchedCiCConversion(const Photon *p, const DecayParticleCol *conversions, 
                                                       Double_t dPhiMin,
                                                       Double_t dEtaMin,
                                                       Double_t dRMin,
                                                       bool print) {
  
  const DecayParticle *match = NULL;

  double minDeta = 999.;
  double minDphi = 999.;
  double minDR   = 999.;

  double phPhi = p->SCluster()->Phi();
  double phEta = p->SCluster()->Eta();

  if(print)
    std::cout<<"  --- conv match photon eta = "<<phEta<<"  phi = "<<phPhi<<std::endl;  


  for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
    const DecayParticle *c = conversions->At(i);

    if(print)
      std::cout<< "   c "<<i+1<<"  pt = "<<c->Pt()<<std::endl;

    if(c->Pt()   < 1. )    continue; // is this refittedPirMomentum?

    //ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
    //ThreeVector dirconvsc = p->CaloPos() - c->Position();

    double convPhi = c->Phi();
    double convEta = c->Eta();
    
    const ThreeVector wrong(0,0,0);
    double Zvertex = c->DzCorrected(wrong);
    // ------------------ FROM GLOBE ----------------------
    //---Definitions for ECAL
    const float R_ECAL           = 136.5;
    const float Z_Endcap         = 328.0;
    const float etaBarrelEndcap  = 1.479; 
    
    //---ETA correction
    float Theta = 0.0  ; 
    float ZEcal = R_ECAL*sinh(convEta)+Zvertex;
    
    if(ZEcal != 0.0) Theta = TMath::ATan(R_ECAL/ZEcal);
    if(Theta<0.0) Theta = Theta+M_PI;
    double fEta = - TMath::Log(TMath::Tan(0.5*Theta));
    
    if( fabs(fEta) > etaBarrelEndcap ) {
      float Zend = Z_Endcap ;
      if(convEta<0.0 )  Zend = -Zend ;
      float Zlen = Zend - Zvertex ;
      float RR = Zlen/TMath::SinH(convEta); 
      Theta = TMath::ATan(RR/Zend);
      if(Theta<0.0) Theta = Theta+M_PI;
      fEta = -TMath::Log(TMath::Tan(0.5*Theta));              
    } 
    // ---------------------------------------------------

    if(print) {
      std::cout<<"          eta bare = "<<convEta<<std::endl;
      std::cout<<"          eta new  = "<<fEta<<std::endl;
      std::cout<<"          phi      = "<<convPhi<<std::endl;
    }

    convEta = fEta;
    
    Double_t dphi = TMath::Abs(phPhi - convPhi);
    if(dphi > M_PI) dphi = 2.*M_PI-dphi;
    //Double_t deta = c->Eta()-dirconvsc.Eta();
    Double_t deta = TMath::Abs(convEta-phEta);
    Double_t dR = TMath::Sqrt(dphi*dphi+deta*deta);
    //if(dphi < minDphi && TMath::Abs(deta)<minDeta) {
    if(dR < minDR) {
      minDR = dR;
      minDphi = dphi;
      minDeta = TMath::Abs(deta);
      match = c;

      if(print)
        std::cout<<" conv "<<i+1<<" matches with dPhi = "<<minDphi<<"   dEta = "<<minDeta<<std::endl;

    } 
  }
  
  //if(minDphi < dPhiMin && minDeta < dEtaMin)
  if(minDR < dRMin)
    return match;
  else 
    return NULL;
}

bool PhotonTools::PassCiCSelection(Photon* ph, const Vertex* vtx, 
                                   const TrackCol* trackCol,
                                   const ElectronCol* eleCol,
                                   const VertexCol* vtxCol,
                                   double rho, double ptmin, 
                                   bool print, float* kin) {
  
  // these values are taken from the H2GGlobe code... (actually from Marco/s mail)
  float cic4_allcuts_temp_sublead[] = { 
    3.8,         2.2,         1.77,        1.29,
    11.7,        3.4,         3.9,         1.84,
    3.5,         2.2,         2.3,         1.45,
    0.0106,      0.0097,      0.028,       0.027,
    0.082,       0.062,       0.065,       0.048,
    0.94,        0.36,        0.94,        0.32,
    1.,          0.062,       0.97,        0.97,
    1.5,         1.5,         1.5,         1.5 };  // the last line is PixelmatchVeto and un-used
  
  // cut on Et instead of Pt???    
  Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
    
  // compute all relevant observables first
  double ecalIso3 = ph->EcalRecHitIsoDr03();
  double ecalIso4 = ph->EcalRecHitIsoDr04();
  double hcalIso4 = ph->HcalTowerSumEtDr04();

  unsigned int wVtxInd = 0;

  double trackIso1 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);

  // track iso only
  double trackIso3 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);

  // track iso worst vtx
  double trackIso2 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd, vtxCol);
  
  double combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*rho;
  double combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*rho;
  
  double tIso1 = (combIso1) *50./ph->Et();
  double tIso2 = (combIso2) *50./(ph->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
  //double tIso2 = (combIso2) *50./ph->Et();
  double tIso3 = (trackIso3)*50./ph->Et();
  
  double covIEtaIEta  =ph->CoviEtaiEta();
  double HoE = ph->HadOverEm();
  
  double R9 = ph->R9();
  
  double dRTrack = PhotonTools::ElectronVetoCiC(ph, eleCol);

  // check which category it is ...
  int _tCat = 1;
  if ( !isbarrel ) _tCat = 3;
  if ( R9 < 0.94 ) _tCat++;
  
  if(print) {
    std::cout<<" -------------------------- "<<std::endl;
    std::cout<<" photon Et  = "<<ph->Et()<<std::endl;
    std::cout<<"        Eta = "<<ph->SCluster()->Eta()<<std::endl;
    std::cout<<"        HoE = "<<HoE<<std::endl;
    std::cout<<"         R9 = "<<R9<<std::endl;
    std::cout<<"         dR = "<<dRTrack<<std::endl;
    std::cout<<"       iso1 = "<<tIso1<<std::endl; 
    std::cout<<"       iso2 = "<<tIso2<<std::endl; 
    std::cout<<"       iso3 = "<<tIso3<<std::endl; 
  }

  if(kin) {
    kin[0] = tIso1;
    kin[1] = tIso2;
    kin[2] = tIso3;
    kin[3] = covIEtaIEta;
    kin[4] = HoE;
    kin[5] = R9;
    kin[6] = dRTrack;
    kin[7] = (float) ph->Pt();
    kin[8] = (float) ph->Eta();
    kin[9] = (float) ph->Phi();

    // iso quantities separate
    kin[10] = ecalIso3;
    kin[11] = ecalIso4;
    kin[12] = hcalIso4;
    kin[13] = trackIso1;
    kin[14] = trackIso2;
    kin[15] = trackIso3;

    kin[16] = (float) ph->Et();
    kin[17] = (float) ph->E();

  }

  float passCuts = 1.;

  if ( ph->Pt()     <= ptmin      ) passCuts = -1.;

  // not needed anymore, do in pre-selection...
  //if (  ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) passCuts = -1.;
  
  if(   ! (    tIso1                          < cic4_allcuts_temp_sublead[_tCat-1+0*4]
               && tIso2                       < cic4_allcuts_temp_sublead[_tCat-1+1*4]
               && tIso3                       < cic4_allcuts_temp_sublead[_tCat-1+2*4]
               && covIEtaIEta                 < cic4_allcuts_temp_sublead[_tCat-1+3*4]
               && HoE                         < cic4_allcuts_temp_sublead[_tCat-1+4*4]
               && R9                          > cic4_allcuts_temp_sublead[_tCat-1+5*4]
               && dRTrack                     > cic4_allcuts_temp_sublead[_tCat-1+6*4]  )
        ) {
    passCuts = -1.;
  }
  
  if(kin) {    
    kin[18] = passCuts;
    kin[19] = (float) _tCat;
  }    

  if(passCuts > 0.) return true;
  return false;
}
Revision:	1.7
Committed:	Mon Jul 4 13:40:46 2011 UTC (13 years, 10 months ago) by fabstoec
Content type:	text/plain
Branch:	MAIN
Changes since 1.6:	+40 -15 lines
Log Message:	* empty log message *
#	User	Rev	Content
1	fabstoec	1.7	// $Id: PhotonTools.cc,v 1.6 2011/06/29 18:28:05 fabstoec Exp $
2	bendavid	1.1
3			#include "MitPhysics/Utils/interface/PhotonTools.h"
4			#include "MitPhysics/Utils/interface/ElectronTools.h"
5	fabstoec	1.5	#include "MitPhysics/Utils/interface/IsolationTools.h"
6	bendavid	1.1	#include "MitAna/DataTree/interface/StableData.h"
7			#include <TFile.h>
8
9			ClassImp(mithep::PhotonTools)
10
11			using namespace mithep;
12
13			//--------------------------------------------------------------------------------------------------
14			PhotonTools::PhotonTools()
15			{
16			// Constructor.
17			}
18
19			//--------------------------------------------------------------------------------------------------
20			Bool_t PhotonTools::PassConversionId(const Photon p, const DecayParticle c) {
21
22			if (!c) return kTRUE;
23
24			ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
25			Double_t deta = c->Eta()-dirconvsc.Eta();
26			Double_t dphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
27			Double_t eoverp = p->SCluster()->Energy()/c->P();
28
29			if (p->IsEB() && eoverp>2.0) return kFALSE;
30			if (p->IsEE() && eoverp>3.0) return kFALSE;
31
32			if (p->IsEE() && TMath::Abs(deta)>0.01) return kFALSE;
33			if (p->IsEE() && TMath::Abs(dphi)>0.01) return kFALSE;
34
35			return kTRUE;
36
37			}
38
39			//--------------------------------------------------------------------------------------------------
40			Bool_t PhotonTools::PassElectronVeto(const Photon p, const ElectronCol els) {
41
42			Bool_t pass = kTRUE;
43			for (UInt_t i=0; i<els->GetEntries(); ++i) {
44			const Electron *e = els->At(i);
45			if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0) {
46			pass = kFALSE;
47			}
48			}
49
50			return pass;
51			}
52
53			//--------------------------------------------------------------------------------------------------
54	fabstoec	1.4	Double_t PhotonTools::ElectronVetoCiC(const Photon p, const ElectronCol els) {
55
56			for (UInt_t i=0; i<els->GetEntries(); ++i) {
57			const Electron *e = els->At(i);
58	fabstoec	1.5	if (e->SCluster()==p->SCluster() ) {
59			//if( e->GsfTrk()->NExpectedHitsInner()==0 && e->GsfTrk()->Pt() > 2.5 ) {
60			if( e->GsfTrk()->NExpectedHitsInner()==0 ) {
61			double dEta = e->DeltaEtaSuperClusterTrackAtVtx();
62			double dPhi = e->DeltaPhiSuperClusterTrackAtVtx();
63			double dR = TMath::Sqrt(dEtadEta+dPhidPhi);
64			return dR;
65			}
66			}
67			}
68			return 99.;
69	fabstoec	1.4	}
70
71			//--------------------------------------------------------------------------------------------------
72	bendavid	1.1	Bool_t PhotonTools::PassElectronVetoConvRecovery(const Photon p, const ElectronCol els, const DecayParticleCol conversions, const BaseVertex v) {
73
74			Bool_t pass = kTRUE;
75			for (UInt_t i=0; i<els->GetEntries(); ++i) {
76			const Electron *e = els->At(i);
77			if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0 && ElectronTools::PassConversionFilter(e, conversions,
78			v, 0, 1e-6, 2.0, kFALSE, kFALSE) ) {
79			pass = kFALSE;
80			}
81			}
82
83			return pass;
84			}
85
86			//--------------------------------------------------------------------------------------------------
87			Bool_t PhotonTools::PassTriggerMatching(const Photon p, const TriggerObjectCol trigobjs)
88			{
89
90			for (UInt_t i=0; i<trigobjs->GetEntries(); ++i) {
91			const TriggerObject *trigobj = trigobjs->At(i);
92			if (trigobj->TriggerType()==TriggerObject::TriggerCluster \|\| trigobj->TriggerType()==TriggerObject::TriggerElectron \|\| trigobj->TriggerType()==TriggerObject::TriggerPhoton) {
93			if (MathUtils::DeltaR(p->SCluster(),trigobj)<0.3) {
94			return kTRUE;
95			}
96			}
97			}
98
99			return kFALSE;
100
101
102			}
103
104			//--------------------------------------------------------------------------------------------------
105			const DecayParticle PhotonTools::MatchedConversion(const Photon p, const DecayParticleCol *conversions,
106			const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
107			Double_t lxyMin, Double_t dRMin) {
108
109			const DecayParticle *match = 0;
110	bendavid	1.3	Double_t rhosmallest = 999.;
111	bendavid	1.1	for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
112			const DecayParticle *c = conversions->At(i);
113			ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
114			Double_t dr = MathUtils::DeltaR(*c,dirconvsc);
115	bendavid	1.3	Double_t rho = c->Position().Rho();
116			if (dr<dRMin && rho<rhosmallest && c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
117	bendavid	1.1	Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
118			Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
119			if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
120	bendavid	1.3	rhosmallest = rho;
121	bendavid	1.1	match = c;
122			}
123			}
124
125			}
126
127			return match;
128
129			}
130
131	bendavid	1.2	//--------------------------------------------------------------------------------------------------
132			const DecayParticle PhotonTools::MatchedConversion(const Track t, const DecayParticleCol *conversions,
133			const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
134			Double_t lxyMin) {
135
136			for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
137			const DecayParticle *c = conversions->At(i);
138			if (c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
139			Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
140			Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
141			if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
142			const Track ct1 = dynamic_cast<const ChargedParticle>(c->Daughter(0))->Trk();
143			const Track ct2 = dynamic_cast<const ChargedParticle>(c->Daughter(1))->Trk();
144			if (t==ct1 \|\| t==ct2) return c;
145			}
146			}
147
148			}
149
150			return 0;
151
152			}
153
154	bendavid	1.1	PhotonTools::DiphotonR9EtaCats PhotonTools::DiphotonR9EtaCat(const Photon p1, const Photon p2) {
155
156			if (p1->IsEB() && p2->IsEB()) {
157			if (p1->R9()>0.93 && p2->R9()>0.93) return kCat1;
158			else return kCat2;
159
160			}
161			else {
162			if (p1->R9()>0.93 && p2->R9()>0.93) return kCat3;
163			else return kCat4;
164			}
165
166			}
167
168			PhotonTools::DiphotonR9EtaConversionCats PhotonTools::DiphotonR9EtaConversionCat(const Photon p1, const Photon p2, const DecayParticleCol conversions, const BaseVertex v) {
169
170			const DecayParticle *conv1 = MatchedConversion(p1, conversions, v);
171			const DecayParticle *conv2 = MatchedConversion(p2, conversions, v);
172
173			if (p1->IsEB() && p2->IsEB()) {
174			if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat1;
175			else if (conv1\|\|conv2) return kNewCat2;
176			else return kNewCat3;
177
178			}
179			else {
180			if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat4;
181			else if (conv1\|\|conv2) return kNewCat5;
182			else return kNewCat6;
183			}
184
185	fabstoec	1.4	}
186
187			PhotonTools::CiCBaseLineCats PhotonTools::CiCBaseLineCat(const Photon *p) {
188			if( p->IsEB() ) {
189			if ( p->R9() > 0.94 ) return kCiCCat1;
190			else return kCiCCat2;
191			} else {
192			if ( p->R9() > 0.94 ) return kCiCCat3;
193			else return kCiCCat4;
194			}
195			}
196	fabstoec	1.5
197			//--------------------------------------------------------------------------------------------------
198			const DecayParticle PhotonTools::MatchedCiCConversion(const Photon p, const DecayParticleCol *conversions,
199			Double_t dPhiMin,
200	fabstoec	1.7	Double_t dEtaMin,
201			Double_t dRMin,
202			bool print) {
203	fabstoec	1.5
204	fabstoec	1.6	const DecayParticle *match = NULL;
205
206			double minDeta = 999.;
207			double minDphi = 999.;
208	fabstoec	1.7	double minDR = 999.;
209	fabstoec	1.6
210			double phPhi = p->SCluster()->Phi();
211			double phEta = p->SCluster()->Eta();
212	fabstoec	1.7
213			if(print)
214			std::cout<<" --- conv match photon eta = "<<phEta<<" phi = "<<phPhi<<std::endl;
215
216
217	fabstoec	1.5	for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
218			const DecayParticle *c = conversions->At(i);
219	fabstoec	1.6
220	fabstoec	1.7	if(print)
221			std::cout<< " c "<<i+1<<" pt = "<<c->Pt()<<std::endl;
222
223	fabstoec	1.6	if(c->Pt() < 1. ) continue; // is this refittedPirMomentum?
224
225			//ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
226			//ThreeVector dirconvsc = p->CaloPos() - c->Position();
227
228			double convPhi = c->Phi();
229			double convEta = c->Eta();
230	fabstoec	1.7
231	fabstoec	1.6	const ThreeVector wrong(0,0,0);
232			double Zvertex = c->DzCorrected(wrong);
233			// ------------------ FROM GLOBE ----------------------
234			//---Definitions for ECAL
235			const float R_ECAL = 136.5;
236			const float Z_Endcap = 328.0;
237			const float etaBarrelEndcap = 1.479;
238	fabstoec	1.7
239	fabstoec	1.6	//---ETA correction
240			float Theta = 0.0 ;
241			float ZEcal = R_ECAL*sinh(convEta)+Zvertex;
242	fabstoec	1.7
243	fabstoec	1.6	if(ZEcal != 0.0) Theta = TMath::ATan(R_ECAL/ZEcal);
244			if(Theta<0.0) Theta = Theta+M_PI;
245			double fEta = - TMath::Log(TMath::Tan(0.5*Theta));
246	fabstoec	1.5
247	fabstoec	1.6	if( fabs(fEta) > etaBarrelEndcap ) {
248			float Zend = Z_Endcap ;
249			if(convEta<0.0 ) Zend = -Zend ;
250			float Zlen = Zend - Zvertex ;
251			float RR = Zlen/TMath::SinH(convEta);
252			Theta = TMath::ATan(RR/Zend);
253			if(Theta<0.0) Theta = Theta+M_PI;
254			fEta = -TMath::Log(TMath::Tan(0.5*Theta));
255			}
256			// ---------------------------------------------------
257
258	fabstoec	1.7	if(print) {
259			std::cout<<" eta bare = "<<convEta<<std::endl;
260			std::cout<<" eta new = "<<fEta<<std::endl;
261			std::cout<<" phi = "<<convPhi<<std::endl;
262			}
263
264			convEta = fEta;
265
266	fabstoec	1.6	Double_t dphi = TMath::Abs(phPhi - convPhi);
267			if(dphi > M_PI) dphi = 2.*M_PI-dphi;
268			//Double_t deta = c->Eta()-dirconvsc.Eta();
269			Double_t deta = TMath::Abs(convEta-phEta);
270	fabstoec	1.7	Double_t dR = TMath::Sqrt(dphidphi+detadeta);
271			//if(dphi < minDphi && TMath::Abs(deta)<minDeta) {
272			if(dR < minDR) {
273			minDR = dR;
274			minDphi = dphi;
275			minDeta = TMath::Abs(deta);
276			match = c;
277
278			if(print)
279			std::cout<<" conv "<<i+1<<" matches with dPhi = "<<minDphi<<" dEta = "<<minDeta<<std::endl;
280
281	fabstoec	1.6	}
282	fabstoec	1.5	}
283
284	fabstoec	1.7	//if(minDphi < dPhiMin && minDeta < dEtaMin)
285			if(minDR < dRMin)
286			return match;
287			else
288			return NULL;
289	fabstoec	1.5	}
290
291			bool PhotonTools::PassCiCSelection(Photon* ph, const Vertex* vtx,
292			const TrackCol* trackCol,
293			const ElectronCol* eleCol,
294			const VertexCol* vtxCol,
295			double rho, double ptmin,
296			bool print, float* kin) {
297
298			// these values are taken from the H2GGlobe code... (actually from Marco/s mail)
299			float cic4_allcuts_temp_sublead[] = {
300			3.8, 2.2, 1.77, 1.29,
301			11.7, 3.4, 3.9, 1.84,
302			3.5, 2.2, 2.3, 1.45,
303			0.0106, 0.0097, 0.028, 0.027,
304			0.082, 0.062, 0.065, 0.048,
305			0.94, 0.36, 0.94, 0.32,
306			1., 0.062, 0.97, 0.97,
307			1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
308
309			// cut on Et instead of Pt???
310			Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
311
312			// compute all relevant observables first
313			double ecalIso3 = ph->EcalRecHitIsoDr03();
314			double ecalIso4 = ph->EcalRecHitIsoDr04();
315			double hcalIso4 = ph->HcalTowerSumEtDr04();
316
317			unsigned int wVtxInd = 0;
318
319			double trackIso1 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
320
321			// track iso only
322			double trackIso3 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
323
324			// track iso worst vtx
325			double trackIso2 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd, vtxCol);
326
327			double combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*rho;
328			double combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*rho;
329
330			double tIso1 = (combIso1) *50./ph->Et();
331			double tIso2 = (combIso2) *50./(ph->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
332			//double tIso2 = (combIso2) *50./ph->Et();
333			double tIso3 = (trackIso3)*50./ph->Et();
334
335			double covIEtaIEta =ph->CoviEtaiEta();
336			double HoE = ph->HadOverEm();
337
338			double R9 = ph->R9();
339
340			double dRTrack = PhotonTools::ElectronVetoCiC(ph, eleCol);
341
342			// check which category it is ...
343			int _tCat = 1;
344			if ( !isbarrel ) _tCat = 3;
345			if ( R9 < 0.94 ) _tCat++;
346
347			if(print) {
348			std::cout<<" -------------------------- "<<std::endl;
349			std::cout<<" photon Et = "<<ph->Et()<<std::endl;
350			std::cout<<" Eta = "<<ph->SCluster()->Eta()<<std::endl;
351			std::cout<<" HoE = "<<HoE<<std::endl;
352			std::cout<<" R9 = "<<R9<<std::endl;
353			std::cout<<" dR = "<<dRTrack<<std::endl;
354			std::cout<<" iso1 = "<<tIso1<<std::endl;
355			std::cout<<" iso2 = "<<tIso2<<std::endl;
356			std::cout<<" iso3 = "<<tIso3<<std::endl;
357			}
358
359			if(kin) {
360			kin[0] = tIso1;
361			kin[1] = tIso2;
362			kin[2] = tIso3;
363			kin[3] = covIEtaIEta;
364			kin[4] = HoE;
365			kin[5] = R9;
366			kin[6] = dRTrack;
367			kin[7] = (float) ph->Pt();
368			kin[8] = (float) ph->Eta();
369			kin[9] = (float) ph->Phi();
370
371			// iso quantities separate
372			kin[10] = ecalIso3;
373			kin[11] = ecalIso4;
374			kin[12] = hcalIso4;
375			kin[13] = trackIso1;
376			kin[14] = trackIso2;
377			kin[15] = trackIso3;
378
379			kin[16] = (float) ph->Et();
380			kin[17] = (float) ph->E();
381
382			}
383
384			float passCuts = 1.;
385
386			if ( ph->Pt() <= ptmin ) passCuts = -1.;
387
388			// not needed anymore, do in pre-selection...
389			//if ( ph->SCluster()->AbsEta()>=2.5 \|\| (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) passCuts = -1.;
390
391			if( ! ( tIso1 < cic4_allcuts_temp_sublead[_tCat-1+0*4]
392			&& tIso2 < cic4_allcuts_temp_sublead[_tCat-1+1*4]
393			&& tIso3 < cic4_allcuts_temp_sublead[_tCat-1+2*4]
394			&& covIEtaIEta < cic4_allcuts_temp_sublead[_tCat-1+3*4]
395			&& HoE < cic4_allcuts_temp_sublead[_tCat-1+4*4]
396			&& R9 > cic4_allcuts_temp_sublead[_tCat-1+5*4]
397			&& dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4] )
398			) {
399			passCuts = -1.;
400			}
401
402			if(kin) {
403			kin[18] = passCuts;
404			kin[19] = (float) _tCat;
405			}
406
407			if(passCuts > 0.) return true;
408			return false;
409			}