Mods/src/PhotonIDMod.cc

// $Id: PhotonIDMod.cc,v 1.16 2011/02/01 17:02:22 bendavid Exp $

#include "MitPhysics/Mods/interface/PhotonIDMod.h"
#include "MitAna/DataTree/interface/PhotonCol.h"
#include "MitPhysics/Init/interface/ModNames.h"
#include "MitPhysics/Utils/interface/IsolationTools.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     ("Rho"),
  fPhotonIDType("Custom"),
  fPhotonIsoType("Custom"),
  fPhotonPtMin(15.0),
  fHadOverEmMax(0.02),
  fApplySpikeRemoval(kFALSE),
  fApplyPixelSeed(kTRUE),
  fPhotonR9Min(0.5),
  fPhIdType(kIdUndef),
  fPhIsoType(kIsoUndef),
  fFiduciality(kTRUE),
  fEtaWidthEB(0.01),
  fEtaWidthEE(0.028),
  fAbsEtaMax(2.5),
  fApplyR9Min(kFALSE),
  fEffAreaEcal(0.139),
  fEffAreaHcal(0.056),
  fEffAreaTrack(0.296),
  fPhotons(0),
  fTracks(0),
  fBeamspots(0),
  fPileUpDen(0)
  
{
  // Constructor.
}

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

  LoadEventObject(fPhotonBranchName,   fPhotons);
  LoadEventObject(fTrackBranchName,    fTracks);
  LoadEventObject(fBeamspotBranchName, fBeamspots);
  LoadEventObject(fPileUpDenName,      fPileUpDen);

  const BaseVertex *bsp = NULL;
  if(fBeamspots->GetEntries() > 0)
    bsp = dynamic_cast<const BaseVertex*>(fBeamspots->At(0));

  Double_t _tRho = -1.;
  if(fPileUpDen->GetEntries() > 0)
    _tRho = (Double_t) fPileUpDen->At(0)->Rho();
  
  PhotonOArr *GoodPhotons = new PhotonOArr;
  GoodPhotons->SetName(fGoodPhotonsName);

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

    if (ph->Pt() <= fPhotonPtMin) 
      continue;


    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;

    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 EcalCorrISO =   ph->EcalRecHitIsoDr04();
        if(_tRho > -0.5 ) EcalCorrISO -= _tRho * fEffAreaEcal;
        if ( EcalCorrISO > (2.0+0.006*ph->Pt()) ) isocut = kFALSE; 
        if ( isocut ) {
          Double_t HcalCorrISO = ph->HcalTowerSumEtDr04(); 
          if(_tRho > -0.5 ) HcalCorrISO -= _tRho * fEffAreaHcal;
          if ( HcalCorrISO > (2.0+0.0025*ph->Pt()) ) isocut = kFALSE;
        }
        if ( isocut ) {
          Double_t TrackCorrISO = IsolationTools::TrackIsolationNoPV(ph, bsp, 0.4, 0.04, 0.0, 0.015, 0.1, TrackQuality::highPurity, fTracks);
          if(_tRho > -0.5 )
            TrackCorrISO -= _tRho * fEffAreaTrack;
          if ( TrackCorrISO > (1.5 + 0.001*ph->Pt()) ) isocut = kFALSE;
        }
        break;
      }
    default:
      break;
    }
    
    if (!isocut) 
      continue;
    
    if ( fApplyR9Min && ph->R9() <= fPhotonR9Min) 
      continue;

    if (fFiduciality == kTRUE &&
        ph->IsEB() == kFALSE && ph->IsEE() == kFALSE) 
      continue;

    if ((ph->IsEB() == kTRUE && ph->CoviEtaiEta() >= fEtaWidthEB) ||
        (ph->IsEE() == kTRUE && 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);  
}

