ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonCiCMod.cc
Revision: 1.9
Committed: Fri Jul 8 09:24:45 2011 UTC (13 years, 10 months ago) by fabstoec
Content type: text/plain
Branch: MAIN
Changes since 1.8: +59 -15 lines
Log Message: 
*** empty log message ***

File Contents

# User Rev Content
1 fabstoec 1.1 #include "MitPhysics/Mods/interface/PhotonCiCMod.h"
2     #include "MitAna/DataTree/interface/PhotonCol.h"
3     #include "MitPhysics/Init/interface/ModNames.h"
4     #include "MitPhysics/Utils/interface/IsolationTools.h"
5     #include "MitPhysics/Utils/interface/PhotonTools.h"
6 fabstoec 1.4 #include <TNtuple.h>
7 fabstoec 1.6 #include <TRandom3.h>
8 fabstoec 1.1
9     using namespace mithep;
10    
11     ClassImp(mithep::PhotonCiCMod)
12    
13     //--------------------------------------------------------------------------------------------------
14     PhotonCiCMod::PhotonCiCMod(const char *name, const char *title) :
15     BaseMod(name,title),
16     fPhotonBranchName (Names::gkPhotonBrn),
17     fGoodPhotonsName (ModNames::gkGoodPhotonsName),
18     fTrackBranchName (Names::gkTrackBrn),
19     fPileUpDenName (Names::gkPileupEnergyDensityBrn),
20 fabstoec 1.9
21 fabstoec 1.1 fElectronName ("Electrons"),
22 fabstoec 1.9 fPhotonPtMin (20.0),
23     fAbsEtaMax (999.99),
24     fApplySpikeRemoval (kFALSE),
25    
26 fabstoec 1.1 fPhotons(0),
27     fTracks(0),
28     fPileUpDen(0),
29     fElectrons(0),
30    
31 fabstoec 1.4 fPVName(Names::gkPVBeamSpotBrn),
32 fabstoec 1.1 fPV(0),
33 fabstoec 1.4 fPVFromBranch(kTRUE),
34    
35     fConversions(0),
36     fConversionName(Names::gkMvfConversionBrn),
37    
38 fabstoec 1.5 fBeamspot(0),
39    
40 fabstoec 1.6 // May10 ReReco
41 fabstoec 1.8 fDataEnCorr_EB_hR9(0),
42     fDataEnCorr_EB_lR9(0),
43     fDataEnCorr_EE_hR9(0),
44     fDataEnCorr_EE_lR9(0),
45 fabstoec 1.6
46 fabstoec 1.8 fRunStart(0),
47     fRunEnd(0),
48 fabstoec 1.6
49     fMCSmear_EB_hR9(0.0089),
50     fMCSmear_EB_lR9(0.0199),
51     fMCSmear_EE_hR9(0.0409),
52     fMCSmear_EE_lR9(0.0246),
53    
54     fIsData(false),
55    
56     rng(new TRandom3()),
57 fabstoec 1.9
58 fabstoec 1.6 fMCParticleName(Names::gkMCPartBrn),
59     fMCParticles(0),
60     fPileUpName ("PileupInfo"),
61     fPileUp (0)
62 fabstoec 1.1
63     {
64     // Constructor.
65     }
66    
67 fabstoec 1.6 PhotonCiCMod::~PhotonCiCMod(){
68     if(rng) delete rng;
69     }
70    
71 fabstoec 1.1 //--------------------------------------------------------------------------------------------------
72     void PhotonCiCMod::Process()
73     {
74     // Process entries of the tree.
75    
76     LoadEventObject(fPhotonBranchName, fPhotons);
77 bendavid 1.2
78     PhotonOArr *GoodPhotons = new PhotonOArr;
79     GoodPhotons->SetName(fGoodPhotonsName);
80 fabstoec 1.4 GoodPhotons->SetOwner(kTRUE);
81    
82 fabstoec 1.9 Double_t _tRho = 0.;
83 fabstoec 1.4 LoadEventObject(fTrackBranchName, fTracks);
84     LoadEventObject(fPileUpDenName, fPileUpDen);
85     LoadEventObject(fElectronName, fElectrons);
86     LoadEventObject(fPVName, fPV);
87     LoadEventObject(fConversionName, fConversions);
88     LoadBranch("BeamSpot");
89 fabstoec 1.6
90     if( !fIsData ) {
91     LoadBranch(fMCParticleName);
92     LoadBranch(fPileUpName);
93     }
94    
95     Float_t numPU = -1.;
96     if( !fIsData )
97     numPU = (Float_t) fPileUp->At(0)->GetPU_NumInteractions();
98    
99 fabstoec 1.4 if(fPileUpDen->GetEntries() > 0)
100 fabstoec 1.9 _tRho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
101 fabstoec 1.4
102     bool doVtxSelection = true;
103 fabstoec 1.8 bool doMCSmear = true;
104 fabstoec 1.4
105     const EventHeader* evtHead = this->GetEventHeader();
106    
107     unsigned int evtNum = evtHead->EvtNum();
108     Float_t _evtNum1 = (Float_t) ( (int) (evtNum/10000.) );
109     Float_t _evtNum2 = (Float_t) ( (int) (evtNum % 10000) );
110 fabstoec 1.6
111 fabstoec 1.7 //double evtNumTest = (int) ( ( (double) _evtNum1 )*10000. + (double) _evtNum2 );
112 fabstoec 1.4
113 fabstoec 1.8 UInt_t runNumber = evtHead->RunNum();
114     Float_t _runNum = (Float_t) runNumber;
115 fabstoec 1.4 Float_t _lumiSec = (Float_t) evtHead->LumiSec();
116    
117 fabstoec 1.6
118 fabstoec 1.7 //unsigned int numVertices = fPV->GetEntries();
119 fabstoec 1.4
120     const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));
121 fabstoec 1.3
122     PhotonOArr* preselPh = new PhotonOArr;
123 fabstoec 1.1
124 fabstoec 1.3 // 1. we do the pre-selection; but keep the non-passing photons in a secont container...
125     for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
126     const Photon *ph = fPhotons->At(i);
127 fabstoec 1.4
128     if(ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) continue;
129     if(ph->Et() < 20.) continue;
130     if(ph->HadOverEm() > 0.15) continue;
131 fabstoec 1.3 if(ph->AbsEta() < 1.5) {
132 fabstoec 1.4 if(ph->CoviEtaiEta() > 0.013) continue;
133 fabstoec 1.3 } else {
134 fabstoec 1.4 if(ph->CoviEtaiEta() > 0.03) continue;
135 fabstoec 1.3 }
136 fabstoec 1.4
137     Bool_t passSpikeRemovalFilter = kTRUE;
138    
139     if (ph->SCluster() && ph->SCluster()->Seed()) {
140     if(ph->SCluster()->Seed()->Energy() > 5.0 &&
141     ph->SCluster()->Seed()->EMax() / ph->SCluster()->Seed()->E3x3() > 0.95 )
142     passSpikeRemovalFilter = kFALSE;
143     }
144     if (fApplySpikeRemoval && !passSpikeRemovalFilter) continue;
145    
146 fabstoec 1.3 preselPh->Add(ph);
147     }
148 fabstoec 1.9
149 fabstoec 1.4 // sort both by pt... again ;)
150 fabstoec 1.3 preselPh->Sort();
151 fabstoec 1.9
152 fabstoec 1.4 unsigned int numPairs = 0;
153     if( preselPh->GetEntries() > 0)
154     numPairs = (preselPh->GetEntries()-1)*preselPh->GetEntries()/2;
155    
156     // create all possible pairs of pre-selected photons
157    
158     std::vector<unsigned int> idxFst;
159     std::vector<unsigned int> idxSec;
160     std::vector<bool> pairPasses;
161    
162     if(numPairs > 0) {
163     for(unsigned int iFst = 0; iFst <preselPh->GetEntries() - 1; ++iFst) {
164     for(unsigned int iSec = iFst + 1; iSec <preselPh->GetEntries(); ++iSec) {
165     idxFst.push_back(iFst);
166     idxSec.push_back(iSec);
167     pairPasses.push_back(true);
168     }
169     }
170     }
171    
172     // array to store the index of 'chosen Vtx' for each pair
173     const Vertex** theVtx = new const Vertex*[numPairs];
174 fabstoec 1.8 UInt_t* theVtxIdx = new UInt_t[numPairs];
175    
176 fabstoec 1.4 // arays to store the Vtx 'fixed' photons
177     Photon** fixPhFst = new Photon*[numPairs];
178     Photon** fixPhSec = new Photon*[numPairs];
179    
180     // sotr pair-indices for pairs passing the selection
181     std::vector<unsigned int> passPairs;
182     passPairs.resize(0);
183    
184     unsigned int theChosenVtx = 0;
185    
186     float kinPh1[20];
187     float kinPh2[20];
188    
189     for(int i=0; i<10; ++i) {
190     kinPh1[i] =-99.;
191     kinPh2[i] =-99.;
192     }
193    
194     float ptBefore1 = -99.;
195     float ptBefore2 = -99.;
196    
197     bool print = false;
198 fabstoec 1.8 if(evtNum == 17031 && false) {
199 fabstoec 1.4 std::cout<<" ------------------------------------------- "<<std::endl;
200     std::cout<<" printing info for event #"<<evtNum<<std::endl;
201     print = true;
202     }
203    
204     for(unsigned int iPair = 0; iPair < numPairs; ++iPair) {
205    
206     // copy the photons for manipulation
207     fixPhFst[iPair] = new Photon(*preselPh->At(idxFst[iPair]));
208     fixPhSec[iPair] = new Photon(*preselPh->At(idxSec[iPair]));
209 fabstoec 1.6
210     // if this is Data, scale the energy, if MC smear...
211     FourVectorM scMomFst = fixPhFst[iPair]->Mom();
212     FourVectorM scMomSec = fixPhSec[iPair]->Mom();
213     double scaleFac1 = 1.;
214     double scaleFac2 = 1.;
215     if (fIsData) {
216 fabstoec 1.8 if( fRunStart.size() > 0) {
217     // find run in rnage
218     Int_t runRange=-1;
219     for(UInt_t iRun = 0; iRun<fRunStart.size(); ++iRun) {
220     if( runNumber >= fRunStart[iRun] && runNumber <= fRunEnd[iRun]) {
221     runRange = (Int_t) iRun;
222     break;
223     }
224     }
225     if(runRange > -1) {
226     if(fixPhFst[iPair]->IsEB())
227     if(fixPhFst[iPair]->R9() > 0.94)
228     scaleFac1 += fDataEnCorr_EB_hR9[runRange];
229     else
230     scaleFac1 += fDataEnCorr_EB_lR9[runRange];
231     else
232     if(fixPhFst[iPair]->R9() > 0.94)
233     scaleFac1 += fDataEnCorr_EE_hR9[runRange];
234     else
235     scaleFac1 += fDataEnCorr_EE_lR9[runRange];
236     if(fixPhSec[iPair]->IsEB())
237     if(fixPhSec[iPair]->R9() > 0.94)
238     scaleFac2 += fDataEnCorr_EB_hR9[runRange];
239     else
240     scaleFac2 += fDataEnCorr_EB_lR9[runRange];
241     else
242     if(fixPhSec[iPair]->R9() > 0.94)
243     scaleFac2 += fDataEnCorr_EE_hR9[runRange];
244     else
245     scaleFac2 += fDataEnCorr_EE_lR9[runRange];
246     }
247     }
248 fabstoec 1.6 } else {
249     // get the smearing for MC photons..
250     UInt_t seedBase = (UInt_t) evtNum + (UInt_t) _runNum + (UInt_t) _lumiSec;
251     UInt_t seed1 = seedBase + (UInt_t) fixPhFst[iPair]->E() + (UInt_t) (TMath::Abs(10.*fixPhFst[iPair]->SCluster()->Eta()));
252     UInt_t seed2 = seedBase + (UInt_t) fixPhSec[iPair]->E() + (UInt_t) (TMath::Abs(10.*fixPhSec[iPair]->SCluster()->Eta()));
253    
254     double width1 = 0.;
255     double width2 = 0.;
256     if(fixPhFst[iPair]->IsEB())
257     if(fixPhFst[iPair]->R9() > 0.94)
258     width1 = fMCSmear_EB_hR9;
259     else
260     width1 = fMCSmear_EB_lR9;
261     else
262     if(fixPhFst[iPair]->R9() > 0.94)
263     width1 = fMCSmear_EE_hR9;
264     else
265     width1 = fMCSmear_EE_lR9;
266 fabstoec 1.9
267    
268 fabstoec 1.6 if(fixPhSec[iPair]->IsEB())
269     if(fixPhSec[iPair]->R9() > 0.94)
270     width2 = fMCSmear_EB_hR9;
271     else
272     width2 = fMCSmear_EB_lR9;
273     else
274     if(fixPhSec[iPair]->R9() > 0.94)
275     width2 = fMCSmear_EE_hR9;
276     else
277     width2 = fMCSmear_EE_lR9;
278    
279     if(doMCSmear) {
280     rng->SetSeed(seed1);
281     scaleFac1 = rng->Gaus(1.,width1);
282     rng->SetSeed(seed2);
283     scaleFac2 = rng->Gaus(1.,width2);
284     }
285     }
286 fabstoec 1.9
287 fabstoec 1.6 if(iPair==0) {
288     ptBefore1 = fixPhFst[iPair]->Pt();
289     ptBefore2 = fixPhSec[iPair]->Pt();
290     }
291 fabstoec 1.9
292 fabstoec 1.6 if(print && false) {
293     std::cout<<" Photon Pair #"<<iPair+1<<std::endl;
294     std::cout<<" Ph1 px = "<<fixPhFst[iPair]->Mom().X()<<std::endl;
295     std::cout<<" py = "<<fixPhFst[iPair]->Mom().Y()<<std::endl;
296     std::cout<<" pz = "<<fixPhFst[iPair]->Mom().Z()<<std::endl;
297     std::cout<<" E = "<<fixPhFst[iPair]->Mom().E()<<std::endl;
298     std::cout<<" M = "<<fixPhFst[iPair]->Mom().M()<<std::endl;
299     std::cout<<" Ph2 px = "<<fixPhSec[iPair]->Mom().X()<<std::endl;
300     std::cout<<" py = "<<fixPhSec[iPair]->Mom().Y()<<std::endl;
301     std::cout<<" pz = "<<fixPhSec[iPair]->Mom().Z()<<std::endl;
302     std::cout<<" E = "<<fixPhSec[iPair]->Mom().E()<<std::endl;
303     std::cout<<" M = "<<fixPhSec[iPair]->Mom().M()<<std::endl;
304     }
305    
306     fixPhFst[iPair]->SetMom(scaleFac1*scMomFst.X(), scaleFac1*scMomFst.Y(), scaleFac1*scMomFst.Z(), scaleFac1*scMomFst.E());
307     fixPhSec[iPair]->SetMom(scaleFac2*scMomSec.X(), scaleFac2*scMomSec.Y(), scaleFac2*scMomSec.Z(), scaleFac2*scMomSec.E());
308    
309     if(print && false) {
310     std::cout<<" SF 1 = "<<scaleFac1<<std::endl;
311     std::cout<<" SF 2 = "<<scaleFac2<<std::endl;
312    
313     std::cout<<" Photon Pair #"<<iPair+1<<std::endl;
314     std::cout<<" Ph1 px = "<<fixPhFst[iPair]->Mom().X()<<std::endl;
315     std::cout<<" py = "<<fixPhFst[iPair]->Mom().Y()<<std::endl;
316     std::cout<<" pz = "<<fixPhFst[iPair]->Mom().Z()<<std::endl;
317     std::cout<<" E = "<<fixPhFst[iPair]->Mom().E()<<std::endl;
318     std::cout<<" M = "<<fixPhFst[iPair]->Mom().M()<<std::endl;
319     std::cout<<" Ph2 px = "<<fixPhSec[iPair]->Mom().X()<<std::endl;
320     std::cout<<" py = "<<fixPhSec[iPair]->Mom().Y()<<std::endl;
321     std::cout<<" pz = "<<fixPhSec[iPair]->Mom().Z()<<std::endl;
322     std::cout<<" E = "<<fixPhSec[iPair]->Mom().E()<<std::endl;
323     std::cout<<" M = "<<fixPhSec[iPair]->Mom().M()<<std::endl;
324     }
325 fabstoec 1.9
326 fabstoec 1.6 // store the vertex for this pair
327 fabstoec 1.4 if(doVtxSelection) {
328     unsigned int iVtx = findBestVertex(fixPhFst[iPair],fixPhSec[iPair],bsp, print);
329 fabstoec 1.8 theVtx[iPair] = fPV->At(iVtx);
330     theVtxIdx[iPair] = iVtx;
331 fabstoec 1.4 if(iPair == 0) theChosenVtx = iVtx;
332     } else
333     theVtx[iPair] = fPV->At(0);
334 fabstoec 1.6
335 fabstoec 1.4 // fix the kinematics for both events
336     FourVectorM newMomFst = fixPhFst[iPair]->MomVtx(theVtx[iPair]->Position());
337     FourVectorM newMomSec = fixPhSec[iPair]->MomVtx(theVtx[iPair]->Position());
338     fixPhFst[iPair]->SetMom(newMomFst.X(), newMomFst.Y(), newMomFst.Z(), newMomFst.E());
339     fixPhSec[iPair]->SetMom(newMomSec.X(), newMomSec.Y(), newMomSec.Z(), newMomSec.E());
340 fabstoec 1.6
341    
342     if(print && false) {
343     std::cout<<" Vtx = "<<theChosenVtx<<std::endl;
344     std::cout<<" Photon Pair #"<<iPair+1<<std::endl;
345     std::cout<<" Ph1 px = "<<fixPhFst[iPair]->Mom().X()<<std::endl;
346     std::cout<<" py = "<<fixPhFst[iPair]->Mom().Y()<<std::endl;
347     std::cout<<" pz = "<<fixPhFst[iPair]->Mom().Z()<<std::endl;
348     std::cout<<" E = "<<fixPhFst[iPair]->Mom().E()<<std::endl;
349     std::cout<<" M = "<<fixPhFst[iPair]->Mom().M()<<std::endl;
350     std::cout<<" Ph2 px = "<<fixPhSec[iPair]->Mom().X()<<std::endl;
351     std::cout<<" py = "<<fixPhSec[iPair]->Mom().Y()<<std::endl;
352     std::cout<<" pz = "<<fixPhSec[iPair]->Mom().Z()<<std::endl;
353     std::cout<<" E = "<<fixPhSec[iPair]->Mom().E()<<std::endl;
354     std::cout<<" M = "<<fixPhSec[iPair]->Mom().M()<<std::endl;
355     }
356    
357    
358 fabstoec 1.4
359     if(iPair != 0) {
360     // check if both photons pass the CiC selection
361     bool pass1 = PhotonTools::PassCiCSelection(fixPhFst[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 40., false);
362     bool pass2 = PhotonTools::PassCiCSelection(fixPhSec[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 30., false);
363     if( pass1 && pass2 )
364     passPairs.push_back(iPair);
365     } else {
366     bool pass1 = PhotonTools::PassCiCSelection(fixPhFst[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 40., false, kinPh1);
367     bool pass2 = PhotonTools::PassCiCSelection(fixPhSec[iPair], theVtx[iPair], fTracks, fElectrons, fPV, _tRho, 30., false, kinPh2);
368     if( pass1 && pass2 )
369     passPairs.push_back(iPair);
370     }
371     }
372    
373     // loop over all passing pairs and find the one with the highest sum Et
374     Photon* phHard = NULL;
375     Photon* phSoft = NULL;
376 fabstoec 1.8
377     const Vertex* _theVtx = NULL;
378    
379 fabstoec 1.4 double maxSumEt = 0.;
380     for(unsigned int iPair=0; iPair<passPairs.size(); ++iPair){
381     double sumEt = fixPhFst[passPairs[iPair]]->Et();
382     sumEt += fixPhSec[passPairs[iPair]]->Et();
383     if( sumEt > maxSumEt ) {
384     maxSumEt = sumEt;
385     phHard = fixPhFst[passPairs[iPair]];
386     phSoft = fixPhSec[passPairs[iPair]];
387 fabstoec 1.8 _theVtx = theVtx[iPair];
388     theChosenVtx = theVtxIdx[iPair];
389 fabstoec 1.4 }
390     }
391    
392 fabstoec 1.8 Float_t _theVtxZ = -999.;
393 fabstoec 1.9 Float_t catPh1 = 0.;
394     Float_t catPh2 = 0.;
395     Float_t catEvt = -1.;
396    
397     bool ph1IsLowR9 = false;
398     bool ph2IsLowR9 = false;
399    
400 fabstoec 1.4 if(phHard && phSoft) {
401     GoodPhotons->AddOwned(phHard);
402     GoodPhotons->AddOwned(phSoft);
403 fabstoec 1.8 if(_theVtx)
404     _theVtxZ=_theVtx->Position().Z();
405 fabstoec 1.9
406     catPh1=1.;
407     catPh2=1.;
408     catEvt=0.;
409    
410     if(phHard->SCluster()->AbsEta()>1.5)
411     catPh1=3;
412     if(phHard->R9() < 0.94) {
413     catPh1=catPh1+1;
414     ph1IsLowR9 = true;
415     }
416    
417     if(phSoft->SCluster()->AbsEta()>1.5)
418     catPh2=3;
419    
420     if(phSoft->R9() < 0.94) {
421     catPh2=catPh2+1;
422     ph2IsLowR9 = true;
423     }
424    
425     if(catPh1 > 2.5 || catPh2 > 2.5)
426     catEvt=2.;
427     if(ph1IsLowR9 || ph2IsLowR9)
428     catEvt=catEvt+1.;
429 fabstoec 1.4 }
430 fabstoec 1.9
431 fabstoec 1.4 // sort according to pt
432     GoodPhotons->Sort();
433    
434     // add to event for other modules to use
435     AddObjThisEvt(GoodPhotons);
436    
437     delete preselPh;
438    
439     bool doFill = (phHard && phSoft);
440     Float_t _mass = ( doFill ? (phHard->Mom()+phSoft->Mom()).M() : -100.);
441     Float_t _ptgg = ( doFill ? (phHard->Mom()+phSoft->Mom()).Pt() : -100.);
442 fabstoec 1.9 if(_ptgg < 40. && doFill) catEvt = catEvt+4.;
443    
444 fabstoec 1.8 Float_t _pth = -100.;
445     Float_t _decayZ = -100.;
446     if( !fIsData ) findHiggsPtAndZ(_pth, _decayZ);
447 fabstoec 1.4
448 fabstoec 1.5 Float_t fillEvent[] = { _tRho,
449 fabstoec 1.6 _pth,
450 fabstoec 1.8 _decayZ,
451     _theVtxZ,
452 fabstoec 1.6 numPU,
453 fabstoec 1.5 _mass,
454 fabstoec 1.4 _ptgg,
455     _evtNum1,
456     _evtNum2,
457     _runNum,
458     _lumiSec,
459     (float) theChosenVtx,
460     (float) numPairs,
461 fabstoec 1.9 catPh1,
462     catPh2,
463     catEvt,
464 fabstoec 1.4 kinPh1[0],
465     kinPh1[1],
466     kinPh1[2],
467     kinPh1[3],
468     kinPh1[4],
469     kinPh1[5],
470     kinPh1[6],
471     kinPh1[7],
472     kinPh1[8],
473     kinPh1[9],
474     kinPh1[10],
475     kinPh1[11],
476     kinPh1[12],
477     kinPh1[13],
478     kinPh1[14],
479     kinPh1[15],
480     kinPh1[16],
481     kinPh1[17],
482     kinPh1[18],
483     kinPh1[19],
484     kinPh2[0],
485     kinPh2[1],
486     kinPh2[2],
487     kinPh2[3],
488     kinPh2[4],
489     kinPh2[5],
490     kinPh2[6],
491     kinPh2[7],
492     kinPh2[8],
493     kinPh2[9],
494     kinPh2[10],
495     kinPh2[11],
496     kinPh2[12],
497     kinPh2[13],
498     kinPh2[14],
499     kinPh2[15],
500     kinPh2[16],
501     kinPh2[17],
502     kinPh2[18],
503     kinPh2[19],
504     ptBefore1,
505     ptBefore2
506     };
507    
508    
509 fabstoec 1.9
510     if(_mass > 0.) {
511     //std::cout<<catPh1<<" "<<catPh2<<" "<<catEvt<<" "<<_mass<<std::endl;
512     hCiCTuple->Fill(fillEvent);
513     }
514    
515 fabstoec 1.8 delete[] theVtx;
516     delete[] theVtxIdx;
517    
518 fabstoec 1.4 return;
519    
520     }
521    
522     //--------------------------------------------------------------------------------------------------
523     void PhotonCiCMod::SlaveBegin()
524     {
525     // Run startup code on the computer (slave) doing the actual analysis. Here,
526     // we just request the photon collection branch.
527    
528     ReqEventObject(fPhotonBranchName, fPhotons, kTRUE);
529     ReqEventObject(fTrackBranchName, fTracks, kTRUE);
530     ReqEventObject(fElectronName, fElectrons, kTRUE);
531     ReqEventObject(fPileUpDenName, fPileUpDen, kTRUE);
532     ReqEventObject(fPVName, fPV, fPVFromBranch);
533     ReqEventObject(fConversionName, fConversions,kTRUE);
534     ReqBranch("BeamSpot",fBeamspot);
535 fabstoec 1.6
536     if (!fIsData) {
537     ReqBranch(fPileUpName, fPileUp);
538     ReqBranch(Names::gkMCPartBrn,fMCParticles);
539     }
540 fabstoec 1.4
541 fabstoec 1.9 hCiCTuple = new TNtuple("hCiCTuple","hCiCTuple","rho:higgspt:higgsZ:vtxZ:numPU:mass:ptgg:evtnum1:evtnum2:runnum:lumisec:ivtx:npairs:ph1Cat:ph2Cat:evtCat:ph1Iso1:ph1Iso2:ph1Iso3:ph1Cov:ph1HoE:ph1R9:ph1DR:ph1Pt:ph1Eta:ph1Phi:ph1Eiso3:ph1Eiso4:ph1Hiso4:ph1TisoA:ph1TisoW:ph1Tiso:ph1Et:ph1E:ph1Pass:ph1CatDebug:ph2Iso1:ph2Iso2:ph2Iso3:ph2Cov:ph2HoE:ph2R9:ph2DR:ph2Pt:ph2Eta:ph2Phi:ph2Eiso3:ph2Eiso4:ph2Hiso4:ph2TisoA:ph2TisoW:ph2Tiso:ph2Et:ph2E:ph2Pass:ph2CatDebug:ph1UPt:ph2UPt");
542 fabstoec 1.3
543 fabstoec 1.4 AddOutput(hCiCTuple);
544    
545     }
546    
547     // return the index of the bext vertex
548     unsigned int PhotonCiCMod::findBestVertex(Photon* ph1, Photon* ph2, const BaseVertex *bsp, bool print) {
549 fabstoec 1.3
550     // loop over all vertices and assigne the ranks
551     int* ptbal_rank = new int[fPV->GetEntries()];
552     int* ptasym_rank = new int[fPV->GetEntries()];
553     int* total_rank = new int[fPV->GetEntries()];
554     double* ptbal = new double[fPV->GetEntries()];
555     double* ptasym = new double[fPV->GetEntries()];
556 fabstoec 1.4
557     unsigned int numVertices = fPV->GetEntries();
558    
559     double ptgg = 0.; // stored for later in the conversion
560 fabstoec 1.3
561 fabstoec 1.4 if(print) {
562     std::cout<<" --------------------------------- "<<std::endl;
563     std::cout<<" looking for Vtx for photon pair "<<std::endl;
564     std::cout<<" pt 1 = "<<ph1->Et()<<std::endl;
565     std::cout<<" pt 2 = "<<ph2->Et()<<std::endl;
566     std::cout<<" among #"<<numVertices<<" Vertices."<<std::endl;
567     }
568    
569    
570     for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
571     if(print)
572     std::cout<<std::endl<<" Vertex #"<<iVtx<<std::endl;
573    
574     const Vertex* tVtx = fPV->At(iVtx);
575 fabstoec 1.3 ptbal [iVtx] = 0.0;
576     ptasym[iVtx] = 0.0;
577     ptbal_rank [iVtx] = 1;
578     ptasym_rank[iVtx] = 1;
579 fabstoec 1.4
580     // compute the photon momenta with respect to this Vtx
581     FourVectorM newMomFst = ph1->MomVtx(tVtx->Position());
582     FourVectorM newMomSec = ph2->MomVtx(tVtx->Position());
583    
584     FourVectorM higgsMom = newMomFst+newMomSec;
585    
586     double ph1Eta = newMomFst.Eta();
587     double ph2Eta = newMomSec.Eta();
588    
589     double ph1Phi = newMomFst.Phi();
590     double ph2Phi = newMomSec.Phi();
591    
592     if(print && iVtx == 0 && false ) {
593     std::cout<<" new momenta Et1 = "<<newMomFst.Et()<<std::endl;
594     std::cout<<" Eta = "<<newMomFst.Eta()<<std::endl;
595     std::cout<<" Phi = "<<newMomFst.Phi()<<std::endl;
596     std::cout<<" Px = "<<newMomFst.Px()<<std::endl;
597     std::cout<<" Py = "<<newMomFst.Py()<<std::endl;
598     std::cout<<" new momenta Et2 = "<<newMomSec.Et()<<std::endl;
599     std::cout<<" Eta = "<<newMomSec.Eta()<<std::endl;
600     std::cout<<" Phi = "<<newMomSec.Phi()<<std::endl;
601     std::cout<<" Px = "<<newMomSec.Px()<<std::endl;
602     std::cout<<" Py = "<<newMomSec.Py()<<std::endl;
603     }
604    
605     FourVectorM totTrkMom;
606     for(unsigned int iTrk = 0; iTrk < tVtx->NTracks(); ++iTrk) {
607     const Track* tTrk = tVtx->Trk(iTrk);
608     //if(tTrk->Pt()<1.) continue;
609     double tEta = tTrk->Eta();
610     double tPhi = tTrk->Phi();
611     double dEta1 = TMath::Abs(tEta-ph1Eta);
612     double dEta2 = TMath::Abs(tEta-ph2Eta);
613     double dPhi1 = TMath::Abs(tPhi-ph1Phi);
614     double dPhi2 = TMath::Abs(tPhi-ph2Phi);
615     if(dPhi1 > M_PI) dPhi1 = 2*M_PI - dPhi1;
616     if(dPhi2 > M_PI) dPhi2 = 2*M_PI - dPhi2;
617    
618     double dR1 = TMath::Sqrt(dEta1*dEta1+dPhi1*dPhi1);
619     double dR2 = TMath::Sqrt(dEta2*dEta2+dPhi2*dPhi2);
620    
621 fabstoec 1.6 if( ( iVtx == 0 || iVtx == 1 ) && print && false) {
622 fabstoec 1.4 std::cout<<" Track #"<<iTrk<<std::endl;
623     std::cout<<" pt = "<<tTrk->Pt()<<std::endl;
624     std::cout<<" eta = "<<tTrk->Eta()<<std::endl;
625     std::cout<<" phi = "<<tTrk->Phi()<<std::endl;
626     std::cout<<" px = "<<tTrk->Px()<<std::endl;
627     std::cout<<" py = "<<tTrk->Py()<<std::endl;
628     std::cout<<" dR1 = "<<dR1<<std::endl;
629     std::cout<<" dR2 = "<<dR2<<std::endl;
630     }
631    
632     if(dR1 < 0.05 || dR2 < 0.05) continue;
633    
634     if(iTrk == 0) totTrkMom = tTrk->Mom4(0);
635     else totTrkMom = totTrkMom + tTrk->Mom4(0);
636    
637     }
638    
639     if(iVtx ==0 && print && false) {
640     std::cout<<" Trk passes cuts: "<<std::endl;
641     std::cout<<" px tot = "<<totTrkMom.Px()<<std::endl;
642     std::cout<<" py tot = "<<totTrkMom.Py()<<std::endl;
643     }
644    
645     double ptvtx = totTrkMom.Pt();
646     if(iVtx ==0 && print && false)
647     std::cout<<" Total TkPt = "<<ptvtx<<std::endl;
648     double pthiggs = higgsMom.Pt();
649     if(iVtx ==0 && print && false)
650     std::cout<<" Total H Pt = "<<pthiggs<<std::endl;
651     if(iVtx == 0) ptgg = pthiggs;
652     ptbal [iVtx] = (totTrkMom.Px()*(newMomFst.Px()+newMomSec.Px()));
653     ptbal [iVtx] += (totTrkMom.Py()*(newMomFst.Py()+newMomSec.Py()));
654     ptbal [iVtx] = -ptbal[iVtx]/pthiggs;
655     ptasym[iVtx] = (ptvtx - pthiggs)/(ptvtx + pthiggs);
656    
657     if(iVtx ==0 && print && false) {
658     std::cout<<" Results: ptbal = "<<ptbal [iVtx]<<std::endl;
659     std::cout<<" ptasym = "<<ptasym[iVtx]<<std::endl;
660     }
661    
662     for(unsigned int cVtx =0; cVtx < iVtx; ++cVtx) {
663 fabstoec 1.3 if(ptbal [iVtx] > ptbal [cVtx])
664 fabstoec 1.4 ptbal_rank[cVtx]++;
665 fabstoec 1.3 else
666 fabstoec 1.4 ptbal_rank[iVtx]++;
667 fabstoec 1.3 if(ptasym [iVtx] > ptasym [cVtx])
668 fabstoec 1.4 ptasym_rank[cVtx]++;
669 fabstoec 1.3 else
670 fabstoec 1.4 ptasym_rank[iVtx]++;
671 fabstoec 1.3 }
672     }
673 fabstoec 1.4
674 fabstoec 1.1
675 fabstoec 1.3 // compute the total rank
676 fabstoec 1.4 for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
677     if(print) {
678 fabstoec 1.6 std::cout<<" Vertex #"<<iVtx<<" has rank PTB "<<ptbal_rank[iVtx]<<" ("<<ptbal[iVtx]<<")"<<std::endl;
679 fabstoec 1.4 std::cout<<" Vertex #"<<iVtx<<" has rank PTSYM "<<ptasym_rank[iVtx]<<" ("<<ptasym[iVtx]<<")"<<std::endl;
680     }
681     ptasym_rank [iVtx] = ptbal_rank [iVtx]*ptasym_rank [iVtx]*(iVtx+1);
682 fabstoec 1.3 total_rank [iVtx] = 0;
683 fabstoec 1.4 for(unsigned int cVtx =0; cVtx < iVtx; ++cVtx) {
684     if(ptasym_rank [iVtx] > ptasym_rank [cVtx])
685 fabstoec 1.3 total_rank[iVtx]++;
686 fabstoec 1.4 else if(ptasym_rank [iVtx] == ptasym_rank [cVtx]) {
687     if(ptbal_rank [iVtx] > ptbal_rank [cVtx])
688     total_rank[iVtx]++;
689     else
690     total_rank[cVtx]++;
691     }
692 fabstoec 1.3 else
693     total_rank[cVtx]++;
694     }
695 bendavid 1.2 }
696 fabstoec 1.1
697 fabstoec 1.4 if(print) std::cout<<std::endl;
698    
699     unsigned int bestIdx = 0;
700     for(unsigned int iVtx = 0; iVtx < numVertices; ++iVtx) {
701     if(total_rank[iVtx] == 0) bestIdx = iVtx;
702     if(print) {
703     std::cout<<" Vertex #"<<iVtx<<" has rank "<<total_rank[iVtx]<<std::endl;
704     }
705     }
706 fabstoec 1.3
707 fabstoec 1.6 delete[] ptbal_rank ;
708     delete[] ptasym_rank ;
709     delete[] ptbal ;
710     delete[] ptasym ;
711 fabstoec 1.1
712 fabstoec 1.4 //return bestIdx;
713    
714     // check if there's a conversion among the pre-selected photons
715 fabstoec 1.6 const DecayParticle* conv1 = PhotonTools::MatchedCiCConversion(ph1, fConversions, 0.1, 0.1, 0.1, print);
716     const DecayParticle* conv2 = PhotonTools::MatchedCiCConversion(ph2, fConversions, 0.1, 0.1, 0.1, print);
717    
718     if(print) {
719     if (conv1) {
720     std::cout<<" Photon 1 has has conversion with P = "<<conv1->Prob()<<std::endl;
721     std::cout<<" Rho = "<<conv1->Position().Rho()<<std::endl;
722     std::cout<<" Z = "<<conv1->Position().Z()<<std::endl;
723     }
724     if (conv2) {
725     std::cout<<" Photon 2 has has conversion with P = "<<conv2->Prob()<<std::endl;
726     std::cout<<" Rho = "<<conv2->Position().Rho()<<std::endl;
727     std::cout<<" Z = "<<conv2->Position().Z()<<std::endl;
728     }
729     }
730 fabstoec 1.5
731 fabstoec 1.6 if( conv1 && ( conv1->Prob() < 0.0005) ) conv1 = NULL;
732     if( conv2 && ( conv2->Prob() < 0.0005) ) conv2 = NULL;
733    
734 fabstoec 1.4 double zconv = 0.;
735     double dzconv = 0.;
736 fabstoec 1.6
737 fabstoec 1.4 if(conv1 || conv2) {
738 fabstoec 1.6 if( !conv2 ){
739 fabstoec 1.5 //const mithep::ThreeVector caloPos1(ph1->SCluster()->Point());
740     const mithep::ThreeVector caloPos1(ph1->CaloPos());
741 fabstoec 1.4 zconv = conv1->Z0EcalVtx(bsp->Position(), caloPos1);
742 fabstoec 1.5 if( ph1->IsEB() ) {
743 fabstoec 1.4 double rho = conv1->Position().Rho();
744     if ( rho < 15. ) dzconv = 0.06;
745     else if( rho < 60. ) dzconv = 0.67;
746     else dzconv = 2.04;
747     } else {
748     double z = conv1->Position().Z();
749 fabstoec 1.6 if ( TMath::Abs(z) < 50. ) dzconv = 0.18;
750     else if( TMath::Abs(z) < 100.) dzconv = 0.61;
751 fabstoec 1.4 else dzconv = 0.99;
752     }
753     } else if( !conv1 ) {
754 fabstoec 1.5 //const mithep::ThreeVector caloPos2(ph2->SCluster()->Point());
755     const mithep::ThreeVector caloPos2(ph2->CaloPos());
756 fabstoec 1.4 zconv = conv2->Z0EcalVtx(bsp->Position(), caloPos2);
757 fabstoec 1.5 if( ph2->IsEB() ) {
758 fabstoec 1.4 double rho = conv2->Position().Rho();
759     if ( rho < 15. ) dzconv = 0.06;
760     else if( rho < 60. ) dzconv = 0.67;
761     else dzconv = 2.04;
762     } else {
763     double z = conv2->Position().Z();
764 fabstoec 1.6 if ( TMath::Abs(z) < 50. ) dzconv = 0.18;
765     else if( TMath::Abs(z) < 100.) dzconv = 0.61;
766 fabstoec 1.4 else dzconv = 0.99;
767     }
768     } else {
769 fabstoec 1.5 //const mithep::ThreeVector caloPos1(ph1->SCluster()->Point());
770     const mithep::ThreeVector caloPos1(ph1->CaloPos());
771 fabstoec 1.4 double z1 = conv1->Z0EcalVtx(bsp->Position(), caloPos1);
772     double dz1 = 0.;
773 fabstoec 1.5 if( ph1->IsEB() ) {
774 fabstoec 1.4 double rho = conv1->Position().Rho();
775     if ( rho < 15. ) dz1 = 0.06;
776     else if( rho < 60. ) dz1 = 0.67;
777     else dz1 = 2.04;
778     } else {
779     double z = conv1->Position().Z();
780 fabstoec 1.6 if ( TMath::Abs(z) < 50. ) dz1 = 0.18;
781     else if( TMath::Abs(z) < 100.) dz1 = 0.61;
782 fabstoec 1.4 else dz1 = 0.99;
783     }
784 fabstoec 1.5 //const mithep::ThreeVector caloPos2(ph2->SCluster()->Point());
785     const mithep::ThreeVector caloPos2(ph2->CaloPos());
786 fabstoec 1.4 double z2 = conv2->Z0EcalVtx(bsp->Position(), caloPos2);
787     double dz2 = 0.;
788 fabstoec 1.5 if( ph2->IsEB() ) {
789 fabstoec 1.4 double rho = conv2->Position().Rho();
790     if ( rho < 15. ) dz2 = 0.06;
791     else if( rho < 60. ) dz2 = 0.67;
792     else dz2 = 2.04;
793     } else {
794     double z = conv2->Position().Z();
795 fabstoec 1.6 if ( TMath::Abs(z) < 50. ) dz2 = 0.18;
796     else if( TMath::Abs(z) < 100.) dz2 = 0.61;
797 fabstoec 1.4 else dz2 = 0.99;
798     }
799 fabstoec 1.6
800     if(print) {
801     std::cout<<" z1 = "<<z1<<std::endl;
802     std::cout<<" dz1 = "<<dz1<<std::endl;
803     std::cout<<" z2 = "<<z2<<std::endl;
804     std::cout<<" dz2 = "<<dz2<<std::endl;
805     }
806    
807     zconv = ( 1./(1./dz1/dz1 + 1./dz2/dz2 )*(z1/dz1/dz1 + z2/dz2/dz2) ) ; // weighted average
808 fabstoec 1.4 dzconv = TMath::Sqrt( 1./(1./dz1/dz1 + 1./dz2/dz2)) ;
809 fabstoec 1.1 }
810 fabstoec 1.5
811 fabstoec 1.6 if(print) {
812     std::cout<<" Conversion Z = "<<zconv<<std::endl;
813     std::cout<<" dZ = "<<dzconv<<std::endl;
814     }
815    
816    
817 fabstoec 1.8 if(true) {
818 fabstoec 1.6
819     // loop over all ranked Vertices and choose the closest to the Conversion one
820     int maxVertices = ( ptgg > 30 ? 3 : 5);
821     double minDz = -1.;
822    
823     if(print) std::cout<<std::endl<<" looping over vertices... "<<ptgg<<" "<<maxVertices<<std::endl;
824    
825     for(unsigned int iVtx =0; iVtx < numVertices; ++iVtx) {
826    
827     if(print) std::cout<<" "<<iVtx<<" has rank "<<total_rank[iVtx]<<std::endl;
828    
829     if(total_rank[iVtx] < maxVertices) {
830     const Vertex* tVtx = fPV->At(iVtx);
831     double tDz = TMath::Abs(zconv - tVtx->Z());
832     if(print) std::cout<<" is considered with tDz = "<<tDz<<std::endl;
833     if( (minDz < 0. || tDz < minDz) && ( tDz < dzconv ) ) {
834     minDz = tDz;
835     bestIdx = iVtx;
836     if(print) std::cout<<" and is the best now."<<std::endl;
837     }
838     }
839     }
840     } else {
841     unsigned int bestIdxTmp = bestIdx;
842    
843     // loop over all ranked Vertices and choose the closest to the Conversion one
844     double minDz = -1.;
845     int maxVertices = ( ptgg > 30 ? 3 : 5);
846    
847     if(print) std::cout<<std::endl<<" looping over vertices... "<<ptgg<<" "<<maxVertices<<std::endl;
848    
849     for(unsigned int iVtx =0; iVtx < numVertices; ++iVtx) {
850    
851     if(print) std::cout<<" "<<iVtx<<" has rank "<<total_rank[iVtx]<<std::endl;
852    
853 fabstoec 1.5 const Vertex* tVtx = fPV->At(iVtx);
854     double tDz = TMath::Abs(zconv - tVtx->Z());
855 fabstoec 1.6 if(print) std::cout<<" is considered with tDz = "<<tDz<<std::endl;
856     if( (minDz < 0. || tDz < minDz) && ( tDz < dzconv ) ) {
857 fabstoec 1.5 minDz = tDz;
858 fabstoec 1.6 bestIdxTmp = iVtx;
859     if(print) std::cout<<" and is the best now."<<std::endl;
860 fabstoec 1.5 }
861 fabstoec 1.4 }
862 fabstoec 1.6
863     // check if best Vtx is among higest ranked ones
864     if(total_rank[bestIdxTmp] < maxVertices)
865     bestIdx = bestIdxTmp;
866     }
867 fabstoec 1.1 }
868    
869 fabstoec 1.6 delete[] total_rank ;
870 fabstoec 1.4 return bestIdx;
871 fabstoec 1.1 }
872 fabstoec 1.3
873 fabstoec 1.8 void PhotonCiCMod::findHiggsPtAndZ(Float_t& pt, Float_t& decayZ) {
874    
875     pt = -999.;
876     decayZ = -999.;
877 fabstoec 1.6
878     // loop over all GEN particles and look for status 1 photons
879     for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
880     const MCParticle* p = fMCParticles->At(i);
881     if( !(p->Is(MCParticle::kH)) ) continue;
882 fabstoec 1.8 pt=p->Pt();
883     decayZ = p->DecayVertex().Z();
884     break;
885 fabstoec 1.6 }
886 fabstoec 1.8
887     return;
888 fabstoec 1.6 }