ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonCiCMod.cc
Revision: 1.8
Committed: Wed Jul 6 13:59:40 2011 UTC (13 years, 10 months ago) by fabstoec
Content type: text/plain
Branch: MAIN
Changes since 1.7: +73 -46 lines
Log Message: 
update for EPS reload

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