Utils/src/PhotonTools.cc

// $Id: PhotonTools.cc,v 1.8 2011/07/08 17:54:27 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 applyEleVeto, 
                                   bool print, float* kin) {

  const Photon* phpass = ph;
  return PhotonTools::PassCiCSelection(phpass, vtx, trackCol, eleCol, vtxCol, rho, ptmin, applyEleVeto, print, kin);
}

bool PhotonTools::PassCiCSelection(const Photon* ph, const Vertex* vtx, 
                                   const TrackCol* trackCol,
                                   const ElectronCol* eleCol,
                                   const VertexCol* vtxCol,
                                   double rho, double ptmin, 
                                   bool applyEleVeto,
                                   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] || !applyEleVeto ) ) )   passCuts = -1.;
  
  if(kin) {    
    kin[18] = passCuts;
    kin[19] = (float) _tCat;
  }    

  if(passCuts > 0.) return true;
  return false;
}
Revision:	1.9
Committed:	Fri Jul 15 17:24:36 2011 UTC (13 years, 9 months ago) by fabstoec
Content type:	text/plain
Branch:	MAIN
Changes since 1.8:	+16 -5 lines
Log Message:	added/improved stuff for CiC Selection
#	User	Rev	Content
1	fabstoec	1.9	// $Id: PhotonTools.cc,v 1.8 2011/07/08 17:54:27 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	fabstoec	1.9	double rho, double ptmin,
326			bool applyEleVeto,
327			bool print, float* kin) {
328
329			const Photon* phpass = ph;
330			return PhotonTools::PassCiCSelection(phpass, vtx, trackCol, eleCol, vtxCol, rho, ptmin, applyEleVeto, print, kin);
331			}
332
333			bool PhotonTools::PassCiCSelection(const Photon* ph, const Vertex* vtx,
334			const TrackCol* trackCol,
335			const ElectronCol* eleCol,
336			const VertexCol* vtxCol,
337	fabstoec	1.5	double rho, double ptmin,
338	fabstoec	1.9	bool applyEleVeto,
339	fabstoec	1.5	bool print, float* kin) {
340	fabstoec	1.9
341	fabstoec	1.5
342			// these values are taken from the H2GGlobe code... (actually from Marco/s mail)
343			float cic4_allcuts_temp_sublead[] = {
344			3.8, 2.2, 1.77, 1.29,
345			11.7, 3.4, 3.9, 1.84,
346			3.5, 2.2, 2.3, 1.45,
347			0.0106, 0.0097, 0.028, 0.027,
348			0.082, 0.062, 0.065, 0.048,
349			0.94, 0.36, 0.94, 0.32,
350			1., 0.062, 0.97, 0.97,
351			1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
352
353			// cut on Et instead of Pt???
354			Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
355
356			// compute all relevant observables first
357			double ecalIso3 = ph->EcalRecHitIsoDr03();
358			double ecalIso4 = ph->EcalRecHitIsoDr04();
359			double hcalIso4 = ph->HcalTowerSumEtDr04();
360
361			unsigned int wVtxInd = 0;
362
363			double trackIso1 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
364
365			// track iso only
366			double trackIso3 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
367
368			// track iso worst vtx
369			double trackIso2 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd, vtxCol);
370
371			double combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*rho;
372			double combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*rho;
373
374			double tIso1 = (combIso1) *50./ph->Et();
375			double tIso2 = (combIso2) *50./(ph->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
376			//double tIso2 = (combIso2) *50./ph->Et();
377			double tIso3 = (trackIso3)*50./ph->Et();
378
379			double covIEtaIEta =ph->CoviEtaiEta();
380			double HoE = ph->HadOverEm();
381
382			double R9 = ph->R9();
383
384			double dRTrack = PhotonTools::ElectronVetoCiC(ph, eleCol);
385
386			// check which category it is ...
387			int _tCat = 1;
388			if ( !isbarrel ) _tCat = 3;
389			if ( R9 < 0.94 ) _tCat++;
390
391			if(print) {
392			std::cout<<" -------------------------- "<<std::endl;
393			std::cout<<" photon Et = "<<ph->Et()<<std::endl;
394			std::cout<<" Eta = "<<ph->SCluster()->Eta()<<std::endl;
395			std::cout<<" HoE = "<<HoE<<std::endl;
396			std::cout<<" R9 = "<<R9<<std::endl;
397			std::cout<<" dR = "<<dRTrack<<std::endl;
398			std::cout<<" iso1 = "<<tIso1<<std::endl;
399			std::cout<<" iso2 = "<<tIso2<<std::endl;
400			std::cout<<" iso3 = "<<tIso3<<std::endl;
401			}
402
403			if(kin) {
404			kin[0] = tIso1;
405			kin[1] = tIso2;
406			kin[2] = tIso3;
407			kin[3] = covIEtaIEta;
408			kin[4] = HoE;
409			kin[5] = R9;
410			kin[6] = dRTrack;
411			kin[7] = (float) ph->Pt();
412			kin[8] = (float) ph->Eta();
413			kin[9] = (float) ph->Phi();
414
415			// iso quantities separate
416			kin[10] = ecalIso3;
417			kin[11] = ecalIso4;
418			kin[12] = hcalIso4;
419			kin[13] = trackIso1;
420			kin[14] = trackIso2;
421			kin[15] = trackIso3;
422
423			kin[16] = (float) ph->Et();
424			kin[17] = (float) ph->E();
425
426			}
427
428			float passCuts = 1.;
429
430			if ( ph->Pt() <= ptmin ) passCuts = -1.;
431
432			// not needed anymore, do in pre-selection...
433			//if ( ph->SCluster()->AbsEta()>=2.5 \|\| (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) passCuts = -1.;
434
435			if( ! ( tIso1 < cic4_allcuts_temp_sublead[_tCat-1+0*4]
436			&& tIso2 < cic4_allcuts_temp_sublead[_tCat-1+1*4]
437			&& tIso3 < cic4_allcuts_temp_sublead[_tCat-1+2*4]
438			&& covIEtaIEta < cic4_allcuts_temp_sublead[_tCat-1+3*4]
439			&& HoE < cic4_allcuts_temp_sublead[_tCat-1+4*4]
440			&& R9 > cic4_allcuts_temp_sublead[_tCat-1+5*4]
441	fabstoec	1.9	&& ( dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4] \|\| !applyEleVeto ) ) ) passCuts = -1.;
442	fabstoec	1.5
443			if(kin) {
444			kin[18] = passCuts;
445			kin[19] = (float) _tCat;
446			}
447
448			if(passCuts > 0.) return true;
449			return false;
450			}