//--------------------------------------------------------------------------------------------------
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(fPileUpDenName,      fPileUpDen, kTRUE);

  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 {
    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.17
Committed:	Wed Apr 6 18:03:48 2011 UTC (14 years, 1 month ago) by fabstoec
Content type:	text/plain
Branch:	MAIN
Changes since 1.16:	+63 -13 lines
Log Message:	* empty log message *
#	User	Rev	Content
1	fabstoec	1.17	// $Id: PhotonIDMod.cc,v 1.16 2011/02/01 17:02:22 bendavid Exp $
2	sixie	1.1
3			#include "MitPhysics/Mods/interface/PhotonIDMod.h"
4	loizides	1.7	#include "MitAna/DataTree/interface/PhotonCol.h"
5	sixie	1.1	#include "MitPhysics/Init/interface/ModNames.h"
6	fabstoec	1.17	#include "MitPhysics/Utils/interface/IsolationTools.h"
7	sixie	1.1
8			using namespace mithep;
9
10			ClassImp(mithep::PhotonIDMod)
11
12			//--------------------------------------------------------------------------------------------------
13			PhotonIDMod::PhotonIDMod(const char name, const char title) :
14			BaseMod(name,title),
15	fabstoec	1.17	fPhotonBranchName (Names::gkPhotonBrn),
16			fGoodPhotonsName (ModNames::gkGoodPhotonsName),
17			fTrackBranchName (Names::gkTrackBrn),
18			fBeamspotBranchName(Names::gkBeamSpotBrn),
19			fPileUpDenName ("Rho"),
20	ceballos	1.5	fPhotonIDType("Custom"),
21	bendavid	1.16	fPhotonIsoType("Custom"),
22	sixie	1.1	fPhotonPtMin(15.0),
23	bendavid	1.16	fHadOverEmMax(0.02),
24			fApplySpikeRemoval(kFALSE),
25	ceballos	1.2	fApplyPixelSeed(kTRUE),
26	ceballos	1.10	fPhotonR9Min(0.5),
27	sixie	1.1	fPhIdType(kIdUndef),
28	loizides	1.7	fPhIsoType(kIsoUndef),
29	ceballos	1.10	fFiduciality(kTRUE),
30	bendavid	1.16	fEtaWidthEB(0.01),
31			fEtaWidthEE(0.028),
32	ceballos	1.11	fAbsEtaMax(2.5),
33	bendavid	1.16	fApplyR9Min(kFALSE),
34	fabstoec	1.17	fEffAreaEcal(0.139),
35			fEffAreaHcal(0.056),
36			fEffAreaTrack(0.296),
37			fPhotons(0),
38			fTracks(0),
39			fBeamspots(0),
40			fPileUpDen(0)
41
42	sixie	1.1	{
43			// Constructor.
44			}
45
46			//--------------------------------------------------------------------------------------------------
47			void PhotonIDMod::Process()
48			{
49			// Process entries of the tree.
50
51	fabstoec	1.17	LoadEventObject(fPhotonBranchName, fPhotons);
52			LoadEventObject(fTrackBranchName, fTracks);
53			LoadEventObject(fBeamspotBranchName, fBeamspots);
54			LoadEventObject(fPileUpDenName, fPileUpDen);
55
56			const BaseVertex *bsp = NULL;
57			if(fBeamspots->GetEntries() > 0)
58			bsp = dynamic_cast<const BaseVertex*>(fBeamspots->At(0));
59
60			Double_t _tRho = -1.;
61			if(fPileUpDen->GetEntries() > 0)
62			_tRho = (Double_t) fPileUpDen->At(0)->Rho();
63
64	sixie	1.1	PhotonOArr *GoodPhotons = new PhotonOArr;
65			GoodPhotons->SetName(fGoodPhotonsName);
66
67			for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
68			const Photon *ph = fPhotons->At(i);
69
70			if (ph->Pt() <= fPhotonPtMin)
71			continue;
72	sixie	1.13
73
74	sixie	1.1
75	sixie	1.13	Bool_t passSpikeRemovalFilter = kTRUE;
76	sixie	1.15
77			if (ph->SCluster() && ph->SCluster()->Seed()) {
78			if(ph->SCluster()->Seed()->Energy() > 5.0 &&
79			ph->SCluster()->Seed()->EMax() / ph->SCluster()->Seed()->E3x3() > 0.95
80			) {
81			passSpikeRemovalFilter = kFALSE;
82			}
83	sixie	1.13	}
84
85			// For Now Only use the EMax/E3x3 prescription.
86			//if(ph->SCluster()->Seed()->Energy() > 5.0 &&
87			// (1 - (ph->SCluster()->Seed()->E1x3() + ph->SCluster()->Seed()->E3x1() - 2*ph->SCluster()->Seed()->EMax())) > 0.95
88			// ) {
89			// passSpikeRemovalFilter = kFALSE;
90			//}
91			if (fApplySpikeRemoval && !passSpikeRemovalFilter) continue;
92
93
94
95	ceballos	1.2	if (ph->HadOverEm() >= fHadOverEmMax)
96			continue;
97	ceballos	1.14
98	ceballos	1.2	if (fApplyPixelSeed == kTRUE &&
99			ph->HasPixelSeed() == kTRUE)
100			continue;
101	ceballos	1.14
102	sixie	1.1	Bool_t idcut = kFALSE;
103			switch (fPhIdType) {
104			case kTight:
105			idcut = ph->IsTightPhoton();
106			break;
107			case kLoose:
108			idcut = ph->IsLoosePhoton();
109			break;
110			case kLooseEM:
111			idcut = ph->IsLooseEM();
112			break;
113			case kCustomId:
114	ceballos	1.5	idcut = kTRUE;
115	sixie	1.1	default:
116			break;
117			}
118
119			if (!idcut)
120			continue;
121
122			Bool_t isocut = kFALSE;
123			switch (fPhIsoType) {
124			case kNoIso:
125			isocut = kTRUE;
126			break;
127	ceballos	1.2	case kCombinedIso:
128	bendavid	1.12	{
129			Double_t totalIso = ph->HollowConeTrkIsoDr04()+
130			ph->EcalRecHitIsoDr04() +
131			ph->HcalTowerSumEtDr04();
132			if (totalIso/ph->Pt() < 0.25)
133			isocut = kTRUE;
134			}
135	ceballos	1.2	break;
136	sixie	1.1	case kCustomIso:
137	bendavid	1.16	{
138			if ( ph->HollowConeTrkIsoDr04() < (1.5 + 0.001ph->Pt()) && ph->EcalRecHitIsoDr04()<(2.0+0.006ph->Pt()) && ph->HcalTowerSumEtDr04()<(2.0+0.0025*ph->Pt()) )
139			isocut = kTRUE;
140			}
141	fabstoec	1.17	break;
142
143			case kMITPUCorrected:
144			{
145			// compute the PU corrections only if Rho is available
146			// ... otherwise (_tRho = -1.0) it's the std isolation
147			isocut = kTRUE;
148			Double_t EcalCorrISO = ph->EcalRecHitIsoDr04();
149			if(_tRho > -0.5 ) EcalCorrISO -= _tRho * fEffAreaEcal;
150			if ( EcalCorrISO > (2.0+0.006*ph->Pt()) ) isocut = kFALSE;
151			if ( isocut ) {
152			Double_t HcalCorrISO = ph->HcalTowerSumEtDr04();
153			if(_tRho > -0.5 ) HcalCorrISO -= _tRho * fEffAreaHcal;
154			if ( HcalCorrISO > (2.0+0.0025*ph->Pt()) ) isocut = kFALSE;
155			}
156			if ( isocut ) {
157			Double_t TrackCorrISO = IsolationTools::TrackIsolationNoPV(ph, bsp, 0.4, 0.04, 0.0, 0.015, 0.1, TrackQuality::highPurity, fTracks);
158			if(_tRho > -0.5 )
159			TrackCorrISO -= _tRho * fEffAreaTrack;
160			if ( TrackCorrISO > (1.5 + 0.001*ph->Pt()) ) isocut = kFALSE;
161			}
162			break;
163			}
164			default:
165			break;
166	sixie	1.1	}
167	fabstoec	1.17
168	sixie	1.1	if (!isocut)
169			continue;
170	fabstoec	1.17
171	bendavid	1.16	if ( fApplyR9Min && ph->R9() <= fPhotonR9Min)
172	ceballos	1.5	continue;
173
174	ceballos	1.10	if (fFiduciality == kTRUE &&
175			ph->IsEB() == kFALSE && ph->IsEE() == kFALSE)
176			continue;
177
178	bendavid	1.16	if ((ph->IsEB() == kTRUE && ph->CoviEtaiEta() >= fEtaWidthEB) \|\|
179			(ph->IsEE() == kTRUE && ph->CoviEtaiEta() >= fEtaWidthEE))
180	ceballos	1.10	continue;
181
182	ceballos	1.11	if (ph->AbsEta() >= fAbsEtaMax)
183			continue;
184	ceballos	1.14
185	sixie	1.1	// add good electron
186			GoodPhotons->Add(fPhotons->At(i));
187			}
188
189	loizides	1.4	// sort according to pt
190			GoodPhotons->Sort();
191
192	sixie	1.1	// add to event for other modules to use
193			AddObjThisEvt(GoodPhotons);
194			}
195
196			//--------------------------------------------------------------------------------------------------
197			void PhotonIDMod::SlaveBegin()
198			{
199			// Run startup code on the computer (slave) doing the actual analysis. Here,
200	loizides	1.3	// we just request the photon collection branch.
201	sixie	1.1
202	fabstoec	1.17	ReqEventObject(fPhotonBranchName, fPhotons, kTRUE);
203			ReqEventObject(fTrackBranchName, fTracks, kTRUE);
204			ReqEventObject(fBeamspotBranchName, fBeamspots, kTRUE);
205			ReqEventObject(fPileUpDenName, fPileUpDen, kTRUE);
206	sixie	1.1
207			if (fPhotonIDType.CompareTo("Tight") == 0)
208			fPhIdType = kTight;
209			else if (fPhotonIDType.CompareTo("Loose") == 0)
210			fPhIdType = kLoose;
211			else if (fPhotonIDType.CompareTo("LooseEM") == 0)
212			fPhIdType = kLooseEM;
213			else if (fPhotonIDType.CompareTo("Custom") == 0) {
214			fPhIdType = kCustomId;
215			} else {
216			SendError(kAbortAnalysis, "SlaveBegin",
217	ceballos	1.2	"The specified photon identification %s is not defined.",
218	sixie	1.1	fPhotonIDType.Data());
219			return;
220			}
221
222			if (fPhotonIsoType.CompareTo("NoIso") == 0 )
223			fPhIsoType = kNoIso;
224	ceballos	1.2	else if (fPhotonIsoType.CompareTo("CombinedIso") == 0 )
225			fPhIsoType = kCombinedIso;
226	fabstoec	1.17	else if (fPhotonIsoType.CompareTo("Custom") == 0 )
227	sixie	1.1	fPhIsoType = kCustomIso;
228	fabstoec	1.17	else if (fPhotonIsoType.CompareTo("MITPUCorrected") == 0 )
229			fPhIsoType = kMITPUCorrected;
230			else {
231	sixie	1.1	SendError(kAbortAnalysis, "SlaveBegin",
232	ceballos	1.2	"The specified photon isolation %s is not defined.",
233	sixie	1.1	fPhotonIsoType.Data());
234			return;
235			}
236			}