Utils/src/PhotonTools.cc

// $Id: PhotonTools.cc,v 1.9 2011/07/15 17:24:36 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;

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