Utils/src/PhotonTools.cc

// $Id: PhotonTools.cc,v 1.14 2011/10/07 09:56: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;
}

void PhotonTools::SmearPhotonError(Photon* p, Double_t width) {
  
  if( !p  )       return;
  if( width < 0.) return;  
  
  
  Double_t smear = p->EnergySmearing();
  p->SetEnergySmearing(TMath::Sqrt(smear*smear+width*width*p->E()*p->E()));
  
  Double_t err = p->EnergyErrSmeared();
  if (err>=0.0) {
    p->SetEnergyErrSmeared(TMath::Sqrt(err*err + width*width*p->E()*p->E()));
  }
  
}

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

//--------------------------------------------------------------------------------------------------
const Electron *PhotonTools::MatchedElectron(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() ) {
      return e;
    }
  }
  
  return 0;
}

//--------------------------------------------------------------------------------------------------
const Photon *PhotonTools::MatchedPhoton(const Electron *e, const PhotonCol *phs) {

  for (UInt_t i=0; i<phs->GetEntries(); ++i) {
    const Photon *p = phs->At(i);
    if ( p->SCluster()==e->SCluster() ) {
      return p;
    }
  }
  
  return 0;
}

//--------------------------------------------------------------------------------------------------
const SuperCluster *PhotonTools::MatchedSC(const SuperCluster *psc, const SuperClusterCol *scs, Double_t drMin) {

  Double_t drsmallest = 999.;
  const SuperCluster *match = 0;
  for (UInt_t i=0; i<scs->GetEntries(); ++i) {
    const SuperCluster *sc = scs->At(i);
    Double_t dr = MathUtils::DeltaR(*sc,*psc);
    if ( dr<drsmallest && dr<drMin ) {
      drsmallest = dr;
      match = sc;
    }
  }
  
  return match;
}

//--------------------------------------------------------------------------------------------------
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) {
  
  return MatchedConversion(p->SCluster(), conversions, vtx, nWrongHitsMax, probMin, lxyMin, dRMin);
  
}

