Utils/src/PhotonTools.cc

// $Id: PhotonTools.cc,v 1.7 2011/07/04 13:40:46 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>
#include <TRandom3.h>

ClassImp(mithep::PhotonTools)

using namespace mithep;

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


void PhotonTools::ScalePhoton(Photon* p, Double_t scale) {
  if( !p ) return;
  FourVectorM mom = p->Mom();
  p->SetMom(scale*mom.X(), scale*mom.Y(), scale*mom.Z(), scale*mom.E());
  
}

void PhotonTools::SmearPhoton(Photon* p, TRandom3* rng, Double_t width, UInt_t iSeed) {
  
  if( !p  )       return;
  if( !rng)       return;
  if( width < 0.) return;

  if( iSeed ) rng->SetSeed(iSeed);
  FourVectorM mom = p->Mom();
  Double_t scale = rng->Gaus(1.,width);
  if( scale > 0)
    p->SetMom(scale*mom.X(), scale*mom.Y(), scale*mom.Z(), scale*mom.E());

  return;
}

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

  // shoud NEVER happen, but you never know...
  return kCiCNoCat;
}

//--------------------------------------------------------------------------------------------------
const DecayParticle *PhotonTools::MatchedCiCConversion(const Photon *p, const DecayParticleCol *conversions, 
                                                       Double_t dPhiMin,
                                                       Double_t dEtaMin,
                                                       Double_t dRMin,
                                                       bool print) {
  
  // if there are no conversons, return
  if ( !p || !conversions)  return NULL;

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