ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonPairSelector.cc
Revision: 1.3
Committed: Fri Jul 15 19:42:36 2011 UTC (13 years, 9 months ago) by fabstoec
Content type: text/plain
Branch: MAIN
Changes since 1.2: +7 -2 lines
Log Message:
bug fixes in CiCSelection

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