ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonPairSelector.cc
Revision: 1.1
Committed: Fri Jul 8 17:54:27 2011 UTC (13 years, 10 months ago) by fabstoec
Content type: text/plain
Branch: MAIN
Log Message:
*** empty log message ***

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