Utils/src/PhotonTools.cc

// $Id: PhotonTools.cc,v 1.11 2011/07/27 15:17:37 bendavid 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->SCluster()->AbsEta()<1.5 ) {
    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(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();

    kin[18] = (float) ph->SCluster()->Eta();
    kin[19] = (float) ph->SCluster()->Phi();
  }

  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[20] = passCuts;
    kin[21] = (float) _tCat;
  }    

  if(passCuts > 0.) return true;
  return false;
}

const MCParticle *PhotonTools::MatchMC(const Photon *ph, const MCParticleCol *c, Bool_t matchElectrons) {
  
  for (UInt_t i=0; i<c->GetEntries(); ++i) {
    const MCParticle *p = c->At(i);
//     if (p->IsGenerated() && p->AbsPdgId()==11 && p->Mother()->AbsPdgId()==23) {
//       printf("pdgid = %i, status = %i, pt = %5f\n",p->PdgId(),p->Status(),p->Pt());
//     }
    if (matchElectrons && p->AbsPdgId()==11 && p->IsGenerated() && p->Status()==3 && MathUtils::DeltaR(*ph,*p) < 0.3 && p->Mother() && (p->Mother()->AbsPdgId()==23 || p->Mother()->AbsPdgId()==24 || p->Mother()->AbsPdgId()==22)) {
      return p;
    }
    if ( p->AbsPdgId()==22 && p->IsGenerated() && MathUtils::DeltaR(*ph,*p) < 0.3 && p->Mother() && (p->Mother()->AbsPdgId()==25 || p->Mother()->AbsPdgId()<=21) ) {
      return p;
    }
  }
  return 0;
}

