ATLAS Offline Software
Loading...
Searching...
No Matches
TrigTrackSelector.cxx
Go to the documentation of this file.
1
9
10
13
15
16#include <stdexcept>
17
22
23
24TrigTrackSelector::TrigTrackSelector( TrackFilter* selector, double radius, int selectPdgId, int selectParentPdgId ) :
25 TrackSelector(selector), m_id(0), m_xBeam(0), m_yBeam(0), m_zBeam(0),
26 m_correctTrkTracks(false),
27 m_radius(radius), m_selectPdgId(selectPdgId), m_selectParentPdgId(selectParentPdgId) {
28}
29
30
31
38
40 if ( p==nullptr ) return nullptr;
41 if ( MC::isElectron(p) || MC::isMuon(p) ) return nullptr; //don't want light leptons from eg tau decays - they are found directly
42 if ( p->absPdgId()==pdg_id ) {
43 return p; // recursive stopping conditions
44 }
45 auto vertex = p->prodVtx();
46 if ( vertex == nullptr ) {
47 return nullptr; // has no production vertex !!!
48 }
49 if ( vertex->nIncomingParticles() < 1 ) {
50 return nullptr; // recursive stopping conditions
51 }
52 for( unsigned ip = 0; ip < vertex->nIncomingParticles(); ip++ ) {
53 auto* in = vertex->incomingParticle(ip);
54 auto parent = fromAncestor( pdg_id, in);
55 if ( parent!=nullptr ) {
56 if (parent->absPdgId()==pdg_id) return parent;
57 }
58 }
59
60 return nullptr;
61}
62
63
72
73const xAOD::TruthParticle* TrigTrackSelector::fromAncestor( const std::vector<int>& ids, const xAOD::TruthParticle *p ) const {
74 if ( p==nullptr ) return nullptr;
75 if (MC::isElectron(p) || MC::isMuon(p)) return nullptr; //don't want light leptons from eg tau decays - they are found directly
76 for ( size_t i=ids.size() ; i-- ; ) {
77 if ( p->absPdgId()==ids[i] ) return p; // recursive stopping conditions
78 }
79
80 auto vertex = p->prodVtx();
81 if ( vertex == nullptr ) return nullptr; // has no production vertex !!!
82
83 if ( vertex->nIncomingParticles()<1 ) return nullptr; // recursive stopping conditions
84
85 for( unsigned ip = 0; ip < vertex->nIncomingParticles(); ip++ ) {
86 auto* in = vertex->incomingParticle(ip);
87 auto parent = fromAncestor( ids, in);
88 if ( parent!=nullptr ) {
89 for ( size_t i=ids.size() ; i-- ; ) {
90 if ( parent->absPdgId()==ids[i] ) return parent;
91 }
92 }
93 }
94
95 return nullptr;
96}
97
98
99
106
107// const xAOD::TruthParticle* TrigTrackSelector::fromAncestor(const int pdg_id, const xAOD::TruthParticle *p) const {
108// return fromAncestor( std::vector<int>(1,pdg_id), p );
109// }
110
111
112
113
114
115
116
117
118// add a TrackParticle
120
121 // do the track extraction stuff here....
122
123 static const int hpmap[20] = { 0, 1, 2, 7, 8, 9, 3, 4, 5, 6, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 };
124
125 if ( track ) {
126
127#ifdef TRKPARAMETERS_MEASUREDPERIGEE_H
128 const Trk::MeasuredPerigee* measPer = track->measuredPerigee();
129#else
130 const Trk::Perigee* measPer = track->measuredPerigee();
131#endif
132 // CLHEP::HepVector perigeeParams = measPer->parameters();
133
134 double pT = measPer->pT();
135 double eta = measPer->eta();
136 double phi = measPer->parameters()[Trk::phi0];
137 double z0 = measPer->parameters()[Trk::z0] + m_zBeam;
138 double d0 = measPer->parameters()[Trk::d0];
139
140 double theta = measPer->parameters()[Trk::theta];
141 double p = 1/measPer->parameters()[Trk::qOverP];
142
143 // AAARCH!!!!! the TrackParticle pT is NOT SIGNED!!!! ( I ask you! )
144 if ( measPer->parameters()[Trk::qOverP]<0 && pT>0 ) pT *= -1;
145
146#ifdef TRKPARAMETERS_MEASUREDPERIGEE_H
147 const Trk::ErrorMatrix err = measPer->localErrorMatrix();
148 double dtheta = err.error(Trk::theta);
149 double dqovp = err.error(Trk::qOverP);
150 double covthetaOvP = err.covValue(Trk::qOverP,Trk::theta);
151#else
152 double dtheta = std::sqrt((*measPer->covariance())(Trk::theta,Trk::theta));
153 double dqovp = std::sqrt((*measPer->covariance())(Trk::qOverP,Trk::qOverP));
154 double covthetaOvP = (*measPer->covariance())(Trk::qOverP,Trk::theta);
155#endif
156
157
158 double deta = 0.5*dtheta/(std::cos(0.5*theta)*std::cos(0.5*theta)*std::tan(0.5*theta));
159
160#ifdef TRKPARAMETERS_MEASUREDPERIGEE_H
161 double dphi = err.error(Trk::phi0);
162 double dz0 = err.error(Trk::z0);
163 double dd0 = err.error(Trk::d0);
164#else
165 double dphi = std::sqrt((*measPer->covariance())(Trk::phi0,Trk::phi0));
166 double dz0 = std::sqrt((*measPer->covariance())(Trk::z0,Trk::z0));
167 double dd0 = std::sqrt((*measPer->covariance())(Trk::d0,Trk::d0));
168#endif
169
170 double dpT = 0;
171
172
173 double sintheta = std::sin(theta);
174 double costheta = std::cos(theta);
175 double dpt2 = (p*p*sintheta)*(p*p*sintheta)*dqovp*dqovp + (p*costheta)*(p*costheta)*dtheta*dtheta - 2*(p*p*sintheta)*(p*costheta)*covthetaOvP;
176
177 if ( dpt2>0 ) dpT = std::sqrt( dpt2 );
178
179 // Check number of hits
180 // NB: a spacepoint is two offline "hits", so a pixel spacepoint is really
181 // 2 "hits" and an offline SCT "hit" is really a 1D cluster, so two intersecting
182 // stereo clusters making a spacepoint are two "hits"
183 const Trk::TrackSummary *summary = track->trackSummary();
184 int nBlayerHits = 2*summary->get(Trk::numberOfBLayerHits);
185 int nPixelHits = 2*summary->get(Trk::numberOfPixelHits);
186 int nSctHits = summary->get(Trk::numberOfSCTHits);
187 int nStrawHits = summary->get(Trk::numberOfTRTHits);
188 int nTrHits = summary->get(Trk::numberOfTRTHighThresholdHits);
189
190 int nSiHits = nPixelHits + nSctHits;
191 bool expectBL = false; // Not stored for Rec::TrackParticle
192
193 const Trk::FitQuality *quality = track->fitQuality();
194 double chi2 = quality->chiSquared();
195 double dof = quality->numberDoF();
196
197 unsigned bitmap = 0;
198
199
200 unsigned long id = (unsigned long)track;
201
202 for ( int ih=0 ; ih<20 ; ih++ ) {
203 if ( summary->isHit(Trk::DetectorType(ih)) ) bitmap |= ( 1<<hpmap[ih] );
204 }
205
208
209 // std::cout << "fetching author info :" << track->info().trackFitter() << ":"
210 // << track->info().dumpInfo() << ": bm 0x" << std::hex << bitmap << std::dec << std::endl;
211
212 int fitter = track->info().trackFitter();
213 std::string dumpinfo = track->info().dumpInfo();
214
215 int trackAuthor = -1;
216 if ( fitter>0 && fitter<Trk::TrackInfo::NumberOfTrackFitters ) {
217 if ( dumpinfo.find("TRTStandalone")!=std::string::npos) trackAuthor = 2;
218 else if ( dumpinfo.find("TRTSeededTrackFinder")!=std::string::npos) trackAuthor = 1;
219 else trackAuthor = 0;
220 }
221
222#if 0
223 std::cout << "\t\t\tSUTT TP track"
224 << "\teta=" << eta // << " +- " << (*trackitr)->params()->deta()
225 << "\tphi=" << phi // << " +- " << (*trackitr)->params()->dphi()
226 << "\tz0=" << z0
227 << "\tpT=" << pT // << "\t( " << 1/qoverp << ")"
228 << "\td0=" << d0
229 << "\tNsi=" << nSiHits
230 << "\tNtrt=" << nTrHits
231 << "\tNstr=" << nStrawHits
232 << "\tauthor=" << trackAuthor
233 << std::endl;
234#endif
235
236 // Create and save Track
237
238 TIDA::Track* t = new TIDA::Track(eta, phi, z0, d0, pT, chi2, dof,
239 deta, dphi, dz0, dd0, dpT,
240 nBlayerHits, nPixelHits, nSctHits, nSiHits,
241 nStrawHits, nTrHits, bitmap, 0,
242 trackAuthor, false, -1, -1,
243 expectBL, id) ;
244
245 // std::cout << "SUTT TP track " << *t << "\t0x" << std::hex << bitmap << std::dec << std::endl;
246
247 if ( !addTrack( t ) ){
248 delete t;
249 return false;
250 }
251 return true;
252
253 }
254 return false;
255}
256
257
258// extract all the tracks from a TrackParticle collection and add them
260
261 // std::cout << "\t\t\tSUTT \tTrackParticleContainer->size() = " << trigtracks->size() << std::endl;
262
263 Rec::TrackParticleContainer::const_iterator trackitr = trigtracks->begin();
264 Rec::TrackParticleContainer::const_iterator trackend = trigtracks->end();
265
266 while ( trackitr!=trackend ) {
267
268 selectTrack( *trackitr );
269
270 ++trackitr;
271
272 } // loop over tracks
273
274}
275
276
277
278// extract all the tracks from a TruthParticle collection and add them
280 // std::cout << "\t\t\tSUTT \tTrackParticleContainer->size() = " << trigtracks->size() << std::endl;
281
282 TruthParticleContainer::const_iterator trackitr = truthtracks->begin();
283 TruthParticleContainer::const_iterator trackend = truthtracks->end();
284
285 while ( trackitr!=trackend ) {
286
287 selectTrack( *trackitr );
288
289 ++trackitr;
290
291 } // loop over tracks
292
293}
294
295
299
300void TrigTrackSelector::truthBeamline( const xAOD::TruthParticleContainer* truthtracks, double& x0, double& y0 ) {
301
302 x0 = 0;
303 y0 = 0;
304
306
307 int Nx = 300;
308 int Ny = 300;
309
311 std::vector<double> xpos(Nx,0);
312 std::vector<double> ypos(Ny,0);
313
315 std::vector<int> xn(Nx,0);
316 std::vector<int> yn(Ny,0);
317
318 int xoffset = Nx/2;;
319 int yoffset = Ny/2;
320
321 double deltax = 3.0/Nx;
322 double deltay = 3.0/Ny;
323
325
326 xAOD::TruthParticleContainer::const_iterator trackitr = truthtracks->begin();
327 xAOD::TruthParticleContainer::const_iterator trackend = truthtracks->end();
328
329 for ( ; trackitr!=trackend ; ++trackitr ) {
330
331 const xAOD::TruthParticle* track = (*trackitr);
332
333 if ( !MC::isStable(track) || !track->hasProdVtx() ) continue;
334
336
337 double xp[3] = { track->prodVtx()->x(), track->prodVtx()->y(), track->prodVtx()->z() };
338
340
341 int ix = xp[0]/deltax + xoffset;
342 int iy = xp[1]/deltay + yoffset;
343
344 if ( ix<0 || ix>=Nx || iy<0 || iy>=Nx ) continue;
345
346 xpos[ix] += xp[0];
347 ypos[iy] += xp[1];
348
349 xn[ix]++;
350 yn[iy]++;
351
352 } // loop over tracks
353
355
356 int imx = 0;
357 int imy = 0;
358
359 for ( size_t i=0 ; i<xpos.size() ; i++ ) {
360 if ( xn[i]>xn[imx] ) imx = i;
361 if ( yn[i]>yn[imy] ) imy = i;
362 }
363
367 if ( xn[imx]>1 ) x0 = xpos[imx]/xn[imx];
368 if ( yn[imy]>1 ) y0 = ypos[imy]/yn[imy];
369
370}
371
372
373// extract all the tracks from a xAOD::TruthParticle collection and add them
375
376 xAOD::TruthParticleContainer::const_iterator trackitr = truthtracks->begin();
377 xAOD::TruthParticleContainer::const_iterator trackend = truthtracks->end();
378
380 double x0 = 0;
381 double y0 = 0;
382
383 truthBeamline( truthtracks, x0, y0 );
384
385
386 for ( ; trackitr!=trackend; ++trackitr) {
387
388
389 // Only select charged final state particles
390 double q = (*trackitr)->charge();
391
393 static const particleType ptype;
394 if ( q==-999 ) q = ptype.charge( (*trackitr)->pdgId() );
395
396 if (q == 0 || !MC::isStable(*trackitr) ) continue;
397
398 // If looking for tau parents, don't select mu or e children
399
400 // select based on the pdg of final state particle
401 bool gotPdgId = true;
402 if (m_selectPdgId!=0) gotPdgId = (*trackitr)->absPdgId()==m_selectPdgId;
403
404 // select based on the pdg of the parent or ancestor
405 bool gotParentPdgId = true;
406 if ( gotPdgId && m_selectParentPdgId!=0 ) gotParentPdgId = fromAncestor(m_selectParentPdgId, (*trackitr))!=nullptr;
409
410 if ( gotParentPdgId && gotPdgId ) selectTrack( *trackitr, x0, y0);
411
412 } // loop over tracks
413
414}
415
416
417
418// add a TruthParticle from a GenParticle - easy, bet it doesn't work
420
422 if ( !MC::isStable(track) ) return false;
423
425//AV Using memory to get some value is not a good idea. This is not a repruducible/portable way, but I leave it as is.
426 m_id = (unsigned long)(track.get());
427 bool sel;
428 sel = selectTrack( TruthParticle(track) );
429 m_id = 0;
430
431 return sel;
432
433}
434
435
436// add a TruthParticle
438
439 return selectTrack( &track );
440
441}
442
443
444// add a TruthParticle
446 TIDA::Track* t = makeTrack( track, m_id );
447 if ( t == 0 ) return false;
448 if ( !addTrack(t) ) {
449 delete t;
450 return false;
451 }
452 return true;
453}
454
455
456
457
458// add an xAOD::TruthParticle
459bool TrigTrackSelector::selectTrack( const xAOD::TruthParticle* track, double x0, double y0) {
460 if ( track ) {
461
462 if (!MC::isStable(track) ) return false;
463
465 const xAOD::TruthParticle* measPer = track;
466
467 double pT = measPer->pt();
468 double eta = measPer->eta();
469 double phi = measPer->phi();
470
471 // AAARCH!!!!! the TrackParticle pT is NOT SIGNED!!!! ( I ask you! )
472 if ( measPer->charge()<0 && pT>0 ) pT *= -1;
473 double q = track->charge();
474
475 static const particleType ptype;
476
478 if ( q==-999 ) q = ptype.charge( track->pdgId() );
479
481 if ( q==0 ) return 0;
482
486 // double xbeam = getBeamX(); // track->vx();
487 // double ybeam = getBeamY(); // track->vy();
488 // double zbeam = getBeamZ(); // track->vz();
489
490 if ( !track->hasProdVtx() ) return false;
491
493
494 double xp[3] = { measPer->prodVtx()->x(), measPer->prodVtx()->y(), measPer->prodVtx()->z() };
495 double xb[3] = { xp[0]-x0, xp[1]-y0, measPer->prodVtx()->z() };
496 double xd[3] = { 0, 0, 0 };
497
498 if ( track->hasDecayVtx() ) {
499 xd[0] = track->decayVtx()->x();
500 xd[1] = track->decayVtx()->y();
501 xd[2] = track->decayVtx()->z();
502 }
503
504 double rp = std::sqrt( xp[0]*xp[0] + xp[1]*xp[1] );
505 double rd = std::sqrt( xd[0]*xd[0] + xd[1]*xd[1] );
506
509 double theta = 2*std::atan( std::exp( -eta ) );
510 double z0 = xb[2] - (xb[0]*std::cos(phi) + xb[1]*std::sin(phi))/std::tan(theta);
511 double d0 = xb[1]*std::cos(phi) - xb[0]*std::sin(phi);
512
513 bool final_state = false;
514
526 const double inner_radius = m_radius;
527 const double outer_radius = m_radius;
528
529 if ( ( track->hasProdVtx() && rp<=inner_radius ) &&
530 ( !track->hasDecayVtx() || rd>outer_radius ) ) final_state = true;
531
532 if ( !final_state ) return false;
533
534 double deta = 0;
535 double dphi = 0;
536 double dz0 = 0;
537 double dd0 = 0;
538
539 double dpT = 0;
540
541 int nBlayerHits = 0;
542 int nPixelHits = 0;
543 int nSctHits = 0;
544 int nStrawHits = 0;
545 int nTrtHits = 0;
546
547 double chi2 = 0;
548 double dof = 0;
549
550 bool expectBL = false;
551
552 nSctHits += 0;
553 nPixelHits += 0;
554
556 int nSiHits = 0;
557
558 unsigned long id = (unsigned long)track;
559
560 unsigned bitmap = 0;
561
562 int trackAuthor = track->pdgId();
563 int uniqueID = HepMC::uniqueID(track);
564
565#if 0
566 std::cout << "\t\t\tSUTT TP track"
567 << "\teta=" << eta
568 << "\tphi=" << phi
569 << "\tz0=" << z0
570 << "\tpT=" << pT
571 << "\td0=" << d0
572 << "\tauthor=" << trackAuthor
573 << "\tVTX x " << xp[0]<< "\ty " << xp[1] << "\tz " << xp[2]
574 << std::endl;
575#endif
576
577 // Create and save Track
578
579 TIDA::Track* t = new TIDA::Track( eta, phi, z0, d0, pT, chi2, dof,
580 deta, dphi, dz0, dd0, dpT,
581 nBlayerHits, nPixelHits, nSctHits, nSiHits,
582 nStrawHits, nTrtHits, bitmap, 0,
583 trackAuthor, false, uniqueID, -1,
584 expectBL, id) ;
585
587 // std::cout << "SUTT TP track " << *t << "\t0x" << std::hex << bitmap << std::dec << std::endl;
588
589 // addTrack applies additional cuts using the Filter
590 if ( !addTrack( t ) ){
591 delete t;
592 return false;
593 }
594 }
595 return false;
596
597}
598
599
600
601
602// make a TIDA::Track from a GenParticle
604//AV Using memory to get some value is not a good idea. This is not a repruducible/portable way, but I leave it as is.
605 unsigned long id = (unsigned long)(track.get());
606 TruthParticle t = TruthParticle(track);
607 return makeTrack( &t, id );
608}
609
610// make a TIDA::Track from a TruthParticle
611TIDA::Track* TrigTrackSelector::makeTrack( const TruthParticle* track, unsigned long tid ) {
612
613 if ( track==0 ) return 0;
614 if ( !MC::isStable(track)) return 0;
615
616
617 double phi = track->phi();
618 double eta = track->eta();
619
623 double xp[3] = { 0, 0, 0 };
624
625 if ( track->genParticle()->production_vertex() ) {
626 xp[0] = track->genParticle()->production_vertex()->position().x();
627 xp[1] = track->genParticle()->production_vertex()->position().y();
628 xp[2] = track->genParticle()->production_vertex()->position().z();
629 }
630
631 // CHANGED BY JK - z0 with respect to (0,0)
632 // double z0 = xp[2];
633 double theta = 2*std::atan( exp( (-1)*eta ) );
634 double z0 = xp[2] - (xp[0]*std::cos(phi) + xp[1]*std::sin(phi))/std::tan(theta);
635
636 double xd[3] = { 0, 0, 0 };
637
638 if ( track->genParticle()->end_vertex() ) {
639 xd[0] = track->genParticle()->end_vertex()->position().x();
640 xd[1] = track->genParticle()->end_vertex()->position().y();
641 xd[2] = track->genParticle()->end_vertex()->position().z();
642 }
643
644 double rp = std::sqrt( xp[0]*xp[0] + xp[1]*xp[1] );
645 double rd = std::sqrt( xd[0]*xd[0] + xd[1]*xd[1] );
646
647
648 bool final_state = false;
649
661 const double inner_radius = m_radius;
662 const double outer_radius = m_radius;
663 if ( ( track->genParticle()->production_vertex() && rp<=inner_radius ) &&
664 ( track->genParticle()->end_vertex()==0 || rd>outer_radius ) ) final_state = true;
665
666
667 if ( !final_state ) return 0;
668
679
680 double q = track->charge();
681
682 static const particleType ptype;
683
685 if ( q==-999 ) q = ptype.charge( track->pdgId() );
686
688 if ( q==0 ) return 0;
689
690 double pT = q*track->pt();
691
692
693 double d0 = 0;
694
699
701
702
703 // CHANGED BY JK - d0 with respect to (0,0)
704 // d0 = q*rp*std::sin(phi);
705 d0 = xp[1]*std::cos(phi) - xp[0]*std::sin(phi);
706
707
709
710 double dz0 = 0;
711 double dd0 = 0;
712
713 correctToBeamline( z0, dz0, d0, dd0, theta, phi );
714
715
716
718
720 int author = track->pdgId();
721 int uniqueID = HepMC::uniqueID(track);
722
723
724 unsigned long id = (unsigned long)track;
725 if ( tid!=0 ) id = tid;
726
730
731
733
734 // std::cout << "\t\t\tSUTT Truth track"
735 // << "\teta=" << eta // << " +- " << (*trackitr)->params()->deta()
736 // << "\tphi=" << phi // << " +- " << (*trackitr)->params()->dphi()
737 // << "\tz0=" << z0
738 // << "\tpT=" << pT // << "\t( " << 1/qoverp << ")"
739 // << "\td0=" << d0
740 // << "\tauthor=" << author
741 // << std::endl;
742
743
744
745 TIDA::Track* t = new TIDA::Track(eta, phi, z0, d0, pT, 0, 0,
746 0, 0, 0, 0, 0,
747 0, 0, 0, 0,
748 0, 0, 0, 0,
749 author, false, uniqueID, -1,
750 false,
751 id ) ;
752
753 return t;
754
755}
756
757
758
759// add a Trk::Track
761
762 // do the track extraction stuff here....
763
764 static const int hpmap[20] = { 0, 1, 2, 7, 8, 9, 3, 4, 5, 6, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 };
765 // std::cout << __FILE__<< " " <<__LINE__<<std::endl;
766
767 if ( track ) {
768
769 // const Trk::Perigee* startPerigee = track->perigeeParameters();
770
771#ifdef TRKPARAMETERS_MEASUREDPERIGEE_H
772 const Trk::MeasuredPerigee* startPerigee = dynamic_cast<const Trk::MeasuredPerigee*>(track->perigeeParameters());
773 // const Trk::MeasuredPerigee* measPer = startPerigee; // just out of laziness
774#else
775 const Trk::Perigee* startPerigee = track->perigeeParameters();
776 const Trk::Perigee* measPer = startPerigee; // just out of laziness
777#endif
778
779
780 // CLHEP::HepVector perigeeParams = measPer->parameters();
781 // double pT = measPer->pT();
782 // double eta = measPer->eta();
783 // double phi = perigeeParams[Trk::phi0];
784 // double z0 = perigeeParams[Trk::z0];
785 // double d0 = perigeeParams[Trk::d0];
786 // // AAARCH!!!!! the TrackParticle pT is NOT SIGNED!!!! ( I ask you! )
787 // if ( perigeeParams[Trk::qOverP]<0 ) pT *= -1;
788 // std::cout <<pT1<<" pt1vspT "<<pT<<std::endl;
789
790 if (startPerigee){
791
792 double theta = startPerigee->parameters()[Trk::theta];
793 double p = 1/startPerigee->parameters()[Trk::qOverP];
794 double qOverPt = startPerigee->parameters()[Trk::qOverP]/std::sin(theta);
795 double charge = startPerigee->charge();
796 double eta = startPerigee->eta();
797 double phi = startPerigee->parameters()[Trk::phi0];
798 double z0 = startPerigee->parameters()[Trk::z0];
799 double d0 = startPerigee->parameters()[Trk::d0];
800 // double pT = (1./qOverPt)*(charge);
801 double pT = (1./qOverPt); // always use signed PT
802
803 if ( charge<0 && pT>0 ) pT *= -1;
804 if ( charge<0 && p>0 ) p *= -1;
805
806
807
808#ifdef TRKPARAMETERS_MEASUREDPERIGEE_H
809 const Trk::ErrorMatrix err = startPerigee->localErrorMatrix();
810 double dtheta = err.error(Trk::theta);
811 double dqovp = err.error(Trk::qOverP);
812 double covthetaOvP = err.covValue(Trk::qOverP,Trk::theta);
813
814 double dphi = err.error(Trk::phi0);
815 double dz0 = err.error(Trk::z0);
816 double dd0 = err.error(Trk::d0);
817#else
818 double dtheta = std::sqrt((*measPer->covariance())(Trk::theta,Trk::theta));
819 double dqovp = std::sqrt((*measPer->covariance())(Trk::qOverP,Trk::qOverP));
820 double covthetaOvP = (*measPer->covariance())(Trk::qOverP,Trk::theta);
821
822 double dphi = std::sqrt((*measPer->covariance())(Trk::phi0,Trk::phi0));
823 double dz0 = std::sqrt((*measPer->covariance())(Trk::z0,Trk::z0));
824 double dd0 = std::sqrt((*measPer->covariance())(Trk::d0,Trk::d0));
825#endif
826
827 double deta = 0.5*dtheta/(std::cos(0.5*theta)*std::cos(0.5*theta)*std::tan(0.5*theta));
828
829
830 if ( m_correctTrkTracks ) correctToBeamline( z0, dz0, d0, dd0, theta, phi );
831
832 double dpT = 0;
833
834
835 double sintheta = std::sin(theta);
836 double costheta = std::cos(theta);
837 double dpT2 = (p*p*sintheta)*(p*p*sintheta)*dqovp*dqovp + (p*costheta)*(p*costheta)*dtheta*dtheta - 2*(p*p*sintheta)*(p*costheta)*covthetaOvP;
838
839 if ( dpT2>0 ) dpT = std::sqrt( dpT2 );
840
841 // Check number of hits
842 // NB: a spacepoint is two offline "hits", so a pixel spacepoint is really
843 // 2 "hits" and an offline SCT "hit" is really a 1D cluster, so two intersetcting
844 // stereo clusters making a spacepoint are two "hits"
845 // const Trk::TrackSummary *summary = dynamic_cast<const Trk::TrackSummary*>(track->trackSummary());
846 //ToolHandle< Trk::ITrackSummaryTool > m_trackSumTool;
847 //m_trackSumTool = ToolHandle<Trk::ITrackSummaryTool>("Trk::TrackSummaryTool/InDetTrackSummaryTool");
848 //const Trk::TrackSummary* summary = NULL;
849 //summary = m_trackSumTool->createSummary(*track);
850
851 const Trk::TrackSummary * summary = track->trackSummary();
852 int nBlayerHits = 0;
853 int nPixelHits = 0;
854 int nSctHits = 0;
855 int nStrawHits = 0;
856 int nTrHits = 0;
857 int nSiHits = 0;
858 bool expectBL = false; // Not stored for Trk::Track
859 unsigned bitmap = 0;
860
861 if(summary==0){
862 std::cerr << "Could not create TrackSummary - Track will likely fail hits requirements" << std::endl;
863 }
864 else{
865 nBlayerHits = 2*summary->get(Trk::numberOfBLayerHits);
866 nPixelHits = 2*summary->get(Trk::numberOfPixelHits);
867 nSctHits = summary->get(Trk::numberOfSCTHits);
868 nStrawHits = summary->get(Trk::numberOfTRTHits);
869 nTrHits = summary->get(Trk::numberOfTRTHighThresholdHits);
870 nSiHits = nPixelHits + nSctHits;
871
872 for ( int ih=0 ; ih<20 ; ih++ ) {
873 if ( summary->isHit(Trk::DetectorType(ih)) ) bitmap |= ( 1<<hpmap[ih] );
874 }
875 }
876
877 unsigned long id = (unsigned long)track;
878 double chi2 = 0;
879 double dof = 0;
880 //const Trk::FitQuality *quality = dynamic_cast<const Trk::FitQuality*>(track->fitQuality());
881 const Trk::FitQuality *quality = (track->fitQuality());
882 if(quality==0) std::cerr << "Could not create FitQuality - Track will likely fail hits requirements" << std::endl;
883 else{
884 chi2 = quality->chiSquared();
885 dof = quality->numberDoF();
886 }
887
888 int trackAuthor = -1;
889
892
893 // std::cout << "fetching author info :" << track->info().trackFitter() << ":"
894 // << track->info().dumpInfo() << ": bm 0x" << std::hex << bitmap << std::dec << std::endl;
895
896 int fitter = track->info().trackFitter();
897 // std::string dumpinfo = track->info().dumpInfo();
898
899 if ( fitter>0 && fitter<Trk::TrackInfo::NumberOfTrackFitters ) {
900 if ((track->info().dumpInfo()).find("TRTStandalone") != std::string::npos) trackAuthor = 2;
901 else if ((track->info().dumpInfo()).find("TRTSeededTrackFinder") != std::string::npos) trackAuthor = 1;
902 else trackAuthor = 0;
903 }
904
905 #if 0
906 std::cout << "\t\t\tSUTT TP track"
907 << "\teta=" << eta // << " +- " << (*trackitr)->params()->deta()
908 << "\tphi=" << phi // << " +- " << (*trackitr)->params()->dphi()
909 << "\tz0=" << z0
910 << "\tpT=" << pT // << "\t( " << 1/qoverp << ")"
911 << "\td0=" << d0
912 << "\tNsi=" << nSiHits
913 << "\tNtrt=" << nTrHits
914 << "\tNstr=" << nStrawHits
915 << "\tauthor=" << trackAuthor
916 << std::endl;
917 #endif
918 // Create and save Track
919 TIDA::Track* t = new TIDA::Track(eta, phi, z0, d0, pT, chi2, dof,
920 deta, dphi, dz0, dd0, dpT,
921 nBlayerHits, nPixelHits, nSctHits, nSiHits,
922 nStrawHits, nTrHits, bitmap, 0,
923 trackAuthor, false, -1, -1,
924 expectBL, id) ;
925
926 if ( !addTrack( t ) ){
927 delete t;
928 return false;
929 }
930 return true;
931
932 //std::cout << "SUTT TP track " << *t << "\t0x" << std::hex << bitmap << std::dec << std::endl;
933 }
934 }
935
936 return false;
937}
938
939// extract all the tracks from a TrackCollection and add them
941
942 // std::cout << "\t\t\tSUTT \tTrackContainer->size() = " << trigtracks->size() << std::endl;
943
944 TrackCollection::const_iterator trackitr = trigtracks->begin();
945 TrackCollection::const_iterator trackend = trigtracks->end();
946
947 while ( trackitr!=trackend ) {
948 selectTrack( *trackitr );
949 ++trackitr;
950 } // loop over tracks
951
952}
953
954
955
957 // do the track extraction stuff here....
958
959 // static int hpmap[20] = { 0, 1, 2, 7, 8, 9, 3, 4, 5, 6, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 };
960
961 if ( track ) {
962
964 const xAOD::TrackParticle* measPer = track;
965
966 // CLHEP::HepVector perigeeParams = measPer->parameters();
967
968 double pT = measPer->pt();
969 double eta = measPer->eta();
970 double phi = measPer->phi0();
971 double z0 = measPer->z0() + measPer->vz();
972 double d0 = measPer->d0();
973
974 double theta = measPer->theta();
975
980 // if ( measPer->qOverP()==0 ) return false;
981 if ( measPer->qOverP()==0 ) throw std::runtime_error( "probable corrupted track - this should never happen" );
982 double p = 1/measPer->qOverP();
983
984 // AAARCH!!!!! the TrackParticle pT is NOT SIGNED!!!! ( I ask you! )
985 if ( measPer->qOverP()<0 && pT>0 ) pT *= -1;
986
987 double dtheta = std::sqrt(measPer->definingParametersCovMatrix()(Trk::theta,Trk::theta));
988 double dqovp = std::sqrt(measPer->definingParametersCovMatrix()(Trk::qOverP,Trk::qOverP));
989 double covthetaOvP = measPer->definingParametersCovMatrix()(Trk::qOverP,Trk::theta);
990
991 double deta = 0.5*dtheta/(std::cos(0.5*theta)*std::cos(0.5*theta)*std::tan(0.5*theta)); // ???? CHECK THIS <<--
992
993 double dphi = std::sqrt(measPer->definingParametersCovMatrix()(Trk::phi0,Trk::phi0));
994 double dz0 = std::sqrt(measPer->definingParametersCovMatrix()(Trk::z0,Trk::z0));
995 double dd0 = std::sqrt(measPer->definingParametersCovMatrix()(Trk::d0,Trk::d0));
996
997 double dpT = 0;
998
999
1001 // if ( m_xBeam!=0 || m_yBeam!=0 ) correctToBeamline( z0, dz0, d0, dd0, theta, phi );
1002
1003
1004 double sintheta = std::sin(theta);
1005 double costheta = std::cos(theta);
1006 double dpt2 = (p*p*sintheta)*(p*p*sintheta)*dqovp*dqovp + (p*costheta)*(p*costheta)*dtheta*dtheta - 2*(p*p*sintheta)*(p*costheta)*covthetaOvP;
1007
1008 if ( dpt2>0 ) dpT = std::sqrt( dpt2 );
1009
1010 // Check number of hits
1011 // NB: a spacepoint is two offline "hits", so a pixel spacepoint is really
1012 // 2 "hits" and an offline SCT "hit" is really a 1D cluster, so two intersetcting
1013 // stereo clusters making a spacepoint are two "hits"
1014
1015 uint8_t sum_nBlayerHits = 0;
1016 track->summaryValue( sum_nBlayerHits, xAOD::numberOfInnermostPixelLayerHits);
1017 int nBlayerHits = 2*sum_nBlayerHits;
1018
1019 uint8_t sum_nPixelHits = 0;
1020 track->summaryValue( sum_nPixelHits, xAOD::numberOfPixelHits);
1021 int nPixelHits = 2*sum_nPixelHits;
1022
1023 uint8_t sum_nSctHits = 0;
1024 track->summaryValue( sum_nSctHits, xAOD::numberOfSCTHits);
1025 int nSctHits = sum_nSctHits;
1026
1027 uint8_t sum_nStrawHits = 0;
1028 track->summaryValue( sum_nStrawHits, xAOD::numberOfTRTHits);
1029 int nStrawHits = sum_nStrawHits;
1030
1031 uint8_t sum_nTrtHits = 0;
1032 track->summaryValue( sum_nTrtHits, xAOD::numberOfTRTHighThresholdHits);
1033 int nTrtHits = sum_nTrtHits;
1034
1035
1036 uint8_t sum_expectBL = 0;
1037 track->summaryValue( sum_expectBL, xAOD::expectInnermostPixelLayerHit);
1038 bool expectBL = ( sum_expectBL ? true : false );
1039
1041
1042 uint8_t sum_sctholes = 0;
1043 track->summaryValue( sum_sctholes, xAOD::numberOfSCTHoles);
1044
1045 uint8_t sum_pixholes = 0;
1046 track->summaryValue( sum_pixholes, xAOD::numberOfPixelHoles);
1047
1051
1052 nSctHits += 1000*sum_sctholes;
1053 nPixelHits += 1000*sum_pixholes;
1054
1056 int nSiHits = nPixelHits + nSctHits;
1057
1059
1060 double chi2 = track->chiSquared();
1061 double dof = track->numberDoF();
1062
1063 unsigned long id = (unsigned long)track;
1064
1065 unsigned bitmap = track->hitPattern();
1066
1067
1068
1069 double xbeam = track->vx();
1070 double ybeam = track->vy();
1071 double zbeam = track->vz();
1072
1073 if ( xbeam!=getBeamX() || ybeam!=getBeamY() || zbeam!=getBeamZ() ) setBeamline( xbeam, ybeam, zbeam );
1074
1075 int trackAuthor = 0;
1076
1077 int fitter = track->trackFitter();
1078 std::bitset<xAOD::NumberOfTrackRecoInfo> patternrec = track->patternRecoInfo();
1079
1080 //int icount = 0;<- never used if section below is commented
1081 for ( unsigned ipr=patternrec.size() ; ipr-- ; ) {
1082 if ( patternrec[ipr] ) {
1083 //icount++; <- never used if section below is commented
1084 trackAuthor |= (ipr >> 16);
1085 // static bool first = true;
1086 // if ( first && icount>1 ) {
1087 // std::cerr << "more than one pattern rec strategy " << ipr << "\t(suppressing further output)" << std::endl;
1088 // first = false;
1089 // }
1090 }
1091 }
1092
1093 trackAuthor |= fitter;
1094
1095 // if ( fitter>0 && fitter<Trk::TrackInfo::NumberOfTrackFitters ) {
1096 // if ( dumpinfo.find("TRTStandalone")!=std::string::npos) trackAuthor = 2;
1097 // else if ( dumpinfo.find("TRTSeededTrackFinder")!=std::string::npos) trackAuthor = 1;
1098 // else trackAuthor = 0;
1099 // }
1100
1101#if 0
1102 std::cout << "\t\t\tSUTT TP track"
1103 << "\teta=" << eta // << " +- " << (*trackitr)->params()->deta()
1104 << "\tphi=" << phi // << " +- " << (*trackitr)->params()->dphi()
1105 << "\tz0=" << z0
1106 << "\tpT=" << pT // << "\t( " << 1/qoverp << ")"
1107 << "\td0=" << d0
1108 << "\tNsi=" << nSiHits
1109 // << "\tNtrt=" << nTrtHits
1110 // << "\tNstr=" << nStrawHits
1111 << "\tfitter=" << fitter
1112 << "\tauthor=" << trackAuthor
1113 << "\tVTX x " << track->vx() << "\ty " << track->vy() << "\tz " << track->vz()
1114 << std::endl;
1115#endif
1116
1117 // Create and save Track
1118
1119 TIDA::Track* t = new TIDA::Track( eta, phi, z0, d0, pT, chi2, dof,
1120 deta, dphi, dz0, dd0, dpT,
1121 nBlayerHits, nPixelHits, nSctHits, nSiHits,
1122 nStrawHits, nTrtHits, bitmap, 0,
1123 trackAuthor, false, -1, -1,
1124 expectBL, id) ;
1125
1126 // std::cout << "SUTT TP track " << *t << "\t0x" << std::hex << bitmap << std::dec << std::endl;
1127
1128 if ( !addTrack( t ) ){
1129 delete t;
1130 return false;
1131 }
1132 return true;
1133
1134 }
1135 return false;
1136
1137
1138}
1139
1140
1142 // std::cout << "\t\t\tSUTT \tTrackContainer->size() = " << trigtracks->size() << std::endl;
1145 while ( trackitr!=trackend ) {
1146 selectTrack( *trackitr );
1147 ++trackitr;
1148 } // loop over tracks
1149}
1150
1151
1152
1157 // std::cout << "\t\t\tSUTT \tTrackContainer->size() = " << trigtracks->size() << std::endl;
1158 while ( trackitr!=trackend ) {
1159 selectTrack( *trackitr );
1160 ++trackitr;
1161 } // loop over tracks
1162}
1163
1164
1166 for( const auto& track : tracks ) selectTrack( *track );
1167}
1168
1169
1170
1171
1172
1173void TrigTrackSelector::correctToBeamline( double& z0, double& dz0,
1174 double& d0, double& dd0,
1175 double theta, double phi ) {
1176
1178
1179 // if ( m_first ) {
1180 // if ( m_xBeam==0 && m_yBeam==0 ) {
1181 // std::cerr << "TrigTrackSelector::correctToBeamline() WARNING -- Beamline set to (0,0) -- WARNING" << std::endl;
1182 // }
1183 // else {
1184 // std::cout << "TrigTrackSelector::correctToBeamline() Beamline set to " << m_xBeam << " " << m_yBeam << std::endl;
1185 // }
1186 // m_first = false;
1187 // }
1188
1189
1190 // double theta = 2*std::atan( exp( (-1)*eta ) );
1191 double z0t = z0 + ((std::cos(phi)*m_xBeam + std::sin(phi)*m_yBeam)/std::tan(theta));
1192 double a0t = d0 + std::sin(phi)*m_xBeam - std::cos(phi)*m_yBeam;
1193
1195 double dz0t = dz0 + ((std::cos(phi)*m_xBeam + std::sin(phi)*m_yBeam)/std::tan(theta));
1196 double da0t = dd0 + std::sin(phi)*m_xBeam - std::cos(phi)*m_yBeam;
1197
1198 z0 = z0t;
1199 d0 = a0t;
1200
1201 dz0 = dz0t;
1202 dd0 = da0t;
1203}
1204
Scalar eta() const
pseudorapidity method
Scalar phi() const
phi method
Scalar theta() const
theta method
double charge(const T &p)
Definition AtlasPID.h:1003
ATLAS-specific HepMC functions.
ReadCards * rp
DataVector< Trk::Track > TrackCollection
This typedef represents a collection of Trk::Track objects.
DataModel_detail::const_iterator< DataVector > const_iterator
Definition DataVector.h:838
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
virtual bool addTrack(TIDA::Track *t, bool(*f)(const TIDA::Track *)=0)
void setBeamline(double x, double y, double z=0)
const xAOD::TruthParticle * fromAncestor(const int pdg_id, const xAOD::TruthParticle *p) const
recursive functions to identify whether a particle comes from some particle of a specific PDG ID,...
double getBeamY() const
void truthBeamline(const xAOD::TruthParticleContainer *truthtracks, double &x0, double &y0)
extract all the tracks from a xAOD::TruthParticle collection and histogram the x and y production coo...
double getBeamX() const
TrigTrackSelector(TrackFilter *selector)
use a radius of 47 mm corresponding to the Run 1 pixel inner radius For the IBL it should be 32 mm,...
TIDA::Track * makeTrack(HepMC::ConstGenParticlePtr track)
double getBeamZ() const
void correctToBeamline(double &z0, double &dz0, double &d0, double &dd0, double theta, double phi)
bool selectTrack(const Rec::TrackParticle *track)
neater code to make use of vector function also for a single ancestor pdgid, instead of the full code...
static const double s_default_radius
NB: This was 47 for Run 2, but with the addition of the IBL it should be 32 It was kept at 47 for all...
void selectTracks(const Rec::TrackParticleContainer *trigtracks)
std::vector< ElementLink< xAOD::TrackParticleContainer > > TrackParticleLinks_t
Class to represent and store fit qualities from track reconstruction in terms of and number of degre...
Definition FitQuality.h:97
int numberDoF() const
returns the number of degrees of freedom of the overall track or vertex fit as integer
Definition FitQuality.h:60
double chiSquared() const
returns the of the overall track fit
Definition FitQuality.h:56
double eta() const
Access method for pseudorapidity - from momentum.
double charge() const
Returns the charge.
double pT() const
Access method for transverse momentum.
A summary of the information contained by a track.
double charge(int id) const
float z0() const
Returns the parameter.
float theta() const
Returns the parameter, which has range 0 to .
const ParametersCovMatrix_t definingParametersCovMatrix() const
Returns the 5x5 symmetric matrix containing the defining parameters covariance matrix.
float vz() const
The z origin for the parameters.
float d0() const
Returns the parameter.
float qOverP() const
Returns the parameter.
virtual double pt() const override final
The transverse momentum ( ) of the particle.
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
float phi0() const
Returns the parameter, which has range to .
virtual double pt() const override final
The transverse momentum ( ) of the particle.
const TruthVertex_v1 * prodVtx() const
The production vertex of this particle.
double charge() const
Physical charge.
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
virtual double phi() const override final
The azimuthal angle ( ) of the particle.
float z() const
Vertex longitudinal distance along the beam line form the origin.
float y() const
Vertex y displacement.
float x() const
Vertex x displacement.
Definition zbeam.h:24
double chi2(TH1 *h0, TH1 *h1)
int uniqueID(const T &p)
HepMC3::ConstGenParticlePtr ConstGenParticlePtr
Definition GenParticle.h:20
bool isElectron(const T &p)
bool isStable(const T &p)
Identify if the particle is stable, i.e. has not decayed.
bool isMuon(const T &p)
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
Perigee MeasuredPerigee
@ phi0
Definition ParamDefs.h:65
@ theta
Definition ParamDefs.h:66
@ qOverP
perigee
Definition ParamDefs.h:67
@ d0
Definition ParamDefs.h:63
@ z0
Definition ParamDefs.h:64
DetectorType
enumerates the various detector types currently accessible from the isHit() method.
@ numberOfPixelHits
number of pixel layers on track with absence of hits
@ numberOfBLayerHits
these are the hits in the 0th pixel layer?
@ numberOfTRTHighThresholdHits
total number of TRT hits which pass the high threshold
TrackParticle_v1 TrackParticle
Reference the current persistent version:
TruthParticle_v1 TruthParticle
Typedef to implementation.
TrackParticleContainer_v1 TrackParticleContainer
Definition of the current "TrackParticle container version".
@ expectInnermostPixelLayerHit
Do we expect a 0th-layer barrel hit for this track?
@ numberOfPixelHoles
number of pixel layers on track with absence of hits [unit8_t].
@ numberOfTRTHits
number of TRT hits [unit8_t].
@ numberOfSCTHits
number of hits in SCT [unit8_t].
@ numberOfInnermostPixelLayerHits
these are the hits in the 0th pixel barrel layer
@ numberOfPixelHits
these are the pixel hits, including the b-layer [unit8_t].
@ numberOfTRTHighThresholdHits
number of TRT hits which pass the high threshold (only xenon counted) [unit8_t].
@ numberOfSCTHoles
number of SCT holes [unit8_t].
TruthParticleContainer_v1 TruthParticleContainer
Declare the latest version of the truth particle container.