ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/IsolationTools.cc
Revision: 1.26
Committed: Sat Dec 31 23:16:13 2011 UTC (13 years, 4 months ago) by sixie
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_025e, Mit_025d
Changes since 1.25: +13 -3 lines
Log Message: 
add option to select particular pf candidate types

File Contents

# Content
1 // $Id: IsolationTools.cc,v 1.25 2011/12/21 13:05:47 sixie Exp $
2
3 #include "MitPhysics/Utils/interface/IsolationTools.h"
4 #include "MitPhysics/Utils/interface/PhotonTools.h"
5 #include "MitCommon/MathTools/interface/MathUtils.h"
6
7 ClassImp(mithep::IsolationTools)
8
9 using namespace mithep;
10
11 //--------------------------------------------------------------------------------------------------
12 Double_t IsolationTools::TrackIsolation(const Track *p, Double_t extRadius, Double_t intRadius,
13 Double_t ptLow, Double_t maxVtxZDist,
14 const Collection<Track> *tracks)
15 {
16 //Computes the Track Isolation: Summed Transverse Momentum of all tracks inside an
17 //annulus around the electron seed track.
18
19 Double_t ptSum =0.;
20 for (UInt_t i=0; i<tracks->GetEntries();i++) {
21 Double_t tmpPt = tracks->At(i)->Pt();
22 Double_t deltaZ = fabs(p->Z0() - tracks->At(i)->Z0());
23
24 //ignore the track if it is below the pt threshold
25 if (tmpPt < ptLow)
26 continue;
27 //ingore the track if it is too far away in Z
28 if (deltaZ > maxVtxZDist)
29 continue;
30
31 Double_t dr = MathUtils::DeltaR(p->Phi(),p->Eta(),tracks->At(i)->Phi(), tracks->At(i)->Eta());
32 //add the track pt if it is inside the annulus
33 if ( dr < extRadius &&
34 dr >= intRadius ) {
35 ptSum += tmpPt;
36 }
37 }
38 return ptSum;
39 }
40
41 //--------------------------------------------------------------------------------------------------
42 Double_t IsolationTools::EcalIsolation(const SuperCluster *sc, Double_t coneSize, Double_t etLow,
43 const Collection<BasicCluster> *basicClusters)
44 {
45 //Computes the Ecal Isolation: Summed Transverse Energy of all Basic Clusters inside a
46 //cone around the electron, excluding those that are inside the electron super cluster.
47
48 Double_t ecalIsol=0.;
49 const BasicCluster *basicCluster= 0;
50 for (UInt_t i=0; i<basicClusters->GetEntries();i++) {
51 basicCluster = basicClusters->At(i);
52 Double_t basicClusterEnergy = basicCluster->Energy();
53 Double_t basicClusterEta = basicCluster->Eta();
54 Double_t basicClusterEt = basicClusterEnergy*sin(2*atan(exp(basicClusterEta)));
55
56 if (basicClusterEt > etLow) {
57 bool inSuperCluster = false;
58
59 // loop over the basic clusters of the supercluster
60 // to make sure that the basic cluster is not inside
61 // the super cluster. We exclude those.
62 for (UInt_t j=0; j<sc->ClusterSize(); j++) {
63 const BasicCluster *tempBasicClusterInSuperCluster = sc->Cluster(j);
64 if (tempBasicClusterInSuperCluster == basicCluster) {
65 inSuperCluster = true;
66 }
67 }
68
69 if (!inSuperCluster) {
70 Double_t dr = MathUtils::DeltaR(sc->Phi(), sc->Eta(),
71 basicCluster->Phi(),basicCluster->Eta());
72 if(dr < coneSize) {
73 ecalIsol += basicClusterEt;
74 }
75 }
76 }
77 }
78 return ecalIsol;
79 }
80
81 //--------------------------------------------------------------------------------------------------
82 Double_t IsolationTools::CaloTowerHadIsolation(const ThreeVector *p, Double_t extRadius,
83 Double_t intRadius, Double_t etLow,
84 const Collection<CaloTower> *caloTowers)
85 {
86 //Computes the CaloTower Had Et Isolation: Summed Hadronic Transverse Energy of all Calo Towers
87 //inside an annulus around the electron super cluster position.
88
89 Double_t sumEt = 0;
90 for (UInt_t i=0; i<caloTowers->GetEntries();i++) {
91 Double_t caloTowerEt = caloTowers->At(i)->HadEt();
92 Double_t dr = MathUtils::DeltaR(caloTowers->At(i)->Phi(), caloTowers->At(i)->Eta(),
93 p->Phi(), p->Eta());
94 if (dr < extRadius && dr > intRadius && caloTowerEt > etLow) {
95 sumEt += caloTowerEt;
96 }
97 }
98 return sumEt;
99 }
100
101 //--------------------------------------------------------------------------------------------------
102 Double_t IsolationTools::CaloTowerEmIsolation(const ThreeVector *p, Double_t extRadius,
103 Double_t intRadius, Double_t etLow,
104 const Collection<CaloTower> *caloTowers)
105 {
106 //Computes the CaloTower Em Et Isolation: Summed Hadronic Transverse Energy of all Calo Towers
107 //inside an annulus around the electron super cluster position.
108
109 Double_t sumEt = 0;
110 for (UInt_t i=0; i<caloTowers->GetEntries();i++) {
111 Double_t caloTowerEt = caloTowers->At(i)->EmEt();
112 Double_t dr = MathUtils::DeltaR(caloTowers->At(i)->Phi(), caloTowers->At(i)->Eta(),
113 p->Phi(), p->Eta());
114 if (dr < extRadius && dr > intRadius && caloTowerEt > etLow) {
115 sumEt += caloTowerEt;
116 }
117 }
118 return sumEt;
119 }
120
121 //--------------------------------------------------------------------------------------------------
122 Double_t IsolationTools::PFMuonIsolation(const Muon *p, const PFCandidateCol *PFCands,
123 const Vertex *vertex, Double_t delta_z, Double_t ptMin,
124 Double_t extRadius, Double_t intRadiusGamma, Double_t intRadius)
125 {
126 //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
127 //annulus around the particle seed track.
128
129 Double_t zLepton = 0.0;
130 if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
131
132 Double_t ptSum =0.;
133 for (UInt_t i=0; i<PFCands->GetEntries();i++) {
134 const PFCandidate *pf = PFCands->At(i);
135
136 // exclude muon
137 if(pf->TrackerTrk() && p->TrackerTrk() &&
138 pf->TrackerTrk() == p->TrackerTrk()) continue;
139
140 Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
141
142 // pt cut applied to neutrals
143 if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
144
145 // ignore the pf candidate if it is too far away in Z
146 if(pf->HasTrk()) {
147 Double_t deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
148 if (deltaZ >= delta_z) continue;
149 }
150
151 // inner cone veto for gammas
152 if (pf->PFType() == PFCandidate::eGamma && dr < intRadiusGamma) continue;
153
154 // inner cone veto for tracks
155 if (dr < intRadius) continue;
156
157 // add the pf pt if it is inside the extRadius
158 if (dr < extRadius) ptSum += pf->Pt();
159
160 }
161 return ptSum;
162 }
163
164 //--------------------------------------------------------------------------------------------------
165 Double_t IsolationTools::PFMuonIsolation(const Muon *p, const PFCandidateCol *PFCands,
166 const MuonCol *goodMuons, const ElectronCol *goodElectrons,
167 const Vertex *vertex, Double_t delta_z, Double_t ptMin,
168 Double_t extRadius, Double_t intRadiusGamma, Double_t intRadius)
169 {
170 //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
171 //annulus around the particle seed track.
172
173 Double_t zLepton = 0.0;
174 if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
175
176 Double_t ptSum =0.;
177 for (UInt_t i=0; i<PFCands->GetEntries();i++) {
178 const PFCandidate *pf = PFCands->At(i);
179
180 // exclude muon
181 if(pf->TrackerTrk() && p->TrackerTrk() &&
182 pf->TrackerTrk() == p->TrackerTrk()) continue;
183
184 Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
185
186 // pt cut applied to neutrals
187 if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
188
189 // ignore the pf candidate if it is too far away in Z
190 if(pf->HasTrk()) {
191 Double_t deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
192 if (deltaZ >= delta_z) continue;
193 }
194
195 // inner cone veto for gammas
196 if (pf->PFType() == PFCandidate::eGamma && dr < intRadiusGamma) continue;
197
198 // inner cone veto for tracks
199 if (dr < intRadius) continue;
200
201 // add the pf pt if it is inside the extRadius and outside the intRadius
202 if (dr < extRadius ) {
203 Bool_t isLepton = kFALSE;
204 if(goodMuons){
205 for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
206 const Muon *m = goodMuons->At(nl);
207 if(pf->TrackerTrk() && m->TrackerTrk() &&
208 pf->TrackerTrk() == m->TrackerTrk()) {
209 isLepton = kTRUE;
210 break;
211 }
212 }
213 }
214 if(goodElectrons && isLepton == kFALSE){
215 for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
216 const Electron *e = goodElectrons->At(nl);
217 if(pf->TrackerTrk() && e->TrackerTrk() &&
218 pf->TrackerTrk() == e->TrackerTrk()) {
219 isLepton = kTRUE;
220 break;
221 }
222 if(pf->GsfTrk() && e->GsfTrk() &&
223 pf->GsfTrk() == e->GsfTrk()) {
224 isLepton = kTRUE;
225 break;
226 }
227 }
228 }
229 if (isLepton == kTRUE) continue;
230 ptSum += pf->Pt();
231 }
232 }
233 return ptSum;
234 }
235
236 //--------------------------------------------------------------------------------------------------
237 Double_t IsolationTools::PFElectronIsolation(const Electron *p, const PFCandidateCol *PFCands,
238 const Vertex *vertex, Double_t delta_z, Double_t ptMin,
239 Double_t extRadius, Double_t intRadius, Int_t PFCandidateType)
240 {
241
242 //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
243 //annulus around the particle seed track.
244
245 Double_t zLepton = 0.0;
246 if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
247
248 Double_t ptSum =0.;
249 for (UInt_t i=0; i<PFCands->GetEntries();i++) {
250 const PFCandidate *pf = PFCands->At(i);
251
252 //select only specified PFCandidate types
253 if (PFCandidateType >= 0) {
254 if (pf->PFType() != PFCandidateType) continue;
255 }
256
257 // pt cut applied to neutrals
258 if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
259
260 if(pf->TrackerTrk() && p->TrackerTrk() &&
261 pf->TrackerTrk() == p->TrackerTrk()) continue;
262
263 if(pf->GsfTrk() && p->GsfTrk() &&
264 pf->GsfTrk() == p->GsfTrk()) continue;
265
266 // ignore the pf candidate if it is too far away in Z
267 if(pf->BestTrk()) {
268 Double_t deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
269 if (deltaZ >= delta_z) continue;
270 }
271
272
273 Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
274 // add the pf pt if it is inside the extRadius and outside the intRadius
275 if ( dr < extRadius &&
276 dr >= intRadius ) {
277
278 //EtaStrip Veto for Gamma
279 if (pf->PFType() == PFCandidate::eGamma && fabs(p->Eta() - pf->Eta()) < 0.025) continue;
280
281 //InnerCone (One Tower = dR < 0.07) Veto for non-gamma neutrals
282 if (!pf->HasTrk() && pf->PFType() == PFCandidate::eNeutralHadron
283 && MathUtils::DeltaR(p->Mom(), pf->Mom()) < 0.07 ) continue;
284
285 ptSum += pf->Pt();
286
287 }
288 }
289 return ptSum;
290 }
291
292
293 //--------------------------------------------------------------------------------------------------
294 Double_t IsolationTools::PFElectronIsolation(const Electron *p, const PFCandidateCol *PFCands,
295 const MuonCol *goodMuons, const ElectronCol *goodElectrons,
296 const Vertex *vertex, Double_t delta_z, Double_t ptMin,
297 Double_t extRadius, Double_t intRadius, Int_t PFCandidateType)
298 {
299
300 //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
301 //annulus around the particle seed track.
302
303 Double_t zLepton = 0.0;
304 if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
305
306 Double_t ptSum =0.;
307 for (UInt_t i=0; i<PFCands->GetEntries();i++) {
308 const PFCandidate *pf = PFCands->At(i);
309
310 //select only specified PFCandidate types
311 if (PFCandidateType >= 0) {
312 if (pf->PFType() != PFCandidateType) continue;
313 }
314
315 // pt cut applied to neutrals
316 if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
317
318 if(pf->TrackerTrk() && p->TrackerTrk() &&
319 pf->TrackerTrk() == p->TrackerTrk()) continue;
320
321 if(pf->GsfTrk() && p->GsfTrk() &&
322 pf->GsfTrk() == p->GsfTrk()) continue;
323
324 // ignore the pf candidate if it is too far away in Z
325 if(pf->BestTrk()) {
326 Double_t deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
327 if (deltaZ >= delta_z) continue;
328 }
329
330 Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
331 // add the pf pt if it is inside the extRadius and outside the intRadius
332 if ( dr < extRadius &&
333 dr >= intRadius ) {
334
335 //EtaStrip Veto for Gamma
336 if (pf->PFType() == PFCandidate::eGamma && fabs(p->Eta() - pf->Eta()) < 0.025) continue;
337
338 //InnerCone (One Tower = dR < 0.07) Veto for non-gamma neutrals
339 if (!pf->HasTrk() && pf->PFType() == PFCandidate::eNeutralHadron
340 && MathUtils::DeltaR(p->Mom(), pf->Mom()) < 0.07 ) continue;
341
342 Bool_t isLepton = kFALSE;
343 if(goodMuons){
344 for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
345 const Muon *m = goodMuons->At(nl);
346 if(pf->TrackerTrk() && m->TrackerTrk() &&
347 pf->TrackerTrk() == m->TrackerTrk()) {
348 isLepton = kTRUE;
349 break;
350 }
351 }
352 }
353 if(goodElectrons && isLepton == kFALSE){
354 for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
355 const Electron *e = goodElectrons->At(nl);
356 if(pf->TrackerTrk() && e->TrackerTrk() &&
357 pf->TrackerTrk() == e->TrackerTrk()) {
358 isLepton = kTRUE;
359 break;
360 }
361 if(pf->GsfTrk() && e->GsfTrk() &&
362 pf->GsfTrk() == e->GsfTrk()) {
363 isLepton = kTRUE;
364 break;
365 }
366 }
367 }
368
369 if (isLepton == kTRUE) continue;
370 ptSum += pf->Pt();
371
372 }
373 }
374 return ptSum;
375 }
376 //--------------------------------------------------------------------------------------------------
377 Double_t IsolationTools::BetaM(const TrackCol *tracks, const Muon *p, const Vertex *vertex,
378 Double_t ptMin, Double_t delta_z, Double_t extRadius,
379 Double_t intRadius){
380
381 if(!tracks) return 1.0;
382 if(tracks->GetEntries() <= 0) return 1.0;
383
384 double Pt_jets_X = 0. ;
385 double Pt_jets_Y = 0. ;
386 double Pt_jets_X_tot = 0. ;
387 double Pt_jets_Y_tot = 0. ;
388
389 for(int i=0;i<int(tracks->GetEntries());i++){
390 const Track *pTrack = tracks->At(i);
391
392 if(pTrack && p->TrackerTrk() &&
393 pTrack == p->TrackerTrk()) continue;
394
395 if(pTrack->Pt() <= ptMin) continue;
396
397 Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
398 if ( dr < extRadius && dr >= intRadius ) {
399 Pt_jets_X_tot += pTrack->Px();
400 Pt_jets_Y_tot += pTrack->Py();
401 double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
402 if(pDz < delta_z){
403 Pt_jets_X += pTrack->Px();
404 Pt_jets_Y += pTrack->Py();
405 }
406 }
407 }
408
409 if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
410 return sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot);
411
412 return 1.0;
413 }
414
415 //--------------------------------------------------------------------------------------------------
416 Double_t IsolationTools::BetaE(const TrackCol *tracks, const Electron *p, const Vertex *vertex,
417 Double_t ptMin, Double_t delta_z, Double_t extRadius,
418 Double_t intRadius){
419
420 if(!tracks) return 1.0;
421 if(tracks->GetEntries() <= 0) return 1.0;
422
423 double Pt_jets_X = 0. ;
424 double Pt_jets_Y = 0. ;
425 double Pt_jets_X_tot = 0. ;
426 double Pt_jets_Y_tot = 0. ;
427
428 for(int i=0;i<int(tracks->GetEntries());i++){
429 const Track *pTrack = tracks->At(i);
430
431 if(pTrack && p->TrackerTrk() &&
432 pTrack == p->TrackerTrk()) continue;
433
434 if(pTrack && p->GsfTrk() &&
435 pTrack == p->GsfTrk()) continue;
436
437 if(pTrack->Pt() <= ptMin) continue;
438
439 Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
440 if ( dr < extRadius && dr >= intRadius ) {
441 Pt_jets_X_tot += pTrack->Px();
442 Pt_jets_Y_tot += pTrack->Py();
443 double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
444 if(pDz < delta_z){
445 Pt_jets_X += pTrack->Px();
446 Pt_jets_Y += pTrack->Py();
447 }
448 }
449 }
450
451 if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
452 return sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot);
453
454 return 1.0;
455 }
456
457
458 // method added by F.Stoeckli: computes the track isolation with NO constrint on the OV-track compatibility
459 Double_t IsolationTools::TrackIsolationNoPV(const mithep::Particle* p, const BaseVertex* bsp,
460 Double_t extRadius,
461 Double_t intRadius,
462 Double_t ptLow,
463 Double_t etaStrip,
464 Double_t maxD0,
465 mithep::TrackQuality::EQuality quality,
466 const mithep::Collection<mithep::Track> *tracks,
467 UInt_t maxNExpectedHitsInner,
468 const mithep::DecayParticleCol *conversions) {
469
470 // loop over all tracks
471 Double_t tPt = 0.;
472 //std::cout<<" *** TrackIso:"<<std::endl;
473 for(UInt_t i=0; i<tracks->GetEntries(); ++i) {
474 const Track* t = tracks->At(i);
475 if(t->Pt()>1. && false) {
476 std::cout<<" "<<i<<" pt = "<<t->Pt()<<" ("<<ptLow<<")"<<std::endl;
477 std::cout<<" d0 = "<<fabs(t->D0Corrected( *bsp) )<<" ("<<maxD0<<")"<<std::endl;
478 //std::cout<<" conv ? "<<PhotonTools::MatchedConversion(t,conversions,bsp)<<std::endl;
479 std::cout<<" dR = "<<MathUtils::DeltaR(t->Mom(),p->Mom())<<" ("<<extRadius<<","<<intRadius<<")"<<std::endl;
480 std::cout<<" dEta = "<<fabs(t->Eta()-p->Eta())<<" ("<<etaStrip<<")"<<std::endl;
481 }
482 if ( t->Pt() < ptLow ) continue;
483 if ( ! t->Quality().Quality(quality) ) continue;
484 // only check for beamspot if available, otherwise ignore cut
485 if ( bsp && fabs(t->D0Corrected( *bsp) ) > maxD0) continue;
486 if (t->NExpectedHitsInner()>maxNExpectedHitsInner) continue;
487 if (conversions && PhotonTools::MatchedConversion(t,conversions,bsp)) continue;
488 Double_t dR = MathUtils::DeltaR(t->Mom(),p->Mom());
489 Double_t dEta = fabs(t->Eta()-p->Eta());
490 if(dR < extRadius && dR > intRadius && dEta > etaStrip) tPt += t->Pt();
491 }
492 return tPt;
493 }
494
495
496 Double_t IsolationTools::CiCTrackIsolation(const mithep::Photon* p,
497 const BaseVertex* theVtx,
498 Double_t extRadius,
499 Double_t intRadius,
500 Double_t ptLow,
501 Double_t etaStrip,
502 Double_t maxD0,
503 Double_t maxDZ,
504 const mithep::Collection<mithep::Track> *tracks,
505 unsigned int* worstVtxIndex,
506 const mithep::Collection<mithep::Vertex> *vtxs,
507 const mithep::Collection<mithep::Electron> *eles,
508 bool print,
509 double* ptmax, double* dRmax) {
510
511 UInt_t numVtx = 1;
512 const BaseVertex* iVtx = theVtx;
513 if( vtxs ) {
514 numVtx = vtxs->GetEntries();
515 if (numVtx > 0)
516 iVtx = vtxs->At(0);
517 else
518 return 0.;
519 }
520
521 // NEW for Electron T&P: need to remove the electron Gsf Track (applied if eles != NULL)
522 const Track* theGsfTrack = NULL;
523 if ( eles ) {
524 // find the electron that matches the Photon SC
525 for(UInt_t j=0; j<eles->GetEntries(); ++j) {
526 if ( eles->At(j)->SCluster() == p->SCluster() ) {
527 if( eles->At(j)->HasTrackerTrk() )
528 theGsfTrack = eles->At(j)->TrackerTrk();
529 break;
530 }
531 }
532 }
533
534 if(print) {
535 std::cout<<" Testing photon with"<<std::endl;
536 std::cout<<" Et = "<<p->Et()<<std::endl;
537 std::cout<<" Eta = "<<p->Eta()<<std::endl;
538 std::cout<<" Phi = "<<p->Phi()<<std::endl;
539 }
540
541 Double_t iIso = 0.;
542 Double_t maxIso = 0.;
543
544 if(worstVtxIndex)
545 *worstVtxIndex=0;
546
547 double t_ptmax = 0.;
548 double t_dRmax = 0.;
549
550 for(UInt_t i=0; i<numVtx; ++i) {
551
552 if(i>0) iVtx = vtxs->At(i);
553
554
555 if(print) {
556 std::cout<<" Vertex #"<<i<<std::endl;
557 std::cout<<" with X = "<<iVtx->X()<<std::endl;
558 std::cout<<" with Y = "<<iVtx->Y()<<std::endl;
559 std::cout<<" with Z = "<<iVtx->Z()<<std::endl;
560 }
561
562 Photon* phTemp = new Photon(*p);
563
564 // RESET CALO_POS!
565 phTemp->SetCaloPosXYZ(p->SCluster()->Point().X(),p->SCluster()->Point().Y(),p->SCluster()->Point().Z());
566
567 // compute the ph momentum with respect to this Vtx
568 //FourVectorM phMom = p->MomVtx(iVtx->Position());
569 FourVectorM phMom = phTemp->MomVtx(iVtx->Position());
570
571 delete phTemp;
572
573 if(print) {
574 std::cout<<" photon has changed to:"<<std::endl;
575 std::cout<<" Et = "<<phMom.Et()<<std::endl;
576 std::cout<<" eta = "<<phMom.Eta()<<std::endl;
577 std::cout<<" Phi = "<<phMom.Phi()<<std::endl;
578 }
579
580 iIso = 0.;
581 for(UInt_t j=0; j<tracks->GetEntries(); ++j) {
582 const Track* t = tracks->At(j);
583 if( theGsfTrack && t == theGsfTrack ) continue;
584
585 //Double_t dR = MathUtils::DeltaR(t->Mom(),p->Mom());
586 //Double_t dEta = TMath::Abs(t->Eta()-p->Eta());
587
588 Double_t dR = MathUtils::DeltaR(t->Mom(),phMom);
589 Double_t dEta = TMath::Abs(t->Eta()-phMom.Eta());
590
591 if(print && t->Pt()>1. && false) {
592 std::cout<<" passing track #"<<j<<std::endl;
593 std::cout<<" pt = "<<t->Pt()<<std::endl;
594 std::cout<<" eta = "<<t->Eta()<<std::endl;
595 std::cout<<" phi = "<<t->Phi()<<std::endl;
596 std::cout<<" d0 = "<<fabs(t->D0Corrected( *iVtx ))<<std::endl;
597 std::cout<<" dZ = "<<fabs(t->DzCorrected( *iVtx ))<<std::endl;
598 std::cout<<" dR = "<<dR<<std::endl;
599 std::cout<<" dEta = "<<dEta<<std::endl;
600 std::cout<<" vx = "<<t->X0()<<std::endl;
601 std::cout<<" vy = "<<t->Y0()<<std::endl;
602 std::cout<<" vz = "<<t->Z0()<<std::endl;
603 }
604
605 if ( t->Pt() < ptLow ) continue;
606 // only check for beamspot if available, otherwise ignore cut
607 if ( fabs(t->D0Corrected( *iVtx )) > maxD0) continue;
608 if ( fabs(t->DzCorrected( *iVtx )) > maxDZ) continue;
609
610
611 if(dR < extRadius && dR > intRadius && dEta > etaStrip) {
612 iIso += t->Pt();
613
614 if(t->Pt() > t_ptmax) {
615 t_ptmax=t->Pt();
616 t_dRmax=dR;
617 }
618
619 if(print && t->Pt()>1.) {
620 std::cout<<" passing track #"<<j<<std::endl;
621 std::cout<<" pt = "<<t->Pt()<<std::endl;
622 std::cout<<" eta = "<<t->Eta()<<std::endl;
623 std::cout<<" phi = "<<t->Phi()<<std::endl;
624 std::cout<<" d0 = "<<fabs(t->D0Corrected( *iVtx ))<<std::endl;
625 std::cout<<" dZ = "<<fabs(t->DzCorrected( *iVtx ))<<std::endl;
626 std::cout<<" dR = "<<dR<<std::endl;
627 std::cout<<" dEta = "<<dEta<<std::endl;
628 std::cout<<" vx = "<<t->X0()<<std::endl;
629 std::cout<<" vy = "<<t->Y0()<<std::endl;
630 std::cout<<" vz = "<<t->Z0()<<std::endl;
631 std::cout<<" new tIso = "<<iIso<<std::endl;
632 }
633 }
634 }
635 if ( iIso > maxIso ) {
636 maxIso = iIso;
637 if(worstVtxIndex)
638 *worstVtxIndex=i;
639 }
640 }
641
642 if(ptmax) (*ptmax)=t_ptmax;
643 if(dRmax) (*dRmax)=t_dRmax;
644
645 if(print) {
646 if(worstVtxIndex)
647 std::cout<<" max TrkIso is given by Vtx #"<<*worstVtxIndex<<" with an amount of tIso = "<<maxIso<<std::endl;
648 else
649 std::cout<<" max TrkIso is given by Vtx #0 with an amount of tIso = "<<maxIso<<std::endl;
650 }
651 return maxIso;
652 }