Mods/src/PhotonIDMod.cc

// $Id: PhotonIDMod.cc,v 1.28 2011/12/17 20:00:40 bendavid Exp $

#include "TDataMember.h"
#include "TTree.h"
#include "MitPhysics/Mods/interface/PhotonIDMod.h"
#include "MitAna/DataTree/interface/PhotonCol.h"
#include "MitPhysics/Init/interface/ModNames.h"
#include "MitPhysics/Utils/interface/IsolationTools.h"
#include "MitPhysics/Utils/interface/PhotonTools.h"

#include <TSystem.h>

using namespace mithep;

ClassImp(mithep::PhotonIDMod)

//--------------------------------------------------------------------------------------------------
PhotonIDMod::PhotonIDMod(const char *name, const char *title) : 
  BaseMod(name,title),
  fPhotonBranchName  (Names::gkPhotonBrn),
  fGoodPhotonsName   (ModNames::gkGoodPhotonsName),
  fTrackBranchName   (Names::gkTrackBrn),
  fBeamspotBranchName(Names::gkBeamSpotBrn),
  fPileUpDenName     (Names::gkPileupEnergyDensityBrn),
  fConversionName    ("MergedConversions"),
  fElectronName      ("Electrons"),
  fGoodElectronName  ("Electrons"),
  fPVName            (Names::gkPVBeamSpotBrn),
  // MC specific stuff...
  fMCParticleName    (Names::gkMCPartBrn),
  fPileUpName        (Names::gkPileupInfoBrn),

  fPhotonIDType      ("Custom"),
  fPhotonIsoType     ("Custom"),
  fPhotonPtMin       (15.0),
  fHadOverEmMax      (0.02),
  fApplySpikeRemoval (kFALSE),
  fApplyPixelSeed    (kTRUE),
  fApplyElectronVeto (kFALSE),
  fInvertElectronVeto(kFALSE),
  fApplyElectronVetoConvRecovery(kFALSE),
  fApplyConversionId (kFALSE),
  fApplyTriggerMatching(kFALSE),
  fPhotonR9Min       (0.5),
  fPhIdType          (kIdUndef),
  fPhIsoType         (kIsoUndef),
  fFiduciality       (kTRUE),
  fEtaWidthEB        (0.01),
  fEtaWidthEE        (0.028),
  fAbsEtaMax         (999.99),
  fApplyR9Min        (kFALSE),

  fEffAreaEcalEE     (0.089),
  fEffAreaHcalEE     (0.156),
  fEffAreaTrackEE    (0.261),

  fEffAreaEcalEB     (0.183),
  fEffAreaHcalEB     (0.062),
  fEffAreaTrackEB    (0.306),

  fPhotons(0),
  fTracks(0),
  fBeamspots(0),
  fPileUpDen(0),
  fConversions(0),
  fElectrons(0),
  fGoodElectrons(0),
  fPV(0),
  fMCParticles       (0),
  fPileUp            (0),

  // MVA ID Stuff
  fbdtCutBarrel      (0.0744), //cuts give the same effiiciency (relative to preselection) with cic
  fbdtCutEndcap      (0.0959), //cuts give the same effiiciency (relative to preselection) with cic  
  fVariableType      (10), //please use 4 which is the correct type
  fEndcapWeights      (gSystem->Getenv("CMSSW_BASE")+TString("/src/MitPhysics/data/TMVAClassificationPhotonID_Endcap_PassPreSel_Variable_10_BDTnCuts2000_BDT.weights.xml")),
  fBarrelWeights      (gSystem->Getenv("CMSSW_BASE")+TString("/src/MitPhysics/data/TMVAClassificationPhotonID_Barrel_PassPreSel_Variable_10_BDTnCuts2000_BDT.weights.xml")),


  fPVFromBranch      (true),
  fGoodElectronsFromBranch (kTRUE),
  fIsData(false)
  
{
  // Constructor.
}

