Mods/src/PhotonIDMod.cc


// $Id: PhotonIDMod.cc,v 1.29 2011/12/19 23:45:00 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")),


  fDoMCR9Scaling     (kFALSE),
  fMCR9ScaleEB       (1.0),
  fMCR9ScaleEE       (1.0),
  fDoMCSigIEtaIEtaScaling(kFALSE),
  fDoMCWidthScaling(kFALSE),
  fDoMCErrScaling     (kFALSE),
  fMCErrScaleEB       (1.0),
  fMCErrScaleEE       (1.0),

  fPhotonsFromBranch(true),
  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);
  GoodPhotons->SetOwner(kTRUE);

  for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {    
    // need to cpoy the photon in order to be able to scale R9 etc.
    Photon *ph = new Photon(*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) && false) {
      continue;
    }
    
    // -----------------------------------------------------------------------------------
    // Do all the scaling ....
    if (fDoMCErrScaling && !fIsData) {
      if (ph->SCluster()->AbsEta()<1.5) PhotonTools::ScalePhotonError(ph,fMCErrScaleEB);
      else PhotonTools::ScalePhotonError(ph,fMCErrScaleEE);
    }

    if (fDoMCR9Scaling && !fIsData) {
      if (ph->SCluster()->AbsEta()<1.5) PhotonTools::ScalePhotonR9(ph,fMCR9ScaleEB);
      else PhotonTools::ScalePhotonR9(ph,fMCR9ScaleEE);
    }

    if (fDoMCSigIEtaIEtaScaling && !fIsData) {
      if (ph->SCluster()->AbsEta()<1.5) ph->SetCoviEtaiEta(0.87*ph->CoviEtaiEta() + 0.0011);
      else ph->SetCoviEtaiEta(0.99*ph->CoviEtaiEta());
    }

    if (fDoMCWidthScaling && !fIsData) {
      ph->SetEtaWidth(0.99*ph->EtaWidth());
      ph->SetPhiWidth(0.99*ph->PhiWidth());
    }

    
    // ---------------------------------------------------------------------
    // 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->AddOwned(ph);
      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,fInvertElectronVeto) ) {
        GoodPhotons->AddOwned(ph);
      }
      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->AddOwned(ph);
      }
      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->AddOwned(ph);
  }

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