ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonCiCMod.cc
Revision: 1.7
Committed: Mon Jul 4 13:48:50 2011 UTC (13 years, 10 months ago) by fabstoec
Content type: text/plain
Branch: MAIN
Changes since 1.6: +2 -2 lines
Log Message: 
fixed warnings

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