PhotonTools::DiphotonR9EtaPtCats PhotonTools::DiphotonR9EtaPtCat(const Photon *p1, const Photon *p2) {

  PhotonTools::CiCBaseLineCats cat1 = CiCBaseLineCat(p1);
  PhotonTools::CiCBaseLineCats cat2 = CiCBaseLineCat(p2);
  
  PhotonTools::DiphotonR9EtaPtCats evtcat = PhotonTools::kOctCat0;
  
  bool ph1IsEB = (cat1 ==  PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
  bool ph2IsEB = (cat2 ==  PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);

  bool ph1IsHR9 = (cat1 ==  PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
  bool ph2IsHR9 = (cat2 ==  PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
  
  if( ph1IsEB && ph2IsEB )
    evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat0 : PhotonTools::kOctCat1);
  else
    evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat2 : PhotonTools::kOctCat3);
  
  float ptgg = (p1->Mom()+p2->Mom()).Pt();
  if (ptgg<40.0) evtcat = PhotonTools::DiphotonR9EtaPtCats(evtcat + 4);
  
  
  return evtcat;
}
Revision:	1.12
Committed:	Wed Aug 3 17:15:44 2011 UTC (13 years, 9 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_024b, Mit_024a, Mit_024
Changes since 1.11:	+34 -15 lines
Log Message:	Reorganize photon code and add new PhotonTreeWriter class
#	User	Rev	Content
1	bendavid	1.12	// $Id: PhotonTools.cc,v 1.11 2011/07/27 15:17:37 bendavid 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	bendavid	1.12	if( p->SCluster()->AbsEta()<1.5 ) {
214	fabstoec	1.4	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	fabstoec	1.9	bool PhotonTools::PassCiCSelection(const Photon* ph, const Vertex* vtx,
323			const TrackCol* trackCol,
324			const ElectronCol* eleCol,
325			const VertexCol* vtxCol,
326	fabstoec	1.5	double rho, double ptmin,
327	fabstoec	1.9	bool applyEleVeto,
328	fabstoec	1.5	bool print, float* kin) {
329	fabstoec	1.9
330	fabstoec	1.5
331			// these values are taken from the H2GGlobe code... (actually from Marco/s mail)
332			float cic4_allcuts_temp_sublead[] = {
333			3.8, 2.2, 1.77, 1.29,
334			11.7, 3.4, 3.9, 1.84,
335			3.5, 2.2, 2.3, 1.45,
336			0.0106, 0.0097, 0.028, 0.027,
337			0.082, 0.062, 0.065, 0.048,
338			0.94, 0.36, 0.94, 0.32,
339			1., 0.062, 0.97, 0.97,
340			1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
341
342			// cut on Et instead of Pt???
343			Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
344
345			// compute all relevant observables first
346			double ecalIso3 = ph->EcalRecHitIsoDr03();
347			double ecalIso4 = ph->EcalRecHitIsoDr04();
348			double hcalIso4 = ph->HcalTowerSumEtDr04();
349
350			unsigned int wVtxInd = 0;
351
352			double trackIso1 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
353
354			// track iso only
355			double trackIso3 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
356
357			// track iso worst vtx
358			double trackIso2 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd, vtxCol);
359
360			double combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*rho;
361			double combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*rho;
362
363			double tIso1 = (combIso1) *50./ph->Et();
364			double tIso2 = (combIso2) *50./(ph->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
365			//double tIso2 = (combIso2) *50./ph->Et();
366			double tIso3 = (trackIso3)*50./ph->Et();
367
368			double covIEtaIEta =ph->CoviEtaiEta();
369			double HoE = ph->HadOverEm();
370
371			double R9 = ph->R9();
372
373			double dRTrack = PhotonTools::ElectronVetoCiC(ph, eleCol);
374
375			// check which category it is ...
376			int _tCat = 1;
377			if ( !isbarrel ) _tCat = 3;
378			if ( R9 < 0.94 ) _tCat++;
379
380			if(print) {
381			std::cout<<" -------------------------- "<<std::endl;
382			std::cout<<" photon Et = "<<ph->Et()<<std::endl;
383			std::cout<<" Eta = "<<ph->SCluster()->Eta()<<std::endl;
384			std::cout<<" HoE = "<<HoE<<std::endl;
385			std::cout<<" R9 = "<<R9<<std::endl;
386			std::cout<<" dR = "<<dRTrack<<std::endl;
387			std::cout<<" iso1 = "<<tIso1<<std::endl;
388			std::cout<<" iso2 = "<<tIso2<<std::endl;
389			std::cout<<" iso3 = "<<tIso3<<std::endl;
390			}
391
392			if(kin) {
393			kin[0] = tIso1;
394			kin[1] = tIso2;
395			kin[2] = tIso3;
396			kin[3] = covIEtaIEta;
397			kin[4] = HoE;
398			kin[5] = R9;
399			kin[6] = dRTrack;
400			kin[7] = (float) ph->Pt();
401			kin[8] = (float) ph->Eta();
402			kin[9] = (float) ph->Phi();
403
404			// iso quantities separate
405			kin[10] = ecalIso3;
406			kin[11] = ecalIso4;
407			kin[12] = hcalIso4;
408			kin[13] = trackIso1;
409			kin[14] = trackIso2;
410			kin[15] = trackIso3;
411
412			kin[16] = (float) ph->Et();
413			kin[17] = (float) ph->E();
414
415	bendavid	1.11	kin[18] = (float) ph->SCluster()->Eta();
416			kin[19] = (float) ph->SCluster()->Phi();
417	fabstoec	1.5	}
418
419			float passCuts = 1.;
420
421			if ( ph->Pt() <= ptmin ) passCuts = -1.;
422
423			// not needed anymore, do in pre-selection...
424			//if ( ph->SCluster()->AbsEta()>=2.5 \|\| (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) passCuts = -1.;
425
426			if( ! ( tIso1 < cic4_allcuts_temp_sublead[_tCat-1+0*4]
427			&& tIso2 < cic4_allcuts_temp_sublead[_tCat-1+1*4]
428			&& tIso3 < cic4_allcuts_temp_sublead[_tCat-1+2*4]
429			&& covIEtaIEta < cic4_allcuts_temp_sublead[_tCat-1+3*4]
430			&& HoE < cic4_allcuts_temp_sublead[_tCat-1+4*4]
431			&& R9 > cic4_allcuts_temp_sublead[_tCat-1+5*4]
432	fabstoec	1.9	&& ( dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4] \|\| !applyEleVeto ) ) ) passCuts = -1.;
433	fabstoec	1.5
434			if(kin) {
435	bendavid	1.11	kin[20] = passCuts;
436			kin[21] = (float) _tCat;
437	fabstoec	1.5	}
438
439			if(passCuts > 0.) return true;
440			return false;
441			}
442	bendavid	1.11
443	bendavid	1.12	const MCParticle PhotonTools::MatchMC(const Photon ph, const MCParticleCol *c, Bool_t matchElectrons) {
444	bendavid	1.11
445			for (UInt_t i=0; i<c->GetEntries(); ++i) {
446			const MCParticle *p = c->At(i);
447	bendavid	1.12	// if (p->IsGenerated() && p->AbsPdgId()==11 && p->Mother()->AbsPdgId()==23) {
448			// printf("pdgid = %i, status = %i, pt = %5f\n",p->PdgId(),p->Status(),p->Pt());
449			// }
450			if (matchElectrons && p->AbsPdgId()==11 && p->IsGenerated() && p->Status()==3 && MathUtils::DeltaR(ph,p) < 0.3 && p->Mother() && (p->Mother()->AbsPdgId()==23 \|\| p->Mother()->AbsPdgId()==24 \|\| p->Mother()->AbsPdgId()==22)) {
451			return p;
452			}
453	bendavid	1.11	if ( p->AbsPdgId()==22 && p->IsGenerated() && MathUtils::DeltaR(ph,p) < 0.3 && p->Mother() && (p->Mother()->AbsPdgId()==25 \|\| p->Mother()->AbsPdgId()<=21) ) {
454			return p;
455			}
456			}
457			return 0;
458	bendavid	1.12	}
459
460			PhotonTools::DiphotonR9EtaPtCats PhotonTools::DiphotonR9EtaPtCat(const Photon p1, const Photon p2) {
461
462			PhotonTools::CiCBaseLineCats cat1 = CiCBaseLineCat(p1);
463			PhotonTools::CiCBaseLineCats cat2 = CiCBaseLineCat(p2);
464
465			PhotonTools::DiphotonR9EtaPtCats evtcat = PhotonTools::kOctCat0;
466
467			bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 \|\| cat1 == PhotonTools::kCiCCat2);
468			bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 \|\| cat2 == PhotonTools::kCiCCat2);
469
470			bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 \|\| cat1 == PhotonTools::kCiCCat3);
471			bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 \|\| cat2 == PhotonTools::kCiCCat3);
472
473			if( ph1IsEB && ph2IsEB )
474			evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat0 : PhotonTools::kOctCat1);
475			else
476			evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat2 : PhotonTools::kOctCat3);
477
478			float ptgg = (p1->Mom()+p2->Mom()).Pt();
479			if (ptgg<40.0) evtcat = PhotonTools::DiphotonR9EtaPtCats(evtcat + 4);
480
481
482			return evtcat;
483	bendavid	1.11	}