ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonPairSelector.cc
Revision: 1.2
Committed: Fri Jul 15 17:24:37 2011 UTC (13 years, 9 months ago) by fabstoec
Content type: text/plain
Branch: MAIN
Changes since 1.1: +29 -8 lines
Log Message:
added/improved stuff for CiC Selection

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     double width1 = GetMCSmearFac(cat1st[iPair]);
249     double width2 = GetMCSmearFac(cat2nd[iPair]);
250     PhotonTools::SmearPhoton(fixPh1st[iPair], rng, width1, seed1);
251     PhotonTools::SmearPhoton(fixPh2nd[iPair], rng, width2, seed2);
252     }
253     }
254    
255     // store the vertex for this pair
256     switch( fVtxSelType ){
257     case kStdVtxSelection:
258     theVtx[iPair] = fPV->At(0);
259     break;
260     case kCiCVtxSelection:
261     theVtx[iPair] = VertexTools::findVtxBasicRanking(fixPh1st[iPair],fixPh2nd[iPair], bsp, fPV, fConversions);
262     break;
263     case kMITVtxSelection:
264     // need PFCandidate Collection
265     theVtx[iPair] = VertexTools::BestVtx(fPFCands, fPV, bsp, mithep::FourVector((fixPh1st[iPair]->Mom()+fixPh2nd[iPair]->Mom())));
266     break;
267     default:
268     theVtx[iPair] = fPV->At(0);
269 fabstoec 1.2
270 fabstoec 1.1 }
271 fabstoec 1.2
272    
273 fabstoec 1.1 // fix the kinematics for both events
274     FourVectorM newMom1st = fixPh1st[iPair]->MomVtx(theVtx[iPair]->Position());
275     FourVectorM newMom2nd = fixPh2nd[iPair]->MomVtx(theVtx[iPair]->Position());
276     fixPh1st[iPair]->SetMom(newMom1st.X(), newMom1st.Y(), newMom1st.Z(), newMom1st.E());
277     fixPh2nd[iPair]->SetMom(newMom2nd.X(), newMom2nd.Y(), newMom2nd.Z(), newMom2nd.E());
278    
279     // check if both photons pass the CiC selection
280     // FIX-ME: Add other possibilities....
281     bool pass1 = false;
282     bool pass2 = false;
283     switch( fVtxSelType ){
284     case kNoPhSelection:
285     pass1 = ( fixPh1st[iPair]->Pt() > fLeadingPtMin );
286     pass2 = ( fixPh2nd[iPair]->Pt() > fTrailingPtMin );
287     break;
288     case kCiCPhSelection:
289 fabstoec 1.2 pass1 = PhotonTools::PassCiCSelection(fixPh1st[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, fLeadingPtMin, fApplyEleVeto);
290     if(pass1) pass2 = PhotonTools::PassCiCSelection(fixPh2nd[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, fTrailingPtMin, fApplyEleVeto);
291 fabstoec 1.1 break;
292     case kMITPhSelection:
293     // FIX-ME: This is a place-holder.. MIT guys: Please worj hard... ;)
294     pass1 = ( fixPh1st[iPair]->Pt() > fLeadingPtMin );
295     pass2 = ( fixPh2nd[iPair]->Pt() > fTrailingPtMin );
296     break;
297     default:
298     pass1 = true;
299     pass2 = true;
300     }
301     // finally, if both Photons pass the selections, add the pair to the 'passing Pairs)
302     if( pass1 && pass2 ) passPairs.push_back(iPair);
303     }
304    
305    
306     // ---------------------------------------------------------------
307     // ... we're almost done, stau focused...
308     // loop over all passing pairs and find the one with the highest sum Et
309     const Vertex* _theVtx = NULL;
310     Photon* phHard = NULL;
311     Photon* phSoft = NULL;
312    
313     PhotonTools::CiCBaseLineCats catPh1 = PhotonTools::kCiCNoCat;
314     PhotonTools::CiCBaseLineCats catPh2 = PhotonTools::kCiCNoCat;
315    
316     double maxSumEt = 0.;
317     for(unsigned int iPair=0; iPair<passPairs.size(); ++iPair){
318     double sumEt = fixPh1st[passPairs[iPair]]->Et();
319     sumEt += fixPh2nd[passPairs[iPair]]->Et();
320     if( sumEt > maxSumEt ) {
321     maxSumEt = sumEt;
322     phHard = fixPh1st[passPairs[iPair]];
323     phSoft = fixPh2nd[passPairs[iPair]];
324     catPh1 = cat1st[passPairs[iPair]];
325     catPh2 = cat2nd[passPairs[iPair]];
326     _theVtx = theVtx[iPair];
327     }
328     }
329    
330     // ---------------------------------------------------------------
331     // we have the Photons (*PARTY*)... compute some useful qunatities
332    
333     Float_t _theVtxZ = -999.;
334    
335     if(phHard && phSoft) {
336     GoodPhotons->AddOwned(phHard);
337     GoodPhotons->AddOwned(phSoft);
338     }
339    
340     if(_theVtx) _theVtxZ=_theVtx->Position().Z();
341    
342     // sort according to pt
343     GoodPhotons->Sort();
344    
345     // add to event for other modules to use
346     AddObjThisEvt(GoodPhotons);
347    
348     // delete auxiliary photon collection...
349     delete preselPh;
350     delete[] theVtx;
351    
352     // Fill the useful information (mass and di-photon pt)
353     bool doFill = ( phHard && phSoft );
354     Float_t _mass = ( doFill ? (phHard->Mom()+phSoft->Mom()).M() : -100.);
355     Float_t _ptgg = ( doFill ? (phHard->Mom()+phSoft->Mom()).Pt() : -100.);
356    
357     // in case of a MC event, try to find Higgs and Higgs decay Z poisition
358     Float_t _pth = -100.;
359     Float_t _decayZ = -100.;
360     if( !fIsData ) FindHiggsPtAndZ(_pth, _decayZ);
361    
362     Float_t evtCat = -1.;
363     if( doFill ) {
364     evtCat = GetEventCat(catPh1, catPh2);
365     if(_ptgg < 40.) evtCat += 4.;
366     }
367    
368     Float_t fillEvent[] = { _tRho,
369     _pth,
370     _decayZ,
371     _theVtxZ,
372     _numPU,
373     _mass,
374     _ptgg,
375     _evtNum1,
376     _evtNum2,
377     _runNum,
378     _lumiSec,
379     (Float_t) catPh1,
380     (Float_t) catPh2,
381     evtCat
382     };
383    
384     // to keep the Tree slim, only add in case we haqve found a passing pair...
385     if(_mass > 0.)
386     hCiCTuple->Fill(fillEvent);
387    
388     return;
389    
390     }
391    
392     //--------------------------------------------------------------------------------------------------
393     void PhotonPairSelector::SlaveBegin()
394     {
395     // Run startup code on the computer (slave) doing the actual analysis. Here,
396     // we just request the photon collection branch.
397    
398     ReqEventObject(fPhotonBranchName, fPhotons, fPhotonsFromBranch);
399     ReqEventObject(fTrackBranchName, fTracks, true);
400     ReqEventObject(fElectronName, fElectrons, true);
401     ReqEventObject(fPileUpDenName, fPileUpDen, true);
402     ReqEventObject(fPVName, fPV, fPVFromBranch);
403     ReqEventObject(fConversionName, fConversions,true);
404     ReqEventObject(fBeamspotName, fBeamspot, true);
405     ReqEventObject(fPFCandName, fPFCands, true);
406    
407     if (!fIsData) {
408     ReqBranch(fPileUpName, fPileUp);
409     ReqBranch(fMCParticleName, fMCParticles);
410     }
411    
412     if (fPhotonSelType.CompareTo("CiCSelection") == 0)
413     fPhSelType = kCiCPhSelection;
414     else if (fPhotonSelType.CompareTo("MITSelection") == 0)
415     fPhSelType = kMITPhSelection;
416     else
417     fPhSelType = kNoPhSelection;
418    
419     if (fVertexSelType.CompareTo("CiCSelection") == 0)
420     fVtxSelType = kCiCVtxSelection;
421     else if (fVertexSelType.CompareTo("MITSelection") == 0)
422     fVtxSelType = kMITVtxSelection;
423     else
424     fVtxSelType = kStdVtxSelection;
425    
426 fabstoec 1.2 hCiCTuple = new TNtuple(fTupleName.Data(),fTupleName.Data(),"rho:higgspt:higgsZ:vtxZ:numPU:mass:ptgg:evtnum1:evtnum2:runnum:lumisec:ph1Cat:ph2Cat:evtCat");
427 fabstoec 1.1
428     AddOutput(hCiCTuple);
429    
430     }
431    
432     // ----------------------------------------------------------------------------------------
433     // some helpfer functions....
434     void PhotonPairSelector::FindHiggsPtAndZ(Float_t& pt, Float_t& decayZ) {
435    
436     pt = -999.;
437     decayZ = -999.;
438    
439     // loop over all GEN particles and look for status 1 photons
440     for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
441     const MCParticle* p = fMCParticles->At(i);
442     if( !(p->Is(MCParticle::kH)) ) continue;
443     pt=p->Pt();
444     decayZ = p->DecayVertex().Z();
445     break;
446     }
447    
448     return;
449     }
450    
451     // this routine looks for the idx of the run-range
452     Int_t PhotonPairSelector::FindRunRangeIdx(UInt_t run) {
453     Int_t runRange=-1;
454     for(UInt_t iRun = 0; iRun<fRunStart.size(); ++iRun) {
455     if( run >= fRunStart[iRun] && run <= fRunEnd[iRun]) {
456     runRange = (Int_t) iRun;
457     return runRange;
458     }
459     }
460     return runRange;
461     }
462    
463    
464     Double_t PhotonPairSelector::GetDataEnCorr(Int_t runRange, PhotonTools::CiCBaseLineCats cat) {
465     switch( cat ) {
466     case PhotonTools::kCiCCat1:
467     return fDataEnCorr_EB_hR9[runRange];
468     case PhotonTools::kCiCCat2:
469     return fDataEnCorr_EB_lR9[runRange];
470     case PhotonTools::kCiCCat3:
471     return fDataEnCorr_EE_hR9[runRange];
472     case PhotonTools::kCiCCat4:
473     return fDataEnCorr_EE_lR9[runRange];
474     default:
475     return 1.;
476     }
477     }
478    
479    
480     Double_t PhotonPairSelector::GetMCSmearFac(PhotonTools::CiCBaseLineCats cat) {
481     switch( cat ) {
482     case PhotonTools::kCiCCat1:
483     return fMCSmear_EB_hR9;
484     case PhotonTools::kCiCCat2:
485     return fMCSmear_EB_lR9;
486     case PhotonTools::kCiCCat3:
487     return fMCSmear_EE_hR9;
488     case PhotonTools::kCiCCat4:
489     return fMCSmear_EE_lR9;
490     default:
491     return 1.;
492     }
493     }
494    
495     Float_t PhotonPairSelector::GetEventCat(PhotonTools::CiCBaseLineCats cat1, PhotonTools::CiCBaseLineCats cat2) {
496    
497     bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
498     bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
499    
500     bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
501     bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
502    
503     if( ph1IsEB && ph2IsEB )
504     return ( ph1IsHR9 && ph2IsHR9 ? 0. : 1.);
505    
506     return ( ph1IsHR9 && ph2IsHR9 ? 2. : 3.);
507     }