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

# Content
1 #include "MitPhysics/Mods/interface/PhotonPairSelector.h"
2 #include "MitAna/DataTree/interface/PhotonCol.h"
3 #include "MitAna/DataTree/interface/PFCandidateCol.h"
4 #include "MitAna/DataTree/interface/StableData.h"
5 #include "MitAna/DataTree/interface/StableParticle.h"
6 #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 #include "MitPhysics/Utils/interface/PhotonFix.h"
11 #include "TDataMember.h"
12 #include <TNtuple.h>
13 #include <TRandom3.h>
14 #include <TSystem.h>
15
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 fGoodElectronName (Names::gkElectronBrn),
29 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 fGoodElectronsFromBranch (kTRUE),
59
60 // ----------------------------------------
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 fDoVtxSelection (true),
94 fApplyEleVeto (true),
95 fInvertElectronVeto(kFALSE)
96 {
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
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 LoadEventObject(fElectronName, fElectrons);
122 LoadEventObject(fGoodElectronName, fGoodElectrons);
123 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
131
132
133
134
135 // ------------------------------------------------------------
136 // load event based information
137 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 Float_t _tRho = -99.;
142 if( fPileUpDen->GetEntries() > 0 )
143 _tRho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
144
145
146
147 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 //Float_t _evtNum1 = (Float_t) ( (int) (evtNum/10000.) );
154 //Float_t _evtNum2 = (Float_t) ( (int) (evtNum % 10000) );
155 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 }
179
180 if (preselPh->GetEntries()<2) return;
181
182 // 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 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
256 double width1 = GetMCSmearFac(cat1st[iPair]);
257 double width2 = GetMCSmearFac(cat2nd[iPair]);
258
259 PhotonTools::SmearPhoton(fixPh1st[iPair], rng, width1, seed1);
260 PhotonTools::SmearPhoton(fixPh2nd[iPair], rng, width2, seed2);
261
262 }
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
280 }
281
282 //set PV ref in photons
283 fixPh1st[iPair]->SetPV(theVtx[iPair]);
284 fixPh2nd[iPair]->SetPV(theVtx[iPair]);
285
286 // 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 switch( fPhSelType ){
297 case kNoPhSelection:
298 pass1 = ( fixPh1st[iPair]->Pt() > fLeadingPtMin );
299 pass2 = ( fixPh2nd[iPair]->Pt() > fTrailingPtMin );
300 break;
301 case kCiCPhSelection:
302
303
304 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
307 break;
308 case kMITPhSelection:
309 // FIX-ME: This is a place-holder.. MIT guys: Please work hard... ;)
310 pass1 = ( fixPh1st[iPair]->Pt() > fLeadingPtMin );
311 pass2 = ( fixPh2nd[iPair]->Pt() > fTrailingPtMin );
312 break;
313 default:
314 pass1 = true;
315 pass2 = true;
316 }
317
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 // 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
356
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
370 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 ReqEventObject(fGoodElectronName, fGoodElectrons, fGoodElectronsFromBranch);
384 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 void PhotonPairSelector::FindHiggsPtAndZ(Float_t& pt, Float_t& decayZ, Float_t& mass) {
414
415 pt = -999.;
416 decayZ = -999.;
417 mass = -999.;
418
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 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 }
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