//--------------------------------------------------------------------------------------------------
void PhotonIDMod::Process()
{
  // Process entries of the tree. 

  LoadEventObject(fPhotonBranchName,   fPhotons);
  
  const BaseVertex *bsp = NULL;
  Double_t _tRho = -1.;
  const TriggerObjectCol *trigObjs = 0;  
  if (fPhotons->GetEntries()>0) {
    LoadEventObject(fTrackBranchName,    fTracks);
    LoadEventObject(fBeamspotBranchName, fBeamspots);
    LoadEventObject(fPileUpDenName,      fPileUpDen);
    LoadEventObject(fConversionName,     fConversions);
    LoadEventObject(fElectronName,       fElectrons);
    LoadEventObject(fGoodElectronName,   fGoodElectrons);
    LoadEventObject(fPVName,             fPV);
    
    if(fBeamspots->GetEntries() > 0)
      bsp = fBeamspots->At(0);

    if(fPileUpDen->GetEntries() > 0)
      _tRho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();

    
    //get trigger object collection if trigger matching is enabled
    if (fApplyTriggerMatching) {
      trigObjs = GetHLTObjects(fTrigObjectsName);
    }    
  
  }
  
  PhotonOArr *GoodPhotons = new PhotonOArr;
  GoodPhotons->SetName(fGoodPhotonsName);

  for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {    
    const Photon *ph = fPhotons->At(i);        

    
    if (fFiduciality == kTRUE &&
        (ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566) ) ) 
      continue;
    
    if (fInvertElectronVeto && PhotonTools::PassElectronVeto(ph,fGoodElectrons)) {
      continue;
    }
    
    // ---------------------------------------------------------------------
    // check if we use the CiC Selection. If yes, bypass all the below...
    if(fPhIdType == kBaseLineCiC) {
      if( PhotonTools::PassCiCSelection(ph, fPV->At(0), fTracks, fElectrons, fPV, _tRho, fPhotonPtMin, fApplyElectronVeto) )
        GoodPhotons->Add(fPhotons->At(i));
      continue; // go to next Photons
    }
    // ---------------------------------------------------------------------

    //loose photon preselection for subsequent mva
    if(fPhIdType == kMITPhSelection ) {
      if( ph->Pt()>fPhotonPtMin && PhotonTools::PassSinglePhotonPresel(ph,fElectrons,fConversions,bsp,fTracks,fPV->At(0),_tRho,fApplyElectronVeto) ) {
        GoodPhotons->Add(fPhotons->At(i));
      }
      continue;
    }

    // add MingMings MVA ID on single Photon level
    if(fPhIdType == kMITMVAId ) {
      if( ph->Pt()>fPhotonPtMin && PhotonTools::PassSinglePhotonPresel(ph,fElectrons,fConversions,bsp,fTracks,fPV->At(0),_tRho,fApplyElectronVeto) && fTool.PassMVASelection(ph, fPV->At(0) ,fTracks, fPV, _tRho ,fbdtCutBarrel,fbdtCutEndcap, fElectrons, fApplyElectronVeto) ) {
        GoodPhotons->Add(fPhotons->At(i));
      }
      continue;
    }  // go to next Photon

    if (ph->Pt() <= fPhotonPtMin) 
      continue;    // add good electron

    Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
    
    Bool_t passSpikeRemovalFilter = kTRUE;

    if (ph->SCluster() && ph->SCluster()->Seed()) {
      if(ph->SCluster()->Seed()->Energy() > 5.0 && 
         ph->SCluster()->Seed()->EMax() / ph->SCluster()->Seed()->E3x3() > 0.95
        ) {
        passSpikeRemovalFilter = kFALSE;
      }
    }

    // For Now Only use the EMax/E3x3 prescription.
    //if(ph->SCluster()->Seed()->Energy() > 5.0 && 
    //   (1 - (ph->SCluster()->Seed()->E1x3() + ph->SCluster()->Seed()->E3x1() - 2*ph->SCluster()->Seed()->EMax())) > 0.95
    //  ) {
    //  passSpikeRemovalFilter = kFALSE;
    //}
    if (fApplySpikeRemoval && !passSpikeRemovalFilter) continue;
    
    if (ph->HadOverEm() >= fHadOverEmMax) 
      continue;

    if (fApplyPixelSeed == kTRUE &&
        ph->HasPixelSeed() == kTRUE) 
      continue;

    if (fApplyElectronVeto && !PhotonTools::PassElectronVeto(ph,fElectrons) ) continue;

    if (fApplyElectronVetoConvRecovery && !PhotonTools::PassElectronVetoConvRecovery(ph,fElectrons,fConversions,bsp) ) continue;

    if (fApplyConversionId && !PhotonTools::PassConversionId(ph,PhotonTools::MatchedConversion(ph,fConversions,bsp))) continue;

    if (fApplyTriggerMatching && !PhotonTools::PassTriggerMatching(ph,trigObjs)) continue;

    Bool_t idcut = kFALSE;
    switch (fPhIdType) {
      case kTight:
        idcut = ph->IsTightPhoton();
        break;
      case kLoose:
        idcut = ph->IsLoosePhoton();
       break;
      case kLooseEM:
        idcut = ph->IsLooseEM();
        break;
      case kCustomId:
        idcut = kTRUE;
      default:
        break;
    }

    if (!idcut) 
      continue;

    Bool_t isocut = kFALSE;
    switch (fPhIsoType) {      
      case kNoIso:
        isocut = kTRUE;
        break;
      case kCombinedIso:
        {
          Double_t totalIso = ph->HollowConeTrkIsoDr04()+
                              ph->EcalRecHitIsoDr04() +
                              ph->HcalTowerSumEtDr04();
          if (totalIso/ph->Pt() < 0.25)
            isocut = kTRUE;
        }
        break;
      case kCustomIso:
        {
          if ( ph->HollowConeTrkIsoDr04() < (1.5 + 0.001*ph->Pt()) && ph->EcalRecHitIsoDr04()<(2.0+0.006*ph->Pt()) && ph->HcalTowerSumEtDr04()<(2.0+0.0025*ph->Pt()) )
            isocut = kTRUE;
        }
        break;
        
    case kMITPUCorrected:
      {
        // compute the PU corrections only if Rho is available
        // ... otherwise (_tRho = -1.0) it's the std isolation
        isocut = kTRUE;
        Double_t fEffAreaEcal = fEffAreaEcalEB;
        Double_t fEffAreaHcal = fEffAreaHcalEB;
        Double_t fEffAreaTrack = fEffAreaTrackEB;

        if( !isbarrel ) {
          fEffAreaEcal = fEffAreaEcalEE;
          fEffAreaHcal = fEffAreaHcalEE;
          fEffAreaTrack = fEffAreaTrackEE;
        }         

        //std::cout<<"-----------------------------------"<<std::endl;
        //std::cout<<" MIT phooton Et = "<<ph->Et()<<std::endl;           
        Double_t EcalCorrISO =   ph->EcalRecHitIsoDr04();
        //std::cout<<" ecaliso4    = "<<EcalCorrISO<<std::endl;
        if(_tRho > -0.5 ) EcalCorrISO -= _tRho * fEffAreaEcal;
        //std::cout<<" ecaliso4Corr = "<<EcalCorrISO<<"  ("<<2.0+0.006*ph->Pt()<<")"<<std::endl;
        if ( EcalCorrISO > (2.0+0.006*ph->Pt()) ) isocut = kFALSE; 
        if ( isocut || true ) {
          Double_t HcalCorrISO = ph->HcalTowerSumEtDr04(); 
          //std::cout<<" hcaliso4     = "<<HcalCorrISO<<std::endl;
          if(_tRho > -0.5 ) HcalCorrISO -= _tRho * fEffAreaHcal;
          //std::cout<<" hcaliso4Corr = "<<HcalCorrISO<<"  ("<<2.0+0.0025*ph->Pt()<<")"<<std::endl;
          if ( HcalCorrISO > (2.0+0.0025*ph->Pt()) ) isocut = kFALSE;
        }
        if ( isocut || true ) {
          Double_t TrackCorrISO = IsolationTools::TrackIsolationNoPV(ph, bsp, 0.4, 0.04, 0.0, 0.015, 0.1, TrackQuality::highPurity, fTracks);
          //std::cout<<" TrackIso       = "<<TrackCorrISO<<std::endl;     
          if(_tRho > -0.5 )
            TrackCorrISO -= _tRho * fEffAreaTrack;
          //std::cout<<" TrackIsoCorr   = "<<TrackCorrISO<<"  ("<<1.5 + 0.001*ph->Pt()<<")"<<std::endl;   
          if ( TrackCorrISO > (1.5 + 0.001*ph->Pt()) ) isocut = kFALSE;
          //std::cout<<"      passes ? "<<isocut<<std::endl;
        }
        break;
      }
    default:
      break;
    }
    
    if (!isocut) 
      continue;
    
    if ( fApplyR9Min && ph->R9() <= fPhotonR9Min) 
      continue;

    if ((isbarrel && ph->CoviEtaiEta() >= fEtaWidthEB) ||
        (!isbarrel && ph->CoviEtaiEta() >= fEtaWidthEE))
      continue;

    if (ph->AbsEta() >= fAbsEtaMax) 
      continue;
      
    
    // add good electron
    GoodPhotons->Add(fPhotons->At(i));
  }

  // sort according to pt
  GoodPhotons->Sort();

  // add to event for other modules to use
  AddObjThisEvt(GoodPhotons);  
  
  return; 
  
}