//--------------------------------------------------------------------------------------------------
const DecayParticle *PhotonTools::MatchedConversion(const SuperCluster *sc, 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(sc->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<<" trackIso1  = "<<trackIso1<<std::endl;
    std::cout<<" trackIso2  = "<<trackIso2<<std::endl;
    std::cout<<" trackIso3  = "<<trackIso3<<std::endl;
    std::cout<<" ecalIso3   = "<<ecalIso3<<std::endl;
    std::cout<<" ecalIso4   = "<<ecalIso4<<std::endl;
    std::cout<<" hcalIso4   = "<<hcalIso4<<std::endl;

    std::cout<<" photon Et  = "<<ph->Et()<<std::endl;
    std::cout<<"        Eta = "<<ph->SCluster()->Eta()<<std::endl;
    std::cout<<"        HoE = "<<HoE<<"   ("<<cic4_allcuts_temp_sublead[_tCat-1+4*4]<<")"<<std::endl;
    std::cout<<"         R9 = "<<R9<<"   ("<<cic4_allcuts_temp_sublead[_tCat-1+5*4]<<")"<<std::endl;
    std::cout<<"         dR = "<<dRTrack<<"   ("<<cic4_allcuts_temp_sublead[_tCat-1+6*4]<<")"<<std::endl;
    std::cout<<"       iso1 = "<<tIso1<<"   ("<<cic4_allcuts_temp_sublead[_tCat-1+0*4]<<")"<<std::endl; 
    std::cout<<"       iso2 = "<<tIso2<<"   ("<<cic4_allcuts_temp_sublead[_tCat-1+1*4]<<")"<<std::endl; 
    std::cout<<"       iso3 = "<<tIso3<<"   ("<<cic4_allcuts_temp_sublead[_tCat-1+2*4]<<")"<<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(print) std::cout<<"   ---> "<<passCuts<<std::endl;

  if(kin) {    
    kin[20] = passCuts;
    kin[21] = (float) _tCat;
  }    

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

const MCParticle *PhotonTools::MatchMC(const Particle *ph, const MCParticleCol *c, Bool_t matchElectrons) {

//  printf("Start loop\n");
  for (UInt_t i=0; i<c->GetEntries(); ++i) {
    const MCParticle *p = c->At(i);
//     if (p->IsGenerated() && p->AbsPdgId()==11 && (p->DistinctMother()->AbsPdgId()==23|| p->DistinctMother()->AbsPdgId()==24)) {
//       printf("pdgid = %i, status = %i, pt = %5f, mother pdgid = %i\n",p->PdgId(),p->Status(),p->Pt(),p->Mother()->PdgId());
//     }
    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.15
Committed:	Tue Oct 18 11:27:19 2011 UTC (13 years, 6 months ago) by fabstoec
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_025c, Mit_025b, Mit_025a
Changes since 1.14:	+18 -8 lines
Log Message:	Added electron Veto inversion to CiCTrackIso & MVA Tools
#	User	Rev	Content
1	fabstoec	1.15	// $Id: PhotonTools.cc,v 1.14 2011/10/07 09:56: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	bendavid	1.13
41			return;
42			}
43	fabstoec	1.8
44	bendavid	1.13	void PhotonTools::SmearPhotonError(Photon* p, Double_t width) {
45
46			if( !p ) return;
47			if( width < 0.) return;
48
49	bendavid	1.14
50			Double_t smear = p->EnergySmearing();
51			p->SetEnergySmearing(TMath::Sqrt(smearsmear+widthwidthp->E()p->E()));
52
53	bendavid	1.13	Double_t err = p->EnergyErrSmeared();
54			if (err>=0.0) {
55			p->SetEnergyErrSmeared(TMath::Sqrt(errerr + widthwidthp->E()p->E()));
56			}
57
58	fabstoec	1.8	}
59
60	bendavid	1.1	//--------------------------------------------------------------------------------------------------
61			Bool_t PhotonTools::PassConversionId(const Photon p, const DecayParticle c) {
62
63			if (!c) return kTRUE;
64
65			ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
66			Double_t deta = c->Eta()-dirconvsc.Eta();
67			Double_t dphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
68			Double_t eoverp = p->SCluster()->Energy()/c->P();
69
70			if (p->IsEB() && eoverp>2.0) return kFALSE;
71			if (p->IsEE() && eoverp>3.0) return kFALSE;
72
73			if (p->IsEE() && TMath::Abs(deta)>0.01) return kFALSE;
74			if (p->IsEE() && TMath::Abs(dphi)>0.01) return kFALSE;
75
76			return kTRUE;
77
78			}
79
80			//--------------------------------------------------------------------------------------------------
81			Bool_t PhotonTools::PassElectronVeto(const Photon p, const ElectronCol els) {
82
83			Bool_t pass = kTRUE;
84			for (UInt_t i=0; i<els->GetEntries(); ++i) {
85			const Electron *e = els->At(i);
86			if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0) {
87			pass = kFALSE;
88			}
89			}
90
91			return pass;
92			}
93
94			//--------------------------------------------------------------------------------------------------
95	bendavid	1.13	const Electron PhotonTools::MatchedElectron(const Photon p, const ElectronCol *els) {
96
97			for (UInt_t i=0; i<els->GetEntries(); ++i) {
98			const Electron *e = els->At(i);
99			if ( e->SCluster()==p->SCluster() ) {
100			return e;
101			}
102			}
103
104			return 0;
105			}
106
107			//--------------------------------------------------------------------------------------------------
108			const Photon PhotonTools::MatchedPhoton(const Electron e, const PhotonCol *phs) {
109
110			for (UInt_t i=0; i<phs->GetEntries(); ++i) {
111			const Photon *p = phs->At(i);
112			if ( p->SCluster()==e->SCluster() ) {
113			return p;
114			}
115			}
116
117			return 0;
118			}
119
120			//--------------------------------------------------------------------------------------------------
121			const SuperCluster PhotonTools::MatchedSC(const SuperCluster psc, const SuperClusterCol *scs, Double_t drMin) {
122
123			Double_t drsmallest = 999.;
124			const SuperCluster *match = 0;
125			for (UInt_t i=0; i<scs->GetEntries(); ++i) {
126			const SuperCluster *sc = scs->At(i);
127			Double_t dr = MathUtils::DeltaR(sc,psc);
128			if ( dr<drsmallest && dr<drMin ) {
129			drsmallest = dr;
130			match = sc;
131			}
132			}
133
134			return match;
135			}
136
137			//--------------------------------------------------------------------------------------------------
138	fabstoec	1.4	Double_t PhotonTools::ElectronVetoCiC(const Photon p, const ElectronCol els) {
139
140			for (UInt_t i=0; i<els->GetEntries(); ++i) {
141			const Electron *e = els->At(i);
142	fabstoec	1.5	if (e->SCluster()==p->SCluster() ) {
143			//if( e->GsfTrk()->NExpectedHitsInner()==0 && e->GsfTrk()->Pt() > 2.5 ) {
144			if( e->GsfTrk()->NExpectedHitsInner()==0 ) {
145			double dEta = e->DeltaEtaSuperClusterTrackAtVtx();
146			double dPhi = e->DeltaPhiSuperClusterTrackAtVtx();
147			double dR = TMath::Sqrt(dEtadEta+dPhidPhi);
148			return dR;
149			}
150			}
151			}
152			return 99.;
153	fabstoec	1.4	}
154
155			//--------------------------------------------------------------------------------------------------
156	bendavid	1.1	Bool_t PhotonTools::PassElectronVetoConvRecovery(const Photon p, const ElectronCol els, const DecayParticleCol conversions, const BaseVertex v) {
157
158			Bool_t pass = kTRUE;
159			for (UInt_t i=0; i<els->GetEntries(); ++i) {
160			const Electron *e = els->At(i);
161			if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0 && ElectronTools::PassConversionFilter(e, conversions,
162			v, 0, 1e-6, 2.0, kFALSE, kFALSE) ) {
163			pass = kFALSE;
164			}
165			}
166
167			return pass;
168			}
169
170			//--------------------------------------------------------------------------------------------------
171			Bool_t PhotonTools::PassTriggerMatching(const Photon p, const TriggerObjectCol trigobjs)
172			{
173
174			for (UInt_t i=0; i<trigobjs->GetEntries(); ++i) {
175			const TriggerObject *trigobj = trigobjs->At(i);
176			if (trigobj->TriggerType()==TriggerObject::TriggerCluster \|\| trigobj->TriggerType()==TriggerObject::TriggerElectron \|\| trigobj->TriggerType()==TriggerObject::TriggerPhoton) {
177			if (MathUtils::DeltaR(p->SCluster(),trigobj)<0.3) {
178			return kTRUE;
179			}
180			}
181			}
182
183			return kFALSE;
184
185
186			}
187
188			//--------------------------------------------------------------------------------------------------
189			const DecayParticle PhotonTools::MatchedConversion(const Photon p, const DecayParticleCol *conversions,
190			const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
191			Double_t lxyMin, Double_t dRMin) {
192
193	bendavid	1.13	return MatchedConversion(p->SCluster(), conversions, vtx, nWrongHitsMax, probMin, lxyMin, dRMin);
194
195			}
196
197			//--------------------------------------------------------------------------------------------------
198			const DecayParticle PhotonTools::MatchedConversion(const SuperCluster sc, const DecayParticleCol *conversions,
199			const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
200			Double_t lxyMin, Double_t dRMin) {
201
202	bendavid	1.1	const DecayParticle *match = 0;
203	bendavid	1.3	Double_t rhosmallest = 999.;
204	bendavid	1.1	for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
205			const DecayParticle *c = conversions->At(i);
206	bendavid	1.13	ThreeVector dirconvsc = ThreeVector(sc->Point()) - c->Position();
207	bendavid	1.1	Double_t dr = MathUtils::DeltaR(*c,dirconvsc);
208	bendavid	1.3	Double_t rho = c->Position().Rho();
209			if (dr<dRMin && rho<rhosmallest && c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
210	bendavid	1.1	Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
211			Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
212			if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
213	bendavid	1.3	rhosmallest = rho;
214	bendavid	1.1	match = c;
215			}
216			}
217
218			}
219
220			return match;
221
222			}
223
224	bendavid	1.2	//--------------------------------------------------------------------------------------------------
225			const DecayParticle PhotonTools::MatchedConversion(const Track t, const DecayParticleCol *conversions,
226			const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
227			Double_t lxyMin) {
228
229			for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
230			const DecayParticle *c = conversions->At(i);
231			if (c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
232			Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
233			Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
234			if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
235			const Track ct1 = dynamic_cast<const ChargedParticle>(c->Daughter(0))->Trk();
236			const Track ct2 = dynamic_cast<const ChargedParticle>(c->Daughter(1))->Trk();
237			if (t==ct1 \|\| t==ct2) return c;
238			}
239			}
240
241			}
242
243			return 0;
244
245			}
246
247	bendavid	1.1	PhotonTools::DiphotonR9EtaCats PhotonTools::DiphotonR9EtaCat(const Photon p1, const Photon p2) {
248
249			if (p1->IsEB() && p2->IsEB()) {
250			if (p1->R9()>0.93 && p2->R9()>0.93) return kCat1;
251			else return kCat2;
252
253			}
254			else {
255			if (p1->R9()>0.93 && p2->R9()>0.93) return kCat3;
256			else return kCat4;
257			}
258
259			}
260
261			PhotonTools::DiphotonR9EtaConversionCats PhotonTools::DiphotonR9EtaConversionCat(const Photon p1, const Photon p2, const DecayParticleCol conversions, const BaseVertex v) {
262
263			const DecayParticle *conv1 = MatchedConversion(p1, conversions, v);
264			const DecayParticle *conv2 = MatchedConversion(p2, conversions, v);
265
266			if (p1->IsEB() && p2->IsEB()) {
267			if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat1;
268			else if (conv1\|\|conv2) return kNewCat2;
269			else return kNewCat3;
270
271			}
272			else {
273			if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat4;
274			else if (conv1\|\|conv2) return kNewCat5;
275			else return kNewCat6;
276			}
277
278	fabstoec	1.4	}
279
280			PhotonTools::CiCBaseLineCats PhotonTools::CiCBaseLineCat(const Photon *p) {
281	bendavid	1.12	if( p->SCluster()->AbsEta()<1.5 ) {
282	fabstoec	1.4	if ( p->R9() > 0.94 ) return kCiCCat1;
283			else return kCiCCat2;
284			} else {
285			if ( p->R9() > 0.94 ) return kCiCCat3;
286			else return kCiCCat4;
287			}
288	fabstoec	1.8
289			// shoud NEVER happen, but you never know...
290			return kCiCNoCat;
291	fabstoec	1.4	}
292	fabstoec	1.5
293			//--------------------------------------------------------------------------------------------------
294			const DecayParticle PhotonTools::MatchedCiCConversion(const Photon p, const DecayParticleCol *conversions,
295			Double_t dPhiMin,
296	fabstoec	1.7	Double_t dEtaMin,
297			Double_t dRMin,
298			bool print) {
299	fabstoec	1.5
300	fabstoec	1.8	// if there are no conversons, return
301			if ( !p \|\| !conversions) return NULL;
302
303	fabstoec	1.6	const DecayParticle *match = NULL;
304
305			double minDeta = 999.;
306			double minDphi = 999.;
307	fabstoec	1.7	double minDR = 999.;
308	fabstoec	1.6
309			double phPhi = p->SCluster()->Phi();
310			double phEta = p->SCluster()->Eta();
311	fabstoec	1.7
312			if(print)
313			std::cout<<" --- conv match photon eta = "<<phEta<<" phi = "<<phPhi<<std::endl;
314
315
316	fabstoec	1.5	for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
317			const DecayParticle *c = conversions->At(i);
318	fabstoec	1.6
319	fabstoec	1.7	if(print)
320			std::cout<< " c "<<i+1<<" pt = "<<c->Pt()<<std::endl;
321
322	fabstoec	1.6	if(c->Pt() < 1. ) continue; // is this refittedPirMomentum?
323
324			//ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
325			//ThreeVector dirconvsc = p->CaloPos() - c->Position();
326
327			double convPhi = c->Phi();
328			double convEta = c->Eta();
329	fabstoec	1.7
330	fabstoec	1.6	const ThreeVector wrong(0,0,0);
331			double Zvertex = c->DzCorrected(wrong);
332			// ------------------ FROM GLOBE ----------------------
333			//---Definitions for ECAL
334			const float R_ECAL = 136.5;
335			const float Z_Endcap = 328.0;
336			const float etaBarrelEndcap = 1.479;
337	fabstoec	1.7
338	fabstoec	1.6	//---ETA correction
339			float Theta = 0.0 ;
340			float ZEcal = R_ECAL*sinh(convEta)+Zvertex;
341	fabstoec	1.7
342	fabstoec	1.6	if(ZEcal != 0.0) Theta = TMath::ATan(R_ECAL/ZEcal);
343			if(Theta<0.0) Theta = Theta+M_PI;
344			double fEta = - TMath::Log(TMath::Tan(0.5*Theta));
345	fabstoec	1.5
346	fabstoec	1.6	if( fabs(fEta) > etaBarrelEndcap ) {
347			float Zend = Z_Endcap ;
348			if(convEta<0.0 ) Zend = -Zend ;
349			float Zlen = Zend - Zvertex ;
350			float RR = Zlen/TMath::SinH(convEta);
351			Theta = TMath::ATan(RR/Zend);
352			if(Theta<0.0) Theta = Theta+M_PI;
353			fEta = -TMath::Log(TMath::Tan(0.5*Theta));
354			}
355			// ---------------------------------------------------
356
357	fabstoec	1.7	if(print) {
358			std::cout<<" eta bare = "<<convEta<<std::endl;
359			std::cout<<" eta new = "<<fEta<<std::endl;
360			std::cout<<" phi = "<<convPhi<<std::endl;
361			}
362
363			convEta = fEta;
364
365	fabstoec	1.6	Double_t dphi = TMath::Abs(phPhi - convPhi);
366			if(dphi > M_PI) dphi = 2.*M_PI-dphi;
367			//Double_t deta = c->Eta()-dirconvsc.Eta();
368			Double_t deta = TMath::Abs(convEta-phEta);
369	fabstoec	1.7	Double_t dR = TMath::Sqrt(dphidphi+detadeta);
370			//if(dphi < minDphi && TMath::Abs(deta)<minDeta) {
371			if(dR < minDR) {
372			minDR = dR;
373			minDphi = dphi;
374			minDeta = TMath::Abs(deta);
375			match = c;
376
377			if(print)
378			std::cout<<" conv "<<i+1<<" matches with dPhi = "<<minDphi<<" dEta = "<<minDeta<<std::endl;
379
380	fabstoec	1.6	}
381	fabstoec	1.5	}
382
383	fabstoec	1.7	//if(minDphi < dPhiMin && minDeta < dEtaMin)
384			if(minDR < dRMin)
385			return match;
386			else
387			return NULL;
388	fabstoec	1.5	}
389
390	fabstoec	1.9	bool PhotonTools::PassCiCSelection(const Photon* ph, const Vertex* vtx,
391			const TrackCol* trackCol,
392			const ElectronCol* eleCol,
393			const VertexCol* vtxCol,
394	fabstoec	1.5	double rho, double ptmin,
395	fabstoec	1.9	bool applyEleVeto,
396	fabstoec	1.5	bool print, float* kin) {
397	fabstoec	1.9
398	fabstoec	1.5
399			// these values are taken from the H2GGlobe code... (actually from Marco/s mail)
400			float cic4_allcuts_temp_sublead[] = {
401			3.8, 2.2, 1.77, 1.29,
402			11.7, 3.4, 3.9, 1.84,
403			3.5, 2.2, 2.3, 1.45,
404			0.0106, 0.0097, 0.028, 0.027,
405			0.082, 0.062, 0.065, 0.048,
406			0.94, 0.36, 0.94, 0.32,
407			1., 0.062, 0.97, 0.97,
408			1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
409
410			// cut on Et instead of Pt???
411			Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
412
413			// compute all relevant observables first
414			double ecalIso3 = ph->EcalRecHitIsoDr03();
415			double ecalIso4 = ph->EcalRecHitIsoDr04();
416			double hcalIso4 = ph->HcalTowerSumEtDr04();
417
418			unsigned int wVtxInd = 0;
419
420			double trackIso1 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
421
422			// track iso only
423			double trackIso3 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
424
425			// track iso worst vtx
426			double trackIso2 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd, vtxCol);
427
428			double combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*rho;
429			double combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*rho;
430
431			double tIso1 = (combIso1) *50./ph->Et();
432			double tIso2 = (combIso2) *50./(ph->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
433			//double tIso2 = (combIso2) *50./ph->Et();
434			double tIso3 = (trackIso3)*50./ph->Et();
435
436			double covIEtaIEta =ph->CoviEtaiEta();
437			double HoE = ph->HadOverEm();
438
439			double R9 = ph->R9();
440
441			double dRTrack = PhotonTools::ElectronVetoCiC(ph, eleCol);
442
443			// check which category it is ...
444			int _tCat = 1;
445			if ( !isbarrel ) _tCat = 3;
446			if ( R9 < 0.94 ) _tCat++;
447
448			if(print) {
449			std::cout<<" -------------------------- "<<std::endl;
450	fabstoec	1.15
451			std::cout<<" trackIso1 = "<<trackIso1<<std::endl;
452			std::cout<<" trackIso2 = "<<trackIso2<<std::endl;
453			std::cout<<" trackIso3 = "<<trackIso3<<std::endl;
454			std::cout<<" ecalIso3 = "<<ecalIso3<<std::endl;
455			std::cout<<" ecalIso4 = "<<ecalIso4<<std::endl;
456			std::cout<<" hcalIso4 = "<<hcalIso4<<std::endl;
457
458	fabstoec	1.5	std::cout<<" photon Et = "<<ph->Et()<<std::endl;
459			std::cout<<" Eta = "<<ph->SCluster()->Eta()<<std::endl;
460	fabstoec	1.15	std::cout<<" HoE = "<<HoE<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+4*4]<<")"<<std::endl;
461			std::cout<<" R9 = "<<R9<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+5*4]<<")"<<std::endl;
462			std::cout<<" dR = "<<dRTrack<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+6*4]<<")"<<std::endl;
463			std::cout<<" iso1 = "<<tIso1<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+0*4]<<")"<<std::endl;
464			std::cout<<" iso2 = "<<tIso2<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+1*4]<<")"<<std::endl;
465			std::cout<<" iso3 = "<<tIso3<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+2*4]<<")"<<std::endl;
466	fabstoec	1.5	}
467
468			if(kin) {
469			kin[0] = tIso1;
470			kin[1] = tIso2;
471			kin[2] = tIso3;
472			kin[3] = covIEtaIEta;
473			kin[4] = HoE;
474			kin[5] = R9;
475			kin[6] = dRTrack;
476			kin[7] = (float) ph->Pt();
477			kin[8] = (float) ph->Eta();
478			kin[9] = (float) ph->Phi();
479
480			// iso quantities separate
481			kin[10] = ecalIso3;
482			kin[11] = ecalIso4;
483			kin[12] = hcalIso4;
484			kin[13] = trackIso1;
485			kin[14] = trackIso2;
486			kin[15] = trackIso3;
487
488			kin[16] = (float) ph->Et();
489			kin[17] = (float) ph->E();
490
491	bendavid	1.11	kin[18] = (float) ph->SCluster()->Eta();
492			kin[19] = (float) ph->SCluster()->Phi();
493	fabstoec	1.5	}
494
495			float passCuts = 1.;
496
497			if ( ph->Pt() <= ptmin ) passCuts = -1.;
498
499			// not needed anymore, do in pre-selection...
500			//if ( ph->SCluster()->AbsEta()>=2.5 \|\| (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) passCuts = -1.;
501
502			if( ! ( tIso1 < cic4_allcuts_temp_sublead[_tCat-1+0*4]
503			&& tIso2 < cic4_allcuts_temp_sublead[_tCat-1+1*4]
504			&& tIso3 < cic4_allcuts_temp_sublead[_tCat-1+2*4]
505			&& covIEtaIEta < cic4_allcuts_temp_sublead[_tCat-1+3*4]
506			&& HoE < cic4_allcuts_temp_sublead[_tCat-1+4*4]
507			&& R9 > cic4_allcuts_temp_sublead[_tCat-1+5*4]
508	fabstoec	1.9	&& ( dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4] \|\| !applyEleVeto ) ) ) passCuts = -1.;
509	fabstoec	1.5
510	fabstoec	1.15	if(print) std::cout<<" ---> "<<passCuts<<std::endl;
511
512	fabstoec	1.5	if(kin) {
513	bendavid	1.11	kin[20] = passCuts;
514			kin[21] = (float) _tCat;
515	fabstoec	1.5	}
516
517			if(passCuts > 0.) return true;
518			return false;
519			}
520	bendavid	1.11
521	bendavid	1.13	const MCParticle PhotonTools::MatchMC(const Particle ph, const MCParticleCol *c, Bool_t matchElectrons) {
522
523			// printf("Start loop\n");
524	bendavid	1.11	for (UInt_t i=0; i<c->GetEntries(); ++i) {
525			const MCParticle *p = c->At(i);
526	bendavid	1.13	// if (p->IsGenerated() && p->AbsPdgId()==11 && (p->DistinctMother()->AbsPdgId()==23\|\| p->DistinctMother()->AbsPdgId()==24)) {
527			// printf("pdgid = %i, status = %i, pt = %5f, mother pdgid = %i\n",p->PdgId(),p->Status(),p->Pt(),p->Mother()->PdgId());
528	bendavid	1.12	// }
529			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)) {
530			return p;
531			}
532	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) ) {
533			return p;
534			}
535			}
536			return 0;
537	bendavid	1.12	}
538
539			PhotonTools::DiphotonR9EtaPtCats PhotonTools::DiphotonR9EtaPtCat(const Photon p1, const Photon p2) {
540
541			PhotonTools::CiCBaseLineCats cat1 = CiCBaseLineCat(p1);
542			PhotonTools::CiCBaseLineCats cat2 = CiCBaseLineCat(p2);
543
544			PhotonTools::DiphotonR9EtaPtCats evtcat = PhotonTools::kOctCat0;
545
546			bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 \|\| cat1 == PhotonTools::kCiCCat2);
547			bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 \|\| cat2 == PhotonTools::kCiCCat2);
548
549			bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 \|\| cat1 == PhotonTools::kCiCCat3);
550			bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 \|\| cat2 == PhotonTools::kCiCCat3);
551
552			if( ph1IsEB && ph2IsEB )
553			evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat0 : PhotonTools::kOctCat1);
554			else
555			evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat2 : PhotonTools::kOctCat3);
556
557			float ptgg = (p1->Mom()+p2->Mom()).Pt();
558			if (ptgg<40.0) evtcat = PhotonTools::DiphotonR9EtaPtCats(evtcat + 4);
559
560
561			return evtcat;
562	fabstoec	1.15	}