ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonPairSelector.cc
Revision: 1.8
Committed: Wed Aug 3 17:15:43 2011 UTC (13 years, 9 months ago) by bendavid
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_024b, Mit_024a, Mit_024
Changes since 1.7: +42 -300 lines
Log Message:
Reorganize photon code and add new PhotonTreeWriter class

File Contents

# User Rev Content
1 fabstoec 1.1 #include "MitPhysics/Mods/interface/PhotonPairSelector.h"
2     #include "MitAna/DataTree/interface/PhotonCol.h"
3     #include "MitAna/DataTree/interface/PFCandidateCol.h"
4 bendavid 1.8 #include "MitAna/DataTree/interface/StableData.h"
5     #include "MitAna/DataTree/interface/StableParticle.h"
6 fabstoec 1.1 #include "MitPhysics/Init/interface/ModNames.h"
7     #include "MitPhysics/Utils/interface/IsolationTools.h"
8     #include "MitPhysics/Utils/interface/PhotonTools.h"
9     #include "MitPhysics/Utils/interface/VertexTools.h"
10 bendavid 1.6 #include "MitPhysics/Utils/interface/PhotonFix.h"
11 bendavid 1.5 #include "TDataMember.h"
12 fabstoec 1.1 #include <TNtuple.h>
13     #include <TRandom3.h>
14 bendavid 1.6 #include <TSystem.h>
15 fabstoec 1.1
16     using namespace mithep;
17    
18     ClassImp(mithep::PhotonPairSelector)
19    
20     //--------------------------------------------------------------------------------------------------
21     PhotonPairSelector::PhotonPairSelector(const char *name, const char *title) :
22     // Base Module...
23     BaseMod (name,title),
24    
25     // define all the Branches to load
26     fPhotonBranchName (Names::gkPhotonBrn),
27     fElectronName (Names::gkElectronBrn),
28 bendavid 1.8 fGoodElectronName (Names::gkElectronBrn),
29 fabstoec 1.1 fConversionName (Names::gkMvfConversionBrn),
30     fTrackBranchName (Names::gkTrackBrn),
31     fPileUpDenName (Names::gkPileupEnergyDensityBrn),
32     fPVName (Names::gkPVBeamSpotBrn),
33     fBeamspotName (Names::gkBeamSpotBrn),
34     fPFCandName (Names::gkPFCandidatesBrn),
35     // MC specific stuff...
36     fMCParticleName (Names::gkMCPartBrn),
37     fPileUpName (Names::gkPileupInfoBrn),
38    
39     fGoodPhotonsName (ModNames::gkGoodPhotonsName),
40    
41     // ----------------------------------------
42     // Selection Types
43     fPhotonSelType ("NoSelection"),
44     fVertexSelType ("StdSelection"),
45     fPhSelType (kNoPhSelection),
46     fVtxSelType (kStdVtxSelection),
47    
48     // ----------------------------------------
49     fPhotonPtMin (20.0),
50     fPhotonEtaMax (2.5),
51    
52     fLeadingPtMin (40.0),
53     fTrailingPtMin (30.0),
54    
55     fIsData (false),
56     fPhotonsFromBranch (true),
57     fPVFromBranch (true),
58 bendavid 1.8 fGoodElectronsFromBranch (kTRUE),
59    
60 fabstoec 1.1 // ----------------------------------------
61     // collections....
62     fPhotons (0),
63     fElectrons (0),
64     fConversions (0),
65     fTracks (0),
66     fPileUpDen (0),
67     fPV (0),
68     fBeamspot (0),
69     fPFCands (0),
70     fMCParticles (0),
71     fPileUp (0),
72    
73     // ---------------------------------------
74     fDataEnCorr_EB_hR9 (0.),
75     fDataEnCorr_EB_lR9 (0.),
76     fDataEnCorr_EE_hR9 (0.),
77     fDataEnCorr_EE_lR9 (0.),
78    
79     fRunStart (0),
80     fRunEnd (0),
81    
82     fMCSmear_EB_hR9 (0.),
83     fMCSmear_EB_lR9 (0.),
84     fMCSmear_EE_hR9 (0.),
85     fMCSmear_EE_lR9 (0.),
86    
87     // ---------------------------------------
88     rng (new TRandom3()),
89    
90     // ---------------------------------------
91     fDoDataEneCorr (true),
92     fDoMCSmear (true),
93 fabstoec 1.2 fDoVtxSelection (true),
94     fApplyEleVeto (true),
95 bendavid 1.8 fInvertElectronVeto(kFALSE)
96 fabstoec 1.1 {
97     // Constructor.
98     }
99    
100     PhotonPairSelector::~PhotonPairSelector(){
101     if(rng) delete rng;
102     }
103    
104     //--------------------------------------------------------------------------------------------------
105     void PhotonPairSelector::Process()
106     {
107     // ------------------------------------------------------------
108     // Process entries of the tree.
109     LoadEventObject(fPhotonBranchName, fPhotons);
110 bendavid 1.8
111     // -----------------------------------------------------------
112     // OUtput Photon Collection. It will ALWAYS conatrin either 0 or 2 Photons
113     PhotonOArr *GoodPhotons = new PhotonOArr;
114     GoodPhotons->SetName(fGoodPhotonsName);
115     GoodPhotons->SetOwner(kTRUE);
116     // add to event for other modules to use
117     AddObjThisEvt(GoodPhotons);
118    
119     if (fPhotons->GetEntries()<2) return;
120    
121 fabstoec 1.1 LoadEventObject(fElectronName, fElectrons);
122 bendavid 1.8 LoadEventObject(fGoodElectronName, fGoodElectrons);
123 fabstoec 1.1 LoadEventObject(fConversionName, fConversions);
124     LoadEventObject(fTrackBranchName, fTracks);
125     LoadEventObject(fPileUpDenName, fPileUpDen);
126     LoadEventObject(fPVName, fPV);
127     LoadEventObject(fBeamspotName, fBeamspot);
128     LoadEventObject(fPFCandName, fPFCands);
129    
130 bendavid 1.6
131 fabstoec 1.1
132 bendavid 1.8
133 fabstoec 1.1
134    
135     // ------------------------------------------------------------
136     // load event based information
137 bendavid 1.4 Int_t _numPU = -1.; // some sensible default values....
138     Int_t _numPUminus = -1.; // some sensible default values....
139     Int_t _numPUplus = -1.; // some sensible default values....
140    
141 fabstoec 1.1 Float_t _tRho = -99.;
142 fabstoec 1.2 if( fPileUpDen->GetEntries() > 0 )
143     _tRho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
144 bendavid 1.6
145    
146 fabstoec 1.2
147 fabstoec 1.1 const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));
148    
149     // ------------------------------------------------------------
150     // Get Event header for Run info etc.
151     const EventHeader* evtHead = this->GetEventHeader();
152     unsigned int evtNum = evtHead->EvtNum();
153 fabstoec 1.7 //Float_t _evtNum1 = (Float_t) ( (int) (evtNum/10000.) );
154     //Float_t _evtNum2 = (Float_t) ( (int) (evtNum % 10000) );
155 fabstoec 1.1 UInt_t runNumber = evtHead->RunNum();
156     Float_t _runNum = (Float_t) runNumber;
157     Float_t _lumiSec = (Float_t) evtHead->LumiSec();
158    
159     // ------------------------------------------------------------
160     // here we'll store the preselected Photons (and which CiCCategory they are...)
161     PhotonOArr* preselPh = new PhotonOArr;
162     std::vector<PhotonTools::CiCBaseLineCats> preselCat;
163     preselCat.resize(0);
164    
165     // 1. we do the pre-selection; but keep the non-passing photons in a secont container...
166     for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
167     const Photon *ph = fPhotons->At(i);
168    
169     if(ph->SCluster()->AbsEta()>= fPhotonEtaMax || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) continue;
170     if(ph->Et() < fPhotonPtMin) continue;
171     if(ph->HadOverEm() > 0.15) continue;
172     if(ph->IsEB()) {
173     if(ph->CoviEtaiEta() > 0.013) continue;
174     } else {
175     if(ph->CoviEtaiEta() > 0.03) continue;
176     }
177     preselPh->Add(ph);
178 fabstoec 1.2 }
179    
180 bendavid 1.8 if (preselPh->GetEntries()<2) return;
181    
182 fabstoec 1.2 // Sorry... need the second loop here in order to sort & assign the right Categories..
183     preselPh->Sort();
184     for(unsigned int iPh = 0; iPh <preselPh->GetEntries(); ++iPh) {
185     const Photon* ph = preselPh->At(iPh);
186 fabstoec 1.1 preselCat.push_back(PhotonTools::CiCBaseLineCat(ph));
187     }
188    
189     // ------------------------------------------------------------
190     // compute how many pairs there are ...
191     unsigned int numPairs = 0;
192     if( preselPh->GetEntries() > 0) numPairs = (preselPh->GetEntries()-1)*preselPh->GetEntries()/2;
193     // ... and create all possible pairs of pre-selected photons
194     std::vector<unsigned int> idx1st;
195     std::vector<unsigned int> idx2nd;
196     std::vector<PhotonTools::CiCBaseLineCats> cat1st;
197     std::vector<PhotonTools::CiCBaseLineCats> cat2nd;
198     // ... this will be used to store whether a givne pair passes the cuts
199     std::vector<bool> pairPasses;
200    
201     if(numPairs > 0) {
202     for(unsigned int i1st = 0; i1st <preselPh->GetEntries() - 1; ++i1st) {
203     for(unsigned int i2nd = i1st + 1; i2nd <preselPh->GetEntries(); ++i2nd) {
204     idx1st.push_back(i1st);
205     idx2nd.push_back(i2nd);
206     pairPasses.push_back(true);
207     }
208     }
209     }
210    
211     // ------------------------------------------------------------
212     // array to store the index of 'chosen Vtx' for each pair
213     const Vertex** theVtx = new const Vertex*[numPairs]; // holds the 'chosen' Vtx for each Pair
214     Photon** fixPh1st = new Photon*[numPairs]; // holds the 1st Photon for each Pair
215     Photon** fixPh2nd = new Photon*[numPairs]; // holds the 2nd photon for each Pair
216    
217    
218     // store pair-indices for pairs passing the selection
219     std::vector<unsigned int> passPairs;
220     passPairs.resize(0);
221    
222     // ------------------------------------------------------------
223     // Loop over all Pairs and to the 'incredible machine' running....
224     for(unsigned int iPair = 0; iPair < numPairs; ++iPair) {
225    
226     // first we need a hard copy of the incoming photons
227     fixPh1st[iPair] = new Photon(*preselPh->At(idx1st[iPair]));
228     fixPh2nd[iPair] = new Photon(*preselPh->At(idx2nd[iPair]));
229     // we also store the category, so we don't have to ask all the time...
230     cat1st.push_back(preselCat[idx1st[iPair]]);
231     cat2nd.push_back(preselCat[idx2nd[iPair]]);
232    
233     // now we dicide if we either scale (Data) or Smear (MC) the Photons
234     if (fIsData) {
235     if(fDoDataEneCorr) {
236     // statring with scale = 1.
237     double scaleFac1 = 1.;
238     double scaleFac2 = 1.;
239     // checking the run Rangees ...
240     Int_t runRange = FindRunRangeIdx(runNumber);
241     if(runRange > -1) {
242     scaleFac1 += GetDataEnCorr(runRange, cat1st[iPair]);
243     scaleFac2 += GetDataEnCorr(runRange, cat2nd[iPair]);
244     }
245     PhotonTools::ScalePhoton(fixPh1st[iPair], scaleFac1);
246     PhotonTools::ScalePhoton(fixPh2nd[iPair], scaleFac2);
247     }
248     } else {
249     if(fDoMCSmear) {
250     // get the seed to do deterministic smearing...
251     UInt_t seedBase = (UInt_t) evtNum + (UInt_t) _runNum + (UInt_t) _lumiSec;
252     UInt_t seed1 = seedBase + (UInt_t) fixPh1st[iPair]->E() + (UInt_t) (TMath::Abs(10.*fixPh1st[iPair]->SCluster()->Eta()));
253     UInt_t seed2 = seedBase + (UInt_t) fixPh2nd[iPair]->E() + (UInt_t) (TMath::Abs(10.*fixPh2nd[iPair]->SCluster()->Eta()));
254     // get the smearing for MC photons..
255 fabstoec 1.3
256 fabstoec 1.1 double width1 = GetMCSmearFac(cat1st[iPair]);
257     double width2 = GetMCSmearFac(cat2nd[iPair]);
258 fabstoec 1.3
259 fabstoec 1.1 PhotonTools::SmearPhoton(fixPh1st[iPair], rng, width1, seed1);
260     PhotonTools::SmearPhoton(fixPh2nd[iPair], rng, width2, seed2);
261 fabstoec 1.3
262 fabstoec 1.1 }
263     }
264    
265     // store the vertex for this pair
266     switch( fVtxSelType ){
267     case kStdVtxSelection:
268     theVtx[iPair] = fPV->At(0);
269     break;
270     case kCiCVtxSelection:
271     theVtx[iPair] = VertexTools::findVtxBasicRanking(fixPh1st[iPair],fixPh2nd[iPair], bsp, fPV, fConversions);
272     break;
273     case kMITVtxSelection:
274     // need PFCandidate Collection
275     theVtx[iPair] = VertexTools::BestVtx(fPFCands, fPV, bsp, mithep::FourVector((fixPh1st[iPair]->Mom()+fixPh2nd[iPair]->Mom())));
276     break;
277     default:
278     theVtx[iPair] = fPV->At(0);
279 fabstoec 1.2
280 fabstoec 1.1 }
281 bendavid 1.8
282     //set PV ref in photons
283     fixPh1st[iPair]->SetPV(theVtx[iPair]);
284     fixPh2nd[iPair]->SetPV(theVtx[iPair]);
285 fabstoec 1.2
286 fabstoec 1.1 // fix the kinematics for both events
287     FourVectorM newMom1st = fixPh1st[iPair]->MomVtx(theVtx[iPair]->Position());
288     FourVectorM newMom2nd = fixPh2nd[iPair]->MomVtx(theVtx[iPair]->Position());
289     fixPh1st[iPair]->SetMom(newMom1st.X(), newMom1st.Y(), newMom1st.Z(), newMom1st.E());
290     fixPh2nd[iPair]->SetMom(newMom2nd.X(), newMom2nd.Y(), newMom2nd.Z(), newMom2nd.E());
291    
292     // check if both photons pass the CiC selection
293     // FIX-ME: Add other possibilities....
294     bool pass1 = false;
295     bool pass2 = false;
296 bendavid 1.8 switch( fPhSelType ){
297 fabstoec 1.1 case kNoPhSelection:
298     pass1 = ( fixPh1st[iPair]->Pt() > fLeadingPtMin );
299     pass2 = ( fixPh2nd[iPair]->Pt() > fTrailingPtMin );
300     break;
301     case kCiCPhSelection:
302 fabstoec 1.3
303    
304 fabstoec 1.2 pass1 = PhotonTools::PassCiCSelection(fixPh1st[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, fLeadingPtMin, fApplyEleVeto);
305     if(pass1) pass2 = PhotonTools::PassCiCSelection(fixPh2nd[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, fTrailingPtMin, fApplyEleVeto);
306 fabstoec 1.3
307 fabstoec 1.1 break;
308     case kMITPhSelection:
309 fabstoec 1.3 // FIX-ME: This is a place-holder.. MIT guys: Please work hard... ;)
310 fabstoec 1.1 pass1 = ( fixPh1st[iPair]->Pt() > fLeadingPtMin );
311     pass2 = ( fixPh2nd[iPair]->Pt() > fTrailingPtMin );
312     break;
313     default:
314     pass1 = true;
315     pass2 = true;
316     }
317 bendavid 1.8
318     //match to good electrons if requested
319     if (fInvertElectronVeto) {
320     pass1 &= !PhotonTools::PassElectronVeto(fixPh1st[iPair],fGoodElectrons);
321     pass2 &= !PhotonTools::PassElectronVeto(fixPh2nd[iPair],fGoodElectrons);
322     }
323 fabstoec 1.1 // finally, if both Photons pass the selections, add the pair to the 'passing Pairs)
324     if( pass1 && pass2 ) passPairs.push_back(iPair);
325     }
326    
327    
328     // ---------------------------------------------------------------
329     // ... we're almost done, stau focused...
330     // loop over all passing pairs and find the one with the highest sum Et
331     const Vertex* _theVtx = NULL;
332     Photon* phHard = NULL;
333     Photon* phSoft = NULL;
334    
335     PhotonTools::CiCBaseLineCats catPh1 = PhotonTools::kCiCNoCat;
336     PhotonTools::CiCBaseLineCats catPh2 = PhotonTools::kCiCNoCat;
337    
338     double maxSumEt = 0.;
339     for(unsigned int iPair=0; iPair<passPairs.size(); ++iPair){
340     double sumEt = fixPh1st[passPairs[iPair]]->Et();
341     sumEt += fixPh2nd[passPairs[iPair]]->Et();
342     if( sumEt > maxSumEt ) {
343     maxSumEt = sumEt;
344     phHard = fixPh1st[passPairs[iPair]];
345     phSoft = fixPh2nd[passPairs[iPair]];
346     catPh1 = cat1st[passPairs[iPair]];
347     catPh2 = cat2nd[passPairs[iPair]];
348     _theVtx = theVtx[iPair];
349     }
350     }
351    
352     // ---------------------------------------------------------------
353     // we have the Photons (*PARTY*)... compute some useful qunatities
354    
355 bendavid 1.8
356 fabstoec 1.1
357     if(phHard && phSoft) {
358     GoodPhotons->AddOwned(phHard);
359     GoodPhotons->AddOwned(phSoft);
360     }
361    
362    
363     // sort according to pt
364     GoodPhotons->Sort();
365    
366     // delete auxiliary photon collection...
367     delete preselPh;
368     delete[] theVtx;
369 bendavid 1.4
370 fabstoec 1.1 return;
371    
372     }
373    
374     //--------------------------------------------------------------------------------------------------
375     void PhotonPairSelector::SlaveBegin()
376     {
377     // Run startup code on the computer (slave) doing the actual analysis. Here,
378     // we just request the photon collection branch.
379    
380     ReqEventObject(fPhotonBranchName, fPhotons, fPhotonsFromBranch);
381     ReqEventObject(fTrackBranchName, fTracks, true);
382     ReqEventObject(fElectronName, fElectrons, true);
383 bendavid 1.8 ReqEventObject(fGoodElectronName, fGoodElectrons, fGoodElectronsFromBranch);
384 fabstoec 1.1 ReqEventObject(fPileUpDenName, fPileUpDen, true);
385     ReqEventObject(fPVName, fPV, fPVFromBranch);
386     ReqEventObject(fConversionName, fConversions,true);
387     ReqEventObject(fBeamspotName, fBeamspot, true);
388     ReqEventObject(fPFCandName, fPFCands, true);
389    
390     if (!fIsData) {
391     ReqBranch(fPileUpName, fPileUp);
392     ReqBranch(fMCParticleName, fMCParticles);
393     }
394    
395     if (fPhotonSelType.CompareTo("CiCSelection") == 0)
396     fPhSelType = kCiCPhSelection;
397     else if (fPhotonSelType.CompareTo("MITSelection") == 0)
398     fPhSelType = kMITPhSelection;
399     else
400     fPhSelType = kNoPhSelection;
401    
402     if (fVertexSelType.CompareTo("CiCSelection") == 0)
403     fVtxSelType = kCiCVtxSelection;
404     else if (fVertexSelType.CompareTo("MITSelection") == 0)
405     fVtxSelType = kMITVtxSelection;
406     else
407     fVtxSelType = kStdVtxSelection;
408    
409     }
410    
411     // ----------------------------------------------------------------------------------------
412     // some helpfer functions....
413 bendavid 1.8 void PhotonPairSelector::FindHiggsPtAndZ(Float_t& pt, Float_t& decayZ, Float_t& mass) {
414 fabstoec 1.1
415     pt = -999.;
416     decayZ = -999.;
417 bendavid 1.8 mass = -999.;
418 fabstoec 1.1
419     // loop over all GEN particles and look for status 1 photons
420     for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
421     const MCParticle* p = fMCParticles->At(i);
422 bendavid 1.8 if( p->Is(MCParticle::kH) || (!fApplyEleVeto && p->AbsPdgId()==23) ) {
423     pt=p->Pt();
424     decayZ = p->DecayVertex().Z();
425     mass = p->Mass();
426     break;
427     }
428 fabstoec 1.1 }
429    
430     return;
431     }
432    
433     // this routine looks for the idx of the run-range
434     Int_t PhotonPairSelector::FindRunRangeIdx(UInt_t run) {
435     Int_t runRange=-1;
436     for(UInt_t iRun = 0; iRun<fRunStart.size(); ++iRun) {
437     if( run >= fRunStart[iRun] && run <= fRunEnd[iRun]) {
438     runRange = (Int_t) iRun;
439     return runRange;
440     }
441     }
442     return runRange;
443     }
444    
445    
446     Double_t PhotonPairSelector::GetDataEnCorr(Int_t runRange, PhotonTools::CiCBaseLineCats cat) {
447     switch( cat ) {
448     case PhotonTools::kCiCCat1:
449     return fDataEnCorr_EB_hR9[runRange];
450     case PhotonTools::kCiCCat2:
451     return fDataEnCorr_EB_lR9[runRange];
452     case PhotonTools::kCiCCat3:
453     return fDataEnCorr_EE_hR9[runRange];
454     case PhotonTools::kCiCCat4:
455     return fDataEnCorr_EE_lR9[runRange];
456     default:
457     return 1.;
458     }
459     }
460    
461    
462     Double_t PhotonPairSelector::GetMCSmearFac(PhotonTools::CiCBaseLineCats cat) {
463     switch( cat ) {
464     case PhotonTools::kCiCCat1:
465     return fMCSmear_EB_hR9;
466     case PhotonTools::kCiCCat2:
467     return fMCSmear_EB_lR9;
468     case PhotonTools::kCiCCat3:
469     return fMCSmear_EE_hR9;
470     case PhotonTools::kCiCCat4:
471     return fMCSmear_EE_lR9;
472     default:
473     return 1.;
474     }
475     }
476    
477     Float_t PhotonPairSelector::GetEventCat(PhotonTools::CiCBaseLineCats cat1, PhotonTools::CiCBaseLineCats cat2) {
478    
479     bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
480     bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
481    
482     bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
483     bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
484    
485     if( ph1IsEB && ph2IsEB )
486     return ( ph1IsHR9 && ph2IsHR9 ? 0. : 1.);
487    
488     return ( ph1IsHR9 && ph2IsHR9 ? 2. : 3.);
489     }
490 bendavid 1.4