//--------------------------------------------------------------------------------------------------
void PhotonIDMod::SlaveBegin()
{
  // Run startup code on the computer (slave) doing the actual analysis. Here,
  // we just request the photon collection branch.

  ReqEventObject(fPhotonBranchName,   fPhotons,     kTRUE);
  ReqEventObject(fTrackBranchName,    fTracks,      kTRUE);
  ReqEventObject(fBeamspotBranchName, fBeamspots,   kTRUE);
  ReqEventObject(fConversionName,     fConversions, kTRUE);
  ReqEventObject(fElectronName,       fElectrons,   kTRUE);
  ReqEventObject(fGoodElectronName,       fGoodElectrons,   fGoodElectronsFromBranch);  
  ReqEventObject(fPVName, fPV, fPVFromBranch);
  ReqEventObject(fPileUpDenName,      fPileUpDen, kTRUE);
  if (!fIsData) {
    ReqBranch(fPileUpName,            fPileUp);
    ReqBranch(fMCParticleName,        fMCParticles);
  }   
  
  if (fPhotonIDType.CompareTo("Tight") == 0) 
    fPhIdType = kTight;
  else if (fPhotonIDType.CompareTo("Loose") == 0) 
    fPhIdType = kLoose;
  else if (fPhotonIDType.CompareTo("LooseEM") == 0) 
    fPhIdType = kLooseEM;
  else if (fPhotonIDType.CompareTo("Custom") == 0)
    fPhIdType = kCustomId;
  else if (fPhotonIDType.CompareTo("BaseLineCiC") == 0) {
    fPhIdType = kBaseLineCiC;
    fPhotonIsoType = "NoIso";
  }
  else if (fPhotonIDType.CompareTo("MITMVAId") == 0) {
    fPhIdType = kMITMVAId;
    fPhotonIsoType = "NoIso";
    fTool.InitializeMVA(fVariableType,fEndcapWeights,fBarrelWeights);
  }
  else if (fPhotonIDType.CompareTo("MITSelection") == 0)  {
    fPhIdType = kMITPhSelection;
    fPhotonIsoType = "NoIso";    
  }
  else {
    SendError(kAbortAnalysis, "SlaveBegin",
              "The specified photon identification %s is not defined.",
              fPhotonIDType.Data());
    return;
  }

  if (fPhotonIsoType.CompareTo("NoIso") == 0 )
    fPhIsoType = kNoIso;
  else if (fPhotonIsoType.CompareTo("CombinedIso") == 0 )
    fPhIsoType = kCombinedIso;
  else if (fPhotonIsoType.CompareTo("Custom") == 0 )
    fPhIsoType = kCustomIso;
  else if (fPhotonIsoType.CompareTo("MITPUCorrected") == 0 )
    fPhIsoType = kMITPUCorrected;
  else {
    SendError(kAbortAnalysis, "SlaveBegin",
              "The specified photon isolation %s is not defined.",
              fPhotonIsoType.Data());
    return;
  }
    
  
}
Revision:	1.29
Committed:	Mon Dec 19 23:45:00 2011 UTC (13 years, 4 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
Changes since 1.28:	+4 -4 lines
Log Message:	updated id mva
#	User	Rev	Content
1	bendavid	1.29	// $Id: PhotonIDMod.cc,v 1.28 2011/12/17 20:00:40 bendavid Exp $
2	sixie	1.1
3	bendavid	1.24	#include "TDataMember.h"
4			#include "TTree.h"
5	sixie	1.1	#include "MitPhysics/Mods/interface/PhotonIDMod.h"
6	loizides	1.7	#include "MitAna/DataTree/interface/PhotonCol.h"
7	sixie	1.1	#include "MitPhysics/Init/interface/ModNames.h"
8	fabstoec	1.17	#include "MitPhysics/Utils/interface/IsolationTools.h"
9	bendavid	1.19	#include "MitPhysics/Utils/interface/PhotonTools.h"
10	sixie	1.1
11	fabstoec	1.26	#include <TSystem.h>
12
13	sixie	1.1	using namespace mithep;
14
15			ClassImp(mithep::PhotonIDMod)
16
17			//--------------------------------------------------------------------------------------------------
18			PhotonIDMod::PhotonIDMod(const char name, const char title) :
19			BaseMod(name,title),
20	fabstoec	1.17	fPhotonBranchName (Names::gkPhotonBrn),
21			fGoodPhotonsName (ModNames::gkGoodPhotonsName),
22			fTrackBranchName (Names::gkTrackBrn),
23			fBeamspotBranchName(Names::gkBeamSpotBrn),
24	fabstoec	1.18	fPileUpDenName (Names::gkPileupEnergyDensityBrn),
25	bendavid	1.19	fConversionName ("MergedConversions"),
26			fElectronName ("Electrons"),
27	bendavid	1.25	fGoodElectronName ("Electrons"),
28	fabstoec	1.22	fPVName (Names::gkPVBeamSpotBrn),
29	bendavid	1.24	// MC specific stuff...
30			fMCParticleName (Names::gkMCPartBrn),
31			fPileUpName (Names::gkPileupInfoBrn),
32	fabstoec	1.22
33			fPhotonIDType ("Custom"),
34			fPhotonIsoType ("Custom"),
35			fPhotonPtMin (15.0),
36			fHadOverEmMax (0.02),
37			fApplySpikeRemoval (kFALSE),
38			fApplyPixelSeed (kTRUE),
39			fApplyElectronVeto (kFALSE),
40	bendavid	1.25	fInvertElectronVeto(kFALSE),
41	bendavid	1.19	fApplyElectronVetoConvRecovery(kFALSE),
42	fabstoec	1.22	fApplyConversionId (kFALSE),
43	bendavid	1.19	fApplyTriggerMatching(kFALSE),
44	fabstoec	1.22	fPhotonR9Min (0.5),
45			fPhIdType (kIdUndef),
46			fPhIsoType (kIsoUndef),
47			fFiduciality (kTRUE),
48			fEtaWidthEB (0.01),
49			fEtaWidthEE (0.028),
50			fAbsEtaMax (999.99),
51			fApplyR9Min (kFALSE),
52
53			fEffAreaEcalEE (0.089),
54			fEffAreaHcalEE (0.156),
55			fEffAreaTrackEE (0.261),
56
57			fEffAreaEcalEB (0.183),
58			fEffAreaHcalEB (0.062),
59			fEffAreaTrackEB (0.306),
60
61	fabstoec	1.17	fPhotons(0),
62			fTracks(0),
63			fBeamspots(0),
64	bendavid	1.19	fPileUpDen(0),
65			fConversions(0),
66	fabstoec	1.22	fElectrons(0),
67	bendavid	1.25	fGoodElectrons(0),
68	fabstoec	1.22	fPV(0),
69	bendavid	1.24	fMCParticles (0),
70			fPileUp (0),
71	fabstoec	1.22
72	fabstoec	1.26	// MVA ID Stuff
73	bendavid	1.27	fbdtCutBarrel (0.0744), //cuts give the same effiiciency (relative to preselection) with cic
74			fbdtCutEndcap (0.0959), //cuts give the same effiiciency (relative to preselection) with cic
75	bendavid	1.29	fVariableType (10), //please use 4 which is the correct type
76			fEndcapWeights (gSystem->Getenv("CMSSW_BASE")+TString("/src/MitPhysics/data/TMVAClassificationPhotonID_Endcap_PassPreSel_Variable_10_BDTnCuts2000_BDT.weights.xml")),
77			fBarrelWeights (gSystem->Getenv("CMSSW_BASE")+TString("/src/MitPhysics/data/TMVAClassificationPhotonID_Barrel_PassPreSel_Variable_10_BDTnCuts2000_BDT.weights.xml")),
78	bendavid	1.27
79
80
81	fabstoec	1.26
82
83	bendavid	1.24	fPVFromBranch (true),
84	bendavid	1.25	fGoodElectronsFromBranch (kTRUE),
85			fIsData(false)
86	fabstoec	1.17
87	sixie	1.1	{
88			// Constructor.
89			}
90
91			//--------------------------------------------------------------------------------------------------
92			void PhotonIDMod::Process()
93			{
94			// Process entries of the tree.
95
96	fabstoec	1.17	LoadEventObject(fPhotonBranchName, fPhotons);
97	bendavid	1.19
98			const BaseVertex *bsp = NULL;
99			Double_t _tRho = -1.;
100			const TriggerObjectCol *trigObjs = 0;
101			if (fPhotons->GetEntries()>0) {
102			LoadEventObject(fTrackBranchName, fTracks);
103			LoadEventObject(fBeamspotBranchName, fBeamspots);
104	bendavid	1.24	LoadEventObject(fPileUpDenName, fPileUpDen);
105	bendavid	1.19	LoadEventObject(fConversionName, fConversions);
106			LoadEventObject(fElectronName, fElectrons);
107	fabstoec	1.26	LoadEventObject(fGoodElectronName, fGoodElectrons);
108	bendavid	1.24	LoadEventObject(fPVName, fPV);
109	bendavid	1.19
110			if(fBeamspots->GetEntries() > 0)
111			bsp = fBeamspots->At(0);
112	fabstoec	1.17
113	bendavid	1.24	if(fPileUpDen->GetEntries() > 0)
114	fabstoec	1.23	_tRho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
115	fabstoec	1.26
116	bendavid	1.19
117			//get trigger object collection if trigger matching is enabled
118			if (fApplyTriggerMatching) {
119			trigObjs = GetHLTObjects(fTrigObjectsName);
120			}
121
122			}
123	fabstoec	1.17
124	sixie	1.1	PhotonOArr *GoodPhotons = new PhotonOArr;
125			GoodPhotons->SetName(fGoodPhotonsName);
126
127			for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
128			const Photon *ph = fPhotons->At(i);
129
130	bendavid	1.24
131			if (fFiduciality == kTRUE &&
132			(ph->SCluster()->AbsEta()>=2.5 \|\| (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566) ) )
133			continue;
134
135	bendavid	1.25	if (fInvertElectronVeto && PhotonTools::PassElectronVeto(ph,fGoodElectrons)) {
136			continue;
137			}
138
139	fabstoec	1.22	// ---------------------------------------------------------------------
140			// check if we use the CiC Selection. If yes, bypass all the below...
141			if(fPhIdType == kBaseLineCiC) {
142			if( PhotonTools::PassCiCSelection(ph, fPV->At(0), fTracks, fElectrons, fPV, _tRho, fPhotonPtMin, fApplyElectronVeto) )
143			GoodPhotons->Add(fPhotons->At(i));
144			continue; // go to next Photons
145			}
146			// ---------------------------------------------------------------------
147
148	bendavid	1.27	//loose photon preselection for subsequent mva
149			if(fPhIdType == kMITPhSelection ) {
150	bendavid	1.28	if( ph->Pt()>fPhotonPtMin && PhotonTools::PassSinglePhotonPresel(ph,fElectrons,fConversions,bsp,fTracks,fPV->At(0),_tRho,fApplyElectronVeto) ) {
151	bendavid	1.27	GoodPhotons->Add(fPhotons->At(i));
152			}
153			continue;
154			}
155
156	fabstoec	1.26	// add MingMings MVA ID on single Photon level
157			if(fPhIdType == kMITMVAId ) {
158	bendavid	1.28	if( ph->Pt()>fPhotonPtMin && PhotonTools::PassSinglePhotonPresel(ph,fElectrons,fConversions,bsp,fTracks,fPV->At(0),_tRho,fApplyElectronVeto) && fTool.PassMVASelection(ph, fPV->At(0) ,fTracks, fPV, _tRho ,fbdtCutBarrel,fbdtCutEndcap, fElectrons, fApplyElectronVeto) ) {
159	fabstoec	1.26	GoodPhotons->Add(fPhotons->At(i));
160			}
161			continue;
162			} // go to next Photon
163
164	sixie	1.1	if (ph->Pt() <= fPhotonPtMin)
165	fabstoec	1.22	continue; // add good electron
166	sixie	1.13
167	bendavid	1.21	Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
168	sixie	1.1
169	sixie	1.13	Bool_t passSpikeRemovalFilter = kTRUE;
170	sixie	1.15
171			if (ph->SCluster() && ph->SCluster()->Seed()) {
172			if(ph->SCluster()->Seed()->Energy() > 5.0 &&
173			ph->SCluster()->Seed()->EMax() / ph->SCluster()->Seed()->E3x3() > 0.95
174			) {
175			passSpikeRemovalFilter = kFALSE;
176			}
177	sixie	1.13	}
178
179			// For Now Only use the EMax/E3x3 prescription.
180			//if(ph->SCluster()->Seed()->Energy() > 5.0 &&
181			// (1 - (ph->SCluster()->Seed()->E1x3() + ph->SCluster()->Seed()->E3x1() - 2*ph->SCluster()->Seed()->EMax())) > 0.95
182			// ) {
183			// passSpikeRemovalFilter = kFALSE;
184			//}
185			if (fApplySpikeRemoval && !passSpikeRemovalFilter) continue;
186
187	ceballos	1.2	if (ph->HadOverEm() >= fHadOverEmMax)
188			continue;
189	ceballos	1.14
190	ceballos	1.2	if (fApplyPixelSeed == kTRUE &&
191			ph->HasPixelSeed() == kTRUE)
192			continue;
193	ceballos	1.14
194	bendavid	1.19	if (fApplyElectronVeto && !PhotonTools::PassElectronVeto(ph,fElectrons) ) continue;
195
196			if (fApplyElectronVetoConvRecovery && !PhotonTools::PassElectronVetoConvRecovery(ph,fElectrons,fConversions,bsp) ) continue;
197
198			if (fApplyConversionId && !PhotonTools::PassConversionId(ph,PhotonTools::MatchedConversion(ph,fConversions,bsp))) continue;
199
200			if (fApplyTriggerMatching && !PhotonTools::PassTriggerMatching(ph,trigObjs)) continue;
201
202	sixie	1.1	Bool_t idcut = kFALSE;
203			switch (fPhIdType) {
204			case kTight:
205			idcut = ph->IsTightPhoton();
206			break;
207			case kLoose:
208			idcut = ph->IsLoosePhoton();
209			break;
210			case kLooseEM:
211			idcut = ph->IsLooseEM();
212			break;
213			case kCustomId:
214	ceballos	1.5	idcut = kTRUE;
215	sixie	1.1	default:
216			break;
217			}
218
219			if (!idcut)
220			continue;
221
222			Bool_t isocut = kFALSE;
223			switch (fPhIsoType) {
224			case kNoIso:
225			isocut = kTRUE;
226			break;
227	ceballos	1.2	case kCombinedIso:
228	bendavid	1.12	{
229			Double_t totalIso = ph->HollowConeTrkIsoDr04()+
230			ph->EcalRecHitIsoDr04() +
231			ph->HcalTowerSumEtDr04();
232			if (totalIso/ph->Pt() < 0.25)
233			isocut = kTRUE;
234			}
235	ceballos	1.2	break;
236	sixie	1.1	case kCustomIso:
237	bendavid	1.16	{
238			if ( ph->HollowConeTrkIsoDr04() < (1.5 + 0.001ph->Pt()) && ph->EcalRecHitIsoDr04()<(2.0+0.006ph->Pt()) && ph->HcalTowerSumEtDr04()<(2.0+0.0025*ph->Pt()) )
239			isocut = kTRUE;
240			}
241	fabstoec	1.17	break;
242
243			case kMITPUCorrected:
244			{
245			// compute the PU corrections only if Rho is available
246			// ... otherwise (_tRho = -1.0) it's the std isolation
247			isocut = kTRUE;
248	fabstoec	1.20	Double_t fEffAreaEcal = fEffAreaEcalEB;
249			Double_t fEffAreaHcal = fEffAreaHcalEB;
250			Double_t fEffAreaTrack = fEffAreaTrackEB;
251
252	bendavid	1.21	if( !isbarrel ) {
253	fabstoec	1.20	fEffAreaEcal = fEffAreaEcalEE;
254			fEffAreaHcal = fEffAreaHcalEE;
255			fEffAreaTrack = fEffAreaTrackEE;
256			}
257	fabstoec	1.26
258			//std::cout<<"-----------------------------------"<<std::endl;
259			//std::cout<<" MIT phooton Et = "<<ph->Et()<<std::endl;
260	fabstoec	1.17	Double_t EcalCorrISO = ph->EcalRecHitIsoDr04();
261	fabstoec	1.26	//std::cout<<" ecaliso4 = "<<EcalCorrISO<<std::endl;
262	fabstoec	1.17	if(_tRho > -0.5 ) EcalCorrISO -= _tRho * fEffAreaEcal;
263	fabstoec	1.26	//std::cout<<" ecaliso4Corr = "<<EcalCorrISO<<" ("<<2.0+0.006*ph->Pt()<<")"<<std::endl;
264	fabstoec	1.17	if ( EcalCorrISO > (2.0+0.006*ph->Pt()) ) isocut = kFALSE;
265	fabstoec	1.26	if ( isocut \|\| true ) {
266	fabstoec	1.17	Double_t HcalCorrISO = ph->HcalTowerSumEtDr04();
267	fabstoec	1.26	//std::cout<<" hcaliso4 = "<<HcalCorrISO<<std::endl;
268	fabstoec	1.17	if(_tRho > -0.5 ) HcalCorrISO -= _tRho * fEffAreaHcal;
269	fabstoec	1.26	//std::cout<<" hcaliso4Corr = "<<HcalCorrISO<<" ("<<2.0+0.0025*ph->Pt()<<")"<<std::endl;
270	fabstoec	1.17	if ( HcalCorrISO > (2.0+0.0025*ph->Pt()) ) isocut = kFALSE;
271			}
272	fabstoec	1.26	if ( isocut \|\| true ) {
273	fabstoec	1.17	Double_t TrackCorrISO = IsolationTools::TrackIsolationNoPV(ph, bsp, 0.4, 0.04, 0.0, 0.015, 0.1, TrackQuality::highPurity, fTracks);
274	fabstoec	1.26	//std::cout<<" TrackIso = "<<TrackCorrISO<<std::endl;
275	fabstoec	1.17	if(_tRho > -0.5 )
276			TrackCorrISO -= _tRho * fEffAreaTrack;
277	fabstoec	1.26	//std::cout<<" TrackIsoCorr = "<<TrackCorrISO<<" ("<<1.5 + 0.001*ph->Pt()<<")"<<std::endl;
278	fabstoec	1.17	if ( TrackCorrISO > (1.5 + 0.001*ph->Pt()) ) isocut = kFALSE;
279	fabstoec	1.26	//std::cout<<" passes ? "<<isocut<<std::endl;
280	fabstoec	1.17	}
281			break;
282			}
283			default:
284			break;
285	sixie	1.1	}
286	fabstoec	1.17
287	sixie	1.1	if (!isocut)
288			continue;
289	fabstoec	1.17
290	bendavid	1.16	if ( fApplyR9Min && ph->R9() <= fPhotonR9Min)
291	ceballos	1.5	continue;
292
293	bendavid	1.21	if ((isbarrel && ph->CoviEtaiEta() >= fEtaWidthEB) \|\|
294			(!isbarrel && ph->CoviEtaiEta() >= fEtaWidthEE))
295	ceballos	1.10	continue;
296
297	ceballos	1.11	if (ph->AbsEta() >= fAbsEtaMax)
298			continue;
299	bendavid	1.25
300
301	sixie	1.1	// add good electron
302			GoodPhotons->Add(fPhotons->At(i));
303			}
304
305	loizides	1.4	// sort according to pt
306			GoodPhotons->Sort();
307
308	sixie	1.1	// add to event for other modules to use
309			AddObjThisEvt(GoodPhotons);
310	bendavid	1.24
311	bendavid	1.25	return;
312	bendavid	1.24
313	sixie	1.1	}
314
315			//--------------------------------------------------------------------------------------------------
316			void PhotonIDMod::SlaveBegin()
317			{
318			// Run startup code on the computer (slave) doing the actual analysis. Here,
319	loizides	1.3	// we just request the photon collection branch.
320	sixie	1.1
321	fabstoec	1.22	ReqEventObject(fPhotonBranchName, fPhotons, kTRUE);
322			ReqEventObject(fTrackBranchName, fTracks, kTRUE);
323			ReqEventObject(fBeamspotBranchName, fBeamspots, kTRUE);
324	bendavid	1.19	ReqEventObject(fConversionName, fConversions, kTRUE);
325	fabstoec	1.22	ReqEventObject(fElectronName, fElectrons, kTRUE);
326	bendavid	1.25	ReqEventObject(fGoodElectronName, fGoodElectrons, fGoodElectronsFromBranch);
327	bendavid	1.24	ReqEventObject(fPVName, fPV, fPVFromBranch);
328			ReqEventObject(fPileUpDenName, fPileUpDen, kTRUE);
329			if (!fIsData) {
330			ReqBranch(fPileUpName, fPileUp);
331			ReqBranch(fMCParticleName, fMCParticles);
332			}
333
334	sixie	1.1	if (fPhotonIDType.CompareTo("Tight") == 0)
335			fPhIdType = kTight;
336			else if (fPhotonIDType.CompareTo("Loose") == 0)
337			fPhIdType = kLoose;
338			else if (fPhotonIDType.CompareTo("LooseEM") == 0)
339			fPhIdType = kLooseEM;
340	fabstoec	1.22	else if (fPhotonIDType.CompareTo("Custom") == 0)
341	sixie	1.1	fPhIdType = kCustomId;
342	fabstoec	1.22	else if (fPhotonIDType.CompareTo("BaseLineCiC") == 0) {
343			fPhIdType = kBaseLineCiC;
344			fPhotonIsoType = "NoIso";
345	fabstoec	1.26	}
346			else if (fPhotonIDType.CompareTo("MITMVAId") == 0) {
347			fPhIdType = kMITMVAId;
348			fPhotonIsoType = "NoIso";
349			fTool.InitializeMVA(fVariableType,fEndcapWeights,fBarrelWeights);
350	bendavid	1.27	}
351			else if (fPhotonIDType.CompareTo("MITSelection") == 0) {
352			fPhIdType = kMITPhSelection;
353			fPhotonIsoType = "NoIso";
354			}
355			else {
356	sixie	1.1	SendError(kAbortAnalysis, "SlaveBegin",
357	ceballos	1.2	"The specified photon identification %s is not defined.",
358	sixie	1.1	fPhotonIDType.Data());
359			return;
360			}
361
362			if (fPhotonIsoType.CompareTo("NoIso") == 0 )
363			fPhIsoType = kNoIso;
364	ceballos	1.2	else if (fPhotonIsoType.CompareTo("CombinedIso") == 0 )
365			fPhIsoType = kCombinedIso;
366	fabstoec	1.17	else if (fPhotonIsoType.CompareTo("Custom") == 0 )
367	sixie	1.1	fPhIsoType = kCustomIso;
368	fabstoec	1.17	else if (fPhotonIsoType.CompareTo("MITPUCorrected") == 0 )
369			fPhIsoType = kMITPUCorrected;
370			else {
371	sixie	1.1	SendError(kAbortAnalysis, "SlaveBegin",
372	ceballos	1.2	"The specified photon isolation %s is not defined.",
373	sixie	1.1	fPhotonIsoType.Data());
374			return;
375			}
376	bendavid	1.25
377	bendavid	1.24
378	sixie	1.1	}