ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonIDMod.cc
Revision: 1.24
Committed: Wed Jul 27 15:17:37 2011 UTC (13 years, 9 months ago) by bendavid
Content type: text/plain
Branch: MAIN
Changes since 1.23: +134 -12 lines
Log Message:
update photon stuff

File Contents

# User Rev Content
1 bendavid 1.24 // $Id: PhotonIDMod.cc,v 1.23 2011/07/22 10:58:08 fabstoec 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     using namespace mithep;
12    
13     ClassImp(mithep::PhotonIDMod)
14    
15     //--------------------------------------------------------------------------------------------------
16     PhotonIDMod::PhotonIDMod(const char *name, const char *title) :
17     BaseMod(name,title),
18 fabstoec 1.17 fPhotonBranchName (Names::gkPhotonBrn),
19     fGoodPhotonsName (ModNames::gkGoodPhotonsName),
20     fTrackBranchName (Names::gkTrackBrn),
21     fBeamspotBranchName(Names::gkBeamSpotBrn),
22 fabstoec 1.18 fPileUpDenName (Names::gkPileupEnergyDensityBrn),
23 bendavid 1.19 fConversionName ("MergedConversions"),
24     fElectronName ("Electrons"),
25 fabstoec 1.22 fPVName (Names::gkPVBeamSpotBrn),
26 bendavid 1.24 // MC specific stuff...
27     fMCParticleName (Names::gkMCPartBrn),
28     fPileUpName (Names::gkPileupInfoBrn),
29 fabstoec 1.22
30     fPhotonIDType ("Custom"),
31     fPhotonIsoType ("Custom"),
32     fPhotonPtMin (15.0),
33     fHadOverEmMax (0.02),
34     fApplySpikeRemoval (kFALSE),
35     fApplyPixelSeed (kTRUE),
36     fApplyElectronVeto (kFALSE),
37 bendavid 1.19 fApplyElectronVetoConvRecovery(kFALSE),
38 fabstoec 1.22 fApplyConversionId (kFALSE),
39 bendavid 1.19 fApplyTriggerMatching(kFALSE),
40 fabstoec 1.22 fPhotonR9Min (0.5),
41     fPhIdType (kIdUndef),
42     fPhIsoType (kIsoUndef),
43     fFiduciality (kTRUE),
44     fEtaWidthEB (0.01),
45     fEtaWidthEE (0.028),
46     fAbsEtaMax (999.99),
47     fApplyR9Min (kFALSE),
48    
49     fEffAreaEcalEE (0.089),
50     fEffAreaHcalEE (0.156),
51     fEffAreaTrackEE (0.261),
52    
53     fEffAreaEcalEB (0.183),
54     fEffAreaHcalEB (0.062),
55     fEffAreaTrackEB (0.306),
56    
57 fabstoec 1.17 fPhotons(0),
58     fTracks(0),
59     fBeamspots(0),
60 bendavid 1.19 fPileUpDen(0),
61     fConversions(0),
62 fabstoec 1.22 fElectrons(0),
63     fPV(0),
64 bendavid 1.24 fMCParticles (0),
65     fPileUp (0),
66 fabstoec 1.22
67 bendavid 1.24 fPVFromBranch (true),
68     fIsData(false),
69     fWriteTree(false)
70 fabstoec 1.17
71 sixie 1.1 {
72     // Constructor.
73     }
74    
75     //--------------------------------------------------------------------------------------------------
76     void PhotonIDMod::Process()
77     {
78     // Process entries of the tree.
79    
80 fabstoec 1.17 LoadEventObject(fPhotonBranchName, fPhotons);
81 bendavid 1.19
82     const BaseVertex *bsp = NULL;
83     Double_t _tRho = -1.;
84     const TriggerObjectCol *trigObjs = 0;
85     if (fPhotons->GetEntries()>0) {
86     LoadEventObject(fTrackBranchName, fTracks);
87     LoadEventObject(fBeamspotBranchName, fBeamspots);
88 bendavid 1.24 LoadEventObject(fPileUpDenName, fPileUpDen);
89 bendavid 1.19 LoadEventObject(fConversionName, fConversions);
90     LoadEventObject(fElectronName, fElectrons);
91 bendavid 1.24 LoadEventObject(fPVName, fPV);
92 bendavid 1.19
93     if(fBeamspots->GetEntries() > 0)
94     bsp = fBeamspots->At(0);
95 fabstoec 1.17
96 bendavid 1.24 if(fPileUpDen->GetEntries() > 0)
97 fabstoec 1.23 _tRho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
98 bendavid 1.19
99     //get trigger object collection if trigger matching is enabled
100     if (fApplyTriggerMatching) {
101     trigObjs = GetHLTObjects(fTrigObjectsName);
102     }
103    
104     }
105 fabstoec 1.17
106 sixie 1.1 PhotonOArr *GoodPhotons = new PhotonOArr;
107     GoodPhotons->SetName(fGoodPhotonsName);
108    
109     for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
110     const Photon *ph = fPhotons->At(i);
111    
112 bendavid 1.24
113     if (fFiduciality == kTRUE &&
114     (ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566) ) )
115     continue;
116    
117 fabstoec 1.22 // ---------------------------------------------------------------------
118     // check if we use the CiC Selection. If yes, bypass all the below...
119     if(fPhIdType == kBaseLineCiC) {
120     if( PhotonTools::PassCiCSelection(ph, fPV->At(0), fTracks, fElectrons, fPV, _tRho, fPhotonPtMin, fApplyElectronVeto) )
121     GoodPhotons->Add(fPhotons->At(i));
122     continue; // go to next Photons
123     }
124     // ---------------------------------------------------------------------
125    
126 sixie 1.1 if (ph->Pt() <= fPhotonPtMin)
127 fabstoec 1.22 continue; // add good electron
128 sixie 1.13
129 bendavid 1.21 Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
130 sixie 1.1
131 sixie 1.13 Bool_t passSpikeRemovalFilter = kTRUE;
132 sixie 1.15
133     if (ph->SCluster() && ph->SCluster()->Seed()) {
134     if(ph->SCluster()->Seed()->Energy() > 5.0 &&
135     ph->SCluster()->Seed()->EMax() / ph->SCluster()->Seed()->E3x3() > 0.95
136     ) {
137     passSpikeRemovalFilter = kFALSE;
138     }
139 sixie 1.13 }
140    
141     // For Now Only use the EMax/E3x3 prescription.
142     //if(ph->SCluster()->Seed()->Energy() > 5.0 &&
143     // (1 - (ph->SCluster()->Seed()->E1x3() + ph->SCluster()->Seed()->E3x1() - 2*ph->SCluster()->Seed()->EMax())) > 0.95
144     // ) {
145     // passSpikeRemovalFilter = kFALSE;
146     //}
147     if (fApplySpikeRemoval && !passSpikeRemovalFilter) continue;
148    
149 ceballos 1.2 if (ph->HadOverEm() >= fHadOverEmMax)
150     continue;
151 ceballos 1.14
152 ceballos 1.2 if (fApplyPixelSeed == kTRUE &&
153     ph->HasPixelSeed() == kTRUE)
154     continue;
155 ceballos 1.14
156 bendavid 1.19 if (fApplyElectronVeto && !PhotonTools::PassElectronVeto(ph,fElectrons) ) continue;
157    
158     if (fApplyElectronVetoConvRecovery && !PhotonTools::PassElectronVetoConvRecovery(ph,fElectrons,fConversions,bsp) ) continue;
159    
160     if (fApplyConversionId && !PhotonTools::PassConversionId(ph,PhotonTools::MatchedConversion(ph,fConversions,bsp))) continue;
161    
162     if (fApplyTriggerMatching && !PhotonTools::PassTriggerMatching(ph,trigObjs)) continue;
163    
164 sixie 1.1 Bool_t idcut = kFALSE;
165     switch (fPhIdType) {
166     case kTight:
167     idcut = ph->IsTightPhoton();
168     break;
169     case kLoose:
170     idcut = ph->IsLoosePhoton();
171     break;
172     case kLooseEM:
173     idcut = ph->IsLooseEM();
174     break;
175     case kCustomId:
176 ceballos 1.5 idcut = kTRUE;
177 sixie 1.1 default:
178     break;
179     }
180    
181     if (!idcut)
182     continue;
183    
184     Bool_t isocut = kFALSE;
185     switch (fPhIsoType) {
186     case kNoIso:
187     isocut = kTRUE;
188     break;
189 ceballos 1.2 case kCombinedIso:
190 bendavid 1.12 {
191     Double_t totalIso = ph->HollowConeTrkIsoDr04()+
192     ph->EcalRecHitIsoDr04() +
193     ph->HcalTowerSumEtDr04();
194     if (totalIso/ph->Pt() < 0.25)
195     isocut = kTRUE;
196     }
197 ceballos 1.2 break;
198 sixie 1.1 case kCustomIso:
199 bendavid 1.16 {
200     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()) )
201     isocut = kTRUE;
202     }
203 fabstoec 1.17 break;
204    
205     case kMITPUCorrected:
206     {
207     // compute the PU corrections only if Rho is available
208     // ... otherwise (_tRho = -1.0) it's the std isolation
209     isocut = kTRUE;
210 fabstoec 1.20 Double_t fEffAreaEcal = fEffAreaEcalEB;
211     Double_t fEffAreaHcal = fEffAreaHcalEB;
212     Double_t fEffAreaTrack = fEffAreaTrackEB;
213    
214 bendavid 1.21 if( !isbarrel ) {
215 fabstoec 1.20 fEffAreaEcal = fEffAreaEcalEE;
216     fEffAreaHcal = fEffAreaHcalEE;
217     fEffAreaTrack = fEffAreaTrackEE;
218     }
219    
220 fabstoec 1.17 Double_t EcalCorrISO = ph->EcalRecHitIsoDr04();
221     if(_tRho > -0.5 ) EcalCorrISO -= _tRho * fEffAreaEcal;
222     if ( EcalCorrISO > (2.0+0.006*ph->Pt()) ) isocut = kFALSE;
223     if ( isocut ) {
224     Double_t HcalCorrISO = ph->HcalTowerSumEtDr04();
225     if(_tRho > -0.5 ) HcalCorrISO -= _tRho * fEffAreaHcal;
226     if ( HcalCorrISO > (2.0+0.0025*ph->Pt()) ) isocut = kFALSE;
227     }
228     if ( isocut ) {
229     Double_t TrackCorrISO = IsolationTools::TrackIsolationNoPV(ph, bsp, 0.4, 0.04, 0.0, 0.015, 0.1, TrackQuality::highPurity, fTracks);
230     if(_tRho > -0.5 )
231     TrackCorrISO -= _tRho * fEffAreaTrack;
232     if ( TrackCorrISO > (1.5 + 0.001*ph->Pt()) ) isocut = kFALSE;
233     }
234     break;
235     }
236     default:
237     break;
238 sixie 1.1 }
239 fabstoec 1.17
240 sixie 1.1 if (!isocut)
241     continue;
242 fabstoec 1.17
243 bendavid 1.16 if ( fApplyR9Min && ph->R9() <= fPhotonR9Min)
244 ceballos 1.5 continue;
245    
246 bendavid 1.21 if ((isbarrel && ph->CoviEtaiEta() >= fEtaWidthEB) ||
247     (!isbarrel && ph->CoviEtaiEta() >= fEtaWidthEE))
248 ceballos 1.10 continue;
249    
250 ceballos 1.11 if (ph->AbsEta() >= fAbsEtaMax)
251     continue;
252 ceballos 1.14
253 sixie 1.1 // add good electron
254     GoodPhotons->Add(fPhotons->At(i));
255     }
256    
257 loizides 1.4 // sort according to pt
258     GoodPhotons->Sort();
259    
260 sixie 1.1 // add to event for other modules to use
261     AddObjThisEvt(GoodPhotons);
262 bendavid 1.24
263     if (!fWriteTree || !GoodPhotons->GetEntries()) return;
264    
265     Int_t _numPU = -1.; // some sensible default values....
266     Int_t _numPUminus = -1.; // some sensible default values....
267     Int_t _numPUplus = -1.; // some sensible default values....
268    
269     if( !fIsData ) {
270     LoadBranch(fMCParticleName);
271     LoadBranch(fPileUpName);
272     for (UInt_t i=0; i<fPileUp->GetEntries(); ++i) {
273     const PileupInfo *puinfo = fPileUp->At(i);
274     if (puinfo->GetBunchCrossing()==0) _numPU = puinfo->GetPU_NumInteractions();
275     else if (puinfo->GetBunchCrossing() == -1) _numPUminus = puinfo->GetPU_NumInteractions();
276     else if (puinfo->GetBunchCrossing() == 1) _numPUplus = puinfo->GetPU_NumInteractions();
277     }
278     }
279    
280    
281     fDiphotonEvent->rho = _tRho;
282     fDiphotonEvent->genHiggspt = -99.;
283     fDiphotonEvent->genHiggsZ = -99.;
284     fDiphotonEvent->gencostheta = -99.;
285     fDiphotonEvent->nVtx = fPV->GetEntries();
286     fDiphotonEvent->vtxZ = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->Z() : -99.;
287     fDiphotonEvent->numPU = _numPU;
288     fDiphotonEvent->numPUminus = _numPUminus;
289     fDiphotonEvent->numPUplus = _numPUplus;
290     fDiphotonEvent->mass = -99.;
291     fDiphotonEvent->ptgg = -99.;
292     fDiphotonEvent->costheta = -99.;
293     fDiphotonEvent->evt = GetEventHeader()->EvtNum();
294     fDiphotonEvent->run = GetEventHeader()->RunNum();
295     fDiphotonEvent->lumi = GetEventHeader()->LumiSec();
296     fDiphotonEvent->evtcat = -99;
297     fDiphotonEvent->masscor = -99.;
298     fDiphotonEvent->masscorerr = -99.;
299    
300     for (UInt_t i=0; i<GoodPhotons->GetEntries(); ++i) {
301     const Photon *ph = GoodPhotons->At(i);
302     const MCParticle *phgen = PhotonTools::MatchMC(ph,fMCParticles);
303     fSinglePhoton->SetVars(ph,phgen);
304     hPhotonTree->Fill();
305     }
306    
307    
308 sixie 1.1 }
309    
310     //--------------------------------------------------------------------------------------------------
311     void PhotonIDMod::SlaveBegin()
312     {
313     // Run startup code on the computer (slave) doing the actual analysis. Here,
314 loizides 1.3 // we just request the photon collection branch.
315 sixie 1.1
316 fabstoec 1.22 ReqEventObject(fPhotonBranchName, fPhotons, kTRUE);
317     ReqEventObject(fTrackBranchName, fTracks, kTRUE);
318     ReqEventObject(fBeamspotBranchName, fBeamspots, kTRUE);
319 bendavid 1.19 ReqEventObject(fConversionName, fConversions, kTRUE);
320 fabstoec 1.22 ReqEventObject(fElectronName, fElectrons, kTRUE);
321 bendavid 1.24 ReqEventObject(fPVName, fPV, fPVFromBranch);
322     ReqEventObject(fPileUpDenName, fPileUpDen, kTRUE);
323     if (!fIsData) {
324     ReqBranch(fPileUpName, fPileUp);
325     ReqBranch(fMCParticleName, fMCParticles);
326     }
327    
328 sixie 1.1 if (fPhotonIDType.CompareTo("Tight") == 0)
329     fPhIdType = kTight;
330     else if (fPhotonIDType.CompareTo("Loose") == 0)
331     fPhIdType = kLoose;
332     else if (fPhotonIDType.CompareTo("LooseEM") == 0)
333     fPhIdType = kLooseEM;
334 fabstoec 1.22 else if (fPhotonIDType.CompareTo("Custom") == 0)
335 sixie 1.1 fPhIdType = kCustomId;
336 fabstoec 1.22 else if (fPhotonIDType.CompareTo("BaseLineCiC") == 0) {
337     fPhIdType = kBaseLineCiC;
338     fPhotonIsoType = "NoIso";
339 sixie 1.1 } else {
340     SendError(kAbortAnalysis, "SlaveBegin",
341 ceballos 1.2 "The specified photon identification %s is not defined.",
342 sixie 1.1 fPhotonIDType.Data());
343     return;
344     }
345    
346     if (fPhotonIsoType.CompareTo("NoIso") == 0 )
347     fPhIsoType = kNoIso;
348 ceballos 1.2 else if (fPhotonIsoType.CompareTo("CombinedIso") == 0 )
349     fPhIsoType = kCombinedIso;
350 fabstoec 1.17 else if (fPhotonIsoType.CompareTo("Custom") == 0 )
351 sixie 1.1 fPhIsoType = kCustomIso;
352 fabstoec 1.17 else if (fPhotonIsoType.CompareTo("MITPUCorrected") == 0 )
353     fPhIsoType = kMITPUCorrected;
354     else {
355 sixie 1.1 SendError(kAbortAnalysis, "SlaveBegin",
356 ceballos 1.2 "The specified photon isolation %s is not defined.",
357 sixie 1.1 fPhotonIsoType.Data());
358     return;
359     }
360 bendavid 1.21
361 bendavid 1.24 if (fWriteTree) {
362     fDiphotonEvent = new PhotonPairSelectorDiphotonEvent;
363     fSinglePhoton = new PhotonPairSelectorPhoton;
364 bendavid 1.21
365 bendavid 1.24 TString singlename = "hPhotonTree";
366     hPhotonTree = new TTree(singlename,singlename);
367    
368     //make flattish tree from classes so we don't have to rely on dictionaries for reading later
369     TClass *eclass = TClass::GetClass("mithep::PhotonPairSelectorDiphotonEvent");
370     TClass *pclass = TClass::GetClass("mithep::PhotonPairSelectorPhoton");
371     TList *elist = eclass->GetListOfDataMembers();
372     TList *plist = pclass->GetListOfDataMembers();
373    
374     for (int i=0; i<elist->GetEntries(); ++i) {
375     const TDataMember *tdm = static_cast<const TDataMember*>(elist->At(i));
376     if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
377     TString typestring;
378     if (TString(tdm->GetTypeName())=="Char_t") typestring = "B";
379     else if (TString(tdm->GetTypeName())=="UChar_t") typestring = "b";
380     else if (TString(tdm->GetTypeName())=="Short_t") typestring = "S";
381     else if (TString(tdm->GetTypeName())=="UShort_t") typestring = "s";
382     else if (TString(tdm->GetTypeName())=="Int_t") typestring = "I";
383     else if (TString(tdm->GetTypeName())=="UInt_t") typestring = "i";
384     else if (TString(tdm->GetTypeName())=="Float_t") typestring = "F";
385     else if (TString(tdm->GetTypeName())=="Double_t") typestring = "D";
386     else if (TString(tdm->GetTypeName())=="Long64_t") typestring = "L";
387     else if (TString(tdm->GetTypeName())=="ULong64_t") typestring = "l";
388     else if (TString(tdm->GetTypeName())=="Bool_t") typestring = "O";
389     else continue;
390     //printf("%s %s: %i\n",tdm->GetTypeName(),tdm->GetName(),int(tdm->GetOffset()));
391     Char_t *addr = (Char_t*)fDiphotonEvent;
392     assert(sizeof(Char_t)==1);
393     hPhotonTree->Branch(tdm->GetName(),addr + tdm->GetOffset(),TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
394     }
395 bendavid 1.21
396 bendavid 1.24
397     for (int i=0; i<plist->GetEntries(); ++i) {
398     const TDataMember *tdm = static_cast<const TDataMember*>(plist->At(i));
399     if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
400     TString typestring;
401     if (TString(tdm->GetTypeName())=="Char_t") typestring = "B";
402     else if (TString(tdm->GetTypeName())=="UChar_t") typestring = "b";
403     else if (TString(tdm->GetTypeName())=="Short_t") typestring = "S";
404     else if (TString(tdm->GetTypeName())=="UShort_t") typestring = "s";
405     else if (TString(tdm->GetTypeName())=="Int_t") typestring = "I";
406     else if (TString(tdm->GetTypeName())=="UInt_t") typestring = "i";
407     else if (TString(tdm->GetTypeName())=="Float_t") typestring = "F";
408     else if (TString(tdm->GetTypeName())=="Double_t") typestring = "D";
409     else if (TString(tdm->GetTypeName())=="Long64_t") typestring = "L";
410     else if (TString(tdm->GetTypeName())=="ULong64_t") typestring = "l";
411     else if (TString(tdm->GetTypeName())=="Bool_t") typestring = "O";
412     else continue;
413     //printf("%s\n",tdm->GetTypeName());
414    
415     assert(sizeof(Char_t)==1);
416     TString singlename = TString::Format("ph.%s",tdm->GetName());
417     Char_t *addrsingle = (Char_t*)fSinglePhoton;
418     hPhotonTree->Branch(singlename,addrsingle+tdm->GetOffset(),TString::Format("%s/%s",singlename.Data(),typestring.Data()));
419     }
420    
421     AddOutput(hPhotonTree);
422     }
423    
424    
425 sixie 1.1 }