ATLAS Offline Software
Loading...
Searching...
No Matches
SiSpacePointsSeedMaker_ATLxk.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
5
7// Implementation file for class SiSpacePointsSeedMaker_ATLxk
9// AlgTool used for TRT_DriftCircleOnTrack object production
11// Version 1.0 21/04/2004 I.Gavrilenko
13
15
18
19#include "TrkTrack/Track.h"
21
25#include <iomanip>
26#include <ostream>
27
31
33(const std::string& t,const std::string& n,const IInterface* p)
34 : base_class(t, n, p),
35 m_thistSvc("THistSvc", n),
36 m_outputTree(nullptr),
37 m_treeName(""),
38 m_treeFolder("/valNtuples/")
39{
40}
41
45
47{
48 StatusCode sc = AlgTool::initialize();
49
51 ATH_CHECK(m_spacepointsSCT.initialize(m_sct));
53
55 if (not m_beamSpotKey.empty()) {
56 ATH_CHECK(m_beamSpotKey.initialize());
57 }
58
59 ATH_CHECK( m_fieldCondObjInputKey.initialize() );
60
62 ATH_CHECK( m_prdToTrackMap.initialize( !m_prdToTrackMap.key().empty()));
63
67
68 if (msgLvl(MSG::DEBUG)) {
69 EventData data;
71 data.nprint=0;
72 dump(data, msg(MSG::DEBUG));
73 }
74
75 if (m_writeNtuple) {
76 //coverity[MISSING_LOCK]
77 ATH_CHECK( m_thistSvc.retrieve() );
78
79 m_treeName = (std::string("SeedTree_")+name());
80 std::replace( m_treeName.begin(), m_treeName.end(), '.', '_' );
81
82 m_outputTree = new TTree( m_treeName.c_str() , "SeedMakerValTool");
83
84 m_outputTree->Branch("eventNumber", &m_eventNumber);
85 m_outputTree->Branch("d0", &m_d0);
86 m_outputTree->Branch("z0", &m_z0);
87 m_outputTree->Branch("pt", &m_pt);
88 m_outputTree->Branch("eta", &m_eta);
89 m_outputTree->Branch("x1", &m_x1);
90 m_outputTree->Branch("x2", &m_x2);
91 m_outputTree->Branch("x3", &m_x3);
92 m_outputTree->Branch("y1", &m_y1);
93 m_outputTree->Branch("y2", &m_y2);
94 m_outputTree->Branch("y3", &m_y3);
95 m_outputTree->Branch("z1", &m_z1);
96 m_outputTree->Branch("z2", &m_z2);
97 m_outputTree->Branch("z3", &m_z3);
98 m_outputTree->Branch("r1", &m_r1);
99 m_outputTree->Branch("r2", &m_r2);
100 m_outputTree->Branch("r3", &m_r3);
101 m_outputTree->Branch("quality", &m_quality);
102 m_outputTree->Branch("seedType", &m_type);
103 m_outputTree->Branch("givesTrack", &m_givesTrack);
104 m_outputTree->Branch("dzdr_b", &m_dzdr_b);
105 m_outputTree->Branch("dzdr_t", &m_dzdr_t);
106 m_outputTree->Branch("track_pt", &m_trackPt);
107 m_outputTree->Branch("track_eta", &m_trackEta);
108
109 std::string fullTreeName = m_treeFolder + m_treeName;
110
111 ATH_CHECK( m_thistSvc->regTree( fullTreeName, m_outputTree ) );
112
113 }
114
115 return sc;
116}
117
121
123{
124 return AlgTool::finalize();
125}
126
130
131void InDet::SiSpacePointsSeedMaker_ATLxk::newEvent(const EventContext& ctx, EventData& data, int iteration) const
132{
134 if (not data.initialized) initializeEventData(data);
135
136 data.trigger = false;
137 if (!m_pixel && !m_sct) return;
138
140 data.iteration = iteration;
141 if (iteration <=0) data.iteration = 0;
143 erase(data);
144 data.dzdrmin = m_dzdrmin0;
145 data.dzdrmax = m_dzdrmax0;
146 data.maxScore = m_maxScore;
147
149 if (data.iteration == 0) {
150 if (not m_beamSpotKey.empty()) {
151 buildBeamFrameWork(data);
152 }
153
155 double magField[3]{0,0,0};
156 double globalPos[3] ={10.,10.,0.};
157
158 MagField::AtlasFieldCache fieldCache;
161 const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
162 if (fieldCondObj == nullptr) {
163 ATH_MSG_ERROR("SiSpacePointsSeedMaker_ATLxk: Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
164 return;
165 }
166 fieldCondObj->getInitializedCache (fieldCache);
167
168 if (fieldCache.solenoidOn()) {
170 fieldCache.getFieldZR(globalPos,magField);
181 data.K = 2.f/(300.f*float(magField[2]));
182 } else {
183 data.K = 2.f/(300.f* 5.f );
184 }
189 data.ipt2K = m_ipt2/(data.K*data.K);
191 data.ipt2C = m_ipt2*m_COF;
192 data.COFK = m_COF*(data.K*data.K);
193
195 data.i_spforseed = data.l_spforseed.begin();
196 // Set the seed multiplicity strategy of the event data to the one configured
197 // by the user for strip seeds
198 data.maxSeedsPerSP = m_maxOneSizeSSS;
199 data.keepAllConfirmedSeeds = m_alwaysKeepConfirmedStripSeeds;
200
201 }
202 else {
203 data.r_first = 0;
204 // Set the seed multiplicity strategy of the event data to the one configured
205 // by the user for pixel seeds
206 data.maxSeedsPerSP = m_maxOneSizePPP;
207 data.keepAllConfirmedSeeds = m_alwaysKeepConfirmedPixelSeeds;
208
210 fillLists(data);
211 return;
212 }
213
215
216 data.checketa = data.dzdrmin > 1.;
217
219 float oneOverBinSizeR = 1.f/m_binSizeR;
220 int maxBinR = m_nBinsR-1;
221
228 for (int i=0; i<data.nr; ++i) {
229 int n = data.r_index[i];
230 data.r_map[n] = 0;
231 data.r_Sorted[n].clear();
232 }
233 data.ns = data.nr = 0;
234
236 SG::ReadHandle<Trk::PRDtoTrackMap> prd_to_track_map;
237 const Trk::PRDtoTrackMap *prd_to_track_map_cptr = nullptr;
238 if (!m_prdToTrackMap.key().empty()) {
240 if (!prd_to_track_map.isValid()) {
241 ATH_MSG_ERROR("Failed to read PRD to track association map: " << m_prdToTrackMap.key());
242 }
243 prd_to_track_map_cptr = prd_to_track_map.cptr();
244 }
245
249
251 data.r_first = 0;
252
254 if (m_pixel) {
255
257 if (spacepointsPixel.isValid()) {
259 for (const SpacePointCollection* spc: *spacepointsPixel) {
260 for (const Trk::SpacePoint* sp: *spc) {
261
264 if ((prd_to_track_map_cptr && isUsed(sp,*prd_to_track_map_cptr)) || sp->r() > m_r_rmax) continue;
265
269 static_cast<const InDetDD::SiDetectorElement*>(sp->clusterList().first->detectorElement());
270 if (!de || de->isDBM()) continue;
271
278 if (!sps) continue;
279
282 int radiusBin = static_cast<int>(sps->radius()*oneOverBinSizeR);
284 if (radiusBin>maxBinR) radiusBin = maxBinR;
285
287 data.r_Sorted[radiusBin].push_back(sps);
289 ++data.r_map[radiusBin];
292 if (data.r_map[radiusBin]==1) data.r_index[data.nr++] = radiusBin;
294 if (radiusBin > data.r_first) data.r_first = radiusBin;
296 ++data.ns;
297 }
298 }
299 }
300 ++data.r_first;
301 }
302
304 if (m_sct) {
305
307 if (spacepointsSCT.isValid()) {
308
309 for (const SpacePointCollection* spc: *spacepointsSCT) {
310 for (const Trk::SpacePoint* sp: *spc) {
313 if ((prd_to_track_map_cptr && isUsed(sp,*prd_to_track_map_cptr)) || sp->r() > m_r_rmax) continue;
314
317 if (!sps) continue;
318
321 int radiusBin = static_cast<int>(sps->radius()*oneOverBinSizeR);
322 if (radiusBin>maxBinR) radiusBin = maxBinR;
324 data.r_Sorted[radiusBin].push_back(sps);
326 ++data.r_map[radiusBin];
328 if (data.r_map[radiusBin]==1) data.r_index[data.nr++] = radiusBin;
330 ++data.ns;
331 }
332 }
333 }
334
336 if (m_useOverlap && !data.checketa) {
337
339 if (spacepointsOverlap.isValid()) {
340
341 for (const Trk::SpacePoint* sp: *spacepointsOverlap) {
343 if ((prd_to_track_map_cptr && isUsed(sp, *prd_to_track_map_cptr)) || sp->r() > m_r_rmax) continue;
344
347 if (!sps) continue;
348
350 int radiusBin = static_cast<int>(sps->radius()*oneOverBinSizeR);
351 if (radiusBin>maxBinR) radiusBin = maxBinR;
353 data.r_Sorted[radiusBin].push_back(sps);
355 ++data.r_map[radiusBin];
357 if (data.r_map[radiusBin]==1) data.r_index[data.nr++] = radiusBin;
359 ++data.ns;
360 }
361 }
362 }
363 }
364
366 if (iteration < 0) data.r_first = 0;
367
370 fillLists(data);
371}
372
376
379(const EventContext& ctx, EventData& data,
380 const std::vector<IdentifierHash>& vPixel, const std::vector<IdentifierHash>& vSCT) const
381{
382 if (not data.initialized) initializeEventData(data);
383
384 data.iteration = 0;
385 data.trigger = false;
386 erase(data);
387 if (!m_pixel && !m_sct) return;
388
389 data.dzdrmin = m_dzdrmin0;
390 data.dzdrmax = m_dzdrmax0;
391 data.maxScore = m_maxScore;
392
393 if (not m_beamSpotKey.empty()) {
394 buildBeamFrameWork(data);
395 }
396
397 double magField[3]{0,0,0};
398 double globalPos[3] ={10.,10.,0.};
399
400 MagField::AtlasFieldCache fieldCache;
401 // Get field cache object
403 const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
404
405 if (fieldCondObj == nullptr) {
406 ATH_MSG_ERROR("SiSpacePointsSeedMaker_ATLxk: Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
407 return;
408 }
409 fieldCondObj->getInitializedCache (fieldCache);
410
411 // initialise pt cuts and conversion factors
412 if (fieldCache.solenoidOn()) {
413 fieldCache.getFieldZR(globalPos,magField);
414
415 data.K = 2.f/(300.f*float(magField[2]));
416 } else {
417 data.K = 2.f/(300.f* 5.f );
418 }
419
420 data.ipt2K = m_ipt2/(data.K*data.K);
421 data.ipt2C = m_ipt2*m_COF;
422 data.COFK = m_COF*(data.K*data.K);
423
424 data.i_spforseed = data.l_spforseed.begin();
425
426 float oneOverBinSizeR = 1.f/m_binSizeR;
427 int maxBinR = m_nBinsR-1;
428
429 data.r_first = 0;
430 data.checketa = false;
431
432 for (int i=0; i<data.nr; ++i) {
433 int n = data.r_index[i];
434 data.r_map[n] = 0;
435 data.r_Sorted[n].clear();
436 }
437 data.ns = data.nr = 0;
438
439 // Get pixels space points containers from store gate
440 //
441 if (m_pixel && !vPixel.empty()) {
442
444 if ( spacepointsPixel.isValid() ) {
445
446 // Loop through all trigger collections
447 //
448 for (const IdentifierHash& l: vPixel) {
449 const auto *w = spacepointsPixel->indexFindPtr(l);
450 if (w==nullptr) continue;
451 for (const Trk::SpacePoint* sp: *w) {
452 float r = sp->r();
453 if (r > m_r_rmax) continue;
455 int ir = static_cast<int>(sps->radius()*oneOverBinSizeR);
456 if (ir>maxBinR) ir = maxBinR;
457 data.r_Sorted[ir].push_back(sps);
458 ++data.r_map[ir];
459 if (data.r_map[ir]==1) data.r_index[data.nr++] = ir;
460 ++data.ns;
461 }
462 }
463 }
464 }
465
466 // Get sct space points containers from store gate
467 //
468 if (m_sct && !vSCT.empty()) {
469
471 if (spacepointsSCT.isValid()) {
472
473 // Loop through all trigger collections
474 //
475 for (const IdentifierHash& l: vSCT) {
476 const auto *w = spacepointsSCT->indexFindPtr(l);
477 if (w==nullptr) continue;
478 for (const Trk::SpacePoint* sp: *w) {
479 float r = sp->r();
480 if (r > m_r_rmax) continue;
482 int ir = static_cast<int>(sps->radius()*oneOverBinSizeR);
483 if (ir>maxBinR) ir = maxBinR;
484 data.r_Sorted[ir].push_back(sps);
485 ++data.r_map[ir];
486 if (data.r_map[ir]==1) data.r_index[data.nr++] = ir;
487 ++data.ns;
488 }
489 }
490 }
491 }
492 fillLists(data);
493}
494
498
500(const EventContext& ctx, EventData& data,
501 const std::vector<IdentifierHash>& vPixel, const std::vector<IdentifierHash>& vSCT,
502 const IRoiDescriptor& IRD) const
503{
504 constexpr float twoPi = 2.*M_PI;
505
506 if (not data.initialized) initializeEventData(data);
507
508 newRegion(ctx, data, vPixel, vSCT);
509 data.trigger = true;
510
511 double dzdrmin = 1./std::tan(2.*std::atan(std::exp(-IRD.etaMinus())));
512 double dzdrmax = 1./std::tan(2.*std::atan(std::exp(-IRD.etaPlus ())));
513
514 data.zminB = IRD.zedMinus()-data.zbeam[0]; // min bottom Z
515 data.zmaxB = IRD.zedPlus ()-data.zbeam[0]; // max bottom Z
516 data.zminU = data.zminB+550.f*float(dzdrmin);
517 data.zmaxU = data.zmaxB+550.f*float(dzdrmax);
518 float fmax = IRD.phiPlus ();
519 float fmin = IRD.phiMinus();
520 if (fmin > fmax) fmin -= twoPi;
521 data.ftrig = (fmin+fmax)*.5f;
522 data.ftrigW = (fmax-fmin)*.5f;
523}
524
529
530void InDet::SiSpacePointsSeedMaker_ATLxk::find2Sp(EventData& data, const std::list<Trk::Vertex>& lv) const
531{
532 if (not data.initialized) initializeEventData(data);
533
534 data.zminU = m_zmin;
535 data.zmaxU = m_zmax;
536
537 int mode = 0;
538 if (lv.begin()!=lv.end()) mode = 1;
539 bool newv = newVertices(data, lv);
540
541 if (newv || !data.state || data.nspoint!=2 || data.mode!=mode || data.nlist) {
542 data.i_seede_Pro = data.l_seeds_Pro.begin();
543 data.state = 1;
544 data.nspoint = 2;
545 data.nlist = 0;
546 data.mode = mode;
547 data.endlist = true;
548 data.fvNmin = 0;
549 data.fNmin = 0;
550 data.zMin = 0;
551 production2Sp(data);
552 }
553 data.i_seed_Pro = data.l_seeds_Pro.begin();
554
555
556 if (msgLvl(MSG::DEBUG)) {
557 data.nprint=1;
558 dump(data, msg(MSG::DEBUG));
559 }
560}
561
566
567void InDet::SiSpacePointsSeedMaker_ATLxk::find3Sp(const EventContext&, EventData& data, const std::list<Trk::Vertex>& lv) const
568{
569 if (not data.initialized) initializeEventData(data);
570
572 data.zminU = m_zmin;
573 data.zmaxU = m_zmax;
575 int mode = 2;
576 if (lv.begin()!=lv.end()) mode = 3;
581 bool newv = newVertices(data, lv);
583 if (newv || !data.state || data.nspoint!=3 || data.mode!=mode || data.nlist) {
584 data.i_seede_Pro = data.l_seeds_Pro.begin();
585 data.state = 1;
586 data.nspoint = 3;
587 data.nlist = 0;
588 data.mode = mode;
589 data.endlist = true;
590 data.fvNmin = 0;
591 data.fNmin = 0;
592 data.zMin = 0;
593 production3Sp(data);
594 }
597 data.i_seed_Pro = data.l_seeds_Pro.begin();
598
599
600 if (msgLvl(MSG::DEBUG)) {
601 data.nprint=1;
602 dump(data, msg(MSG::DEBUG));
603 }
604}
605
610
611void InDet::SiSpacePointsSeedMaker_ATLxk::find3Sp(const EventContext&, EventData& data, const std::list<Trk::Vertex>& lv, const double* ZVertex) const
612{
613 if (not data.initialized) initializeEventData(data);
614
617 data.zminU = ZVertex[0];
618 if (data.zminU < m_zmin) data.zminU = m_zmin;
619 data.zmaxU = ZVertex[1];
620 if (data.zmaxU > m_zmax) data.zmaxU = m_zmax;
621
623 int mode = 2;
624 if (lv.begin()!=lv.end()) mode = 3;
629 bool newv = newVertices(data, lv);
631 if (newv || !data.state || data.nspoint!=3 || data.mode!=mode || data.nlist) {
632 data.i_seede_Pro = data.l_seeds_Pro.begin();
633 data.state = 1;
634 data.nspoint = 3;
635 data.nlist = 0;
636 data.mode = mode;
637 data.endlist = true;
638 data.fvNmin = 0;
639 data.fNmin = 0;
640 data.zMin = 0;
641 production3Sp(data);
642 }
645 data.i_seed_Pro = data.l_seeds_Pro.begin();
646
647 if (msgLvl(MSG::DEBUG)) {
648 data.nprint=1;
649 dump(data, msg(MSG::DEBUG));
650 }
651}
652
658
659void InDet::SiSpacePointsSeedMaker_ATLxk::findVSp(const EventContext&, EventData& data, const std::list<Trk::Vertex>& lv) const
660{
661 if (not data.initialized) initializeEventData(data);
662
663 data.zminU = m_zmin;
664 data.zmaxU = m_zmax;
665
666 int mode = 5;
667 if (lv.begin()!=lv.end()) mode = 6;
668 bool newv = newVertices(data, lv);
669
670 if (newv || !data.state || data.nspoint!=4 || data.mode!=mode || data.nlist) {
671 data.i_seede_Pro = data.l_seeds_Pro.begin();
672 data.state = 1;
673 data.nspoint = 4;
674 data.nlist = 0;
675 data.mode = mode;
676 data.endlist = true;
677 data.fvNmin = 0;
678 data.fNmin = 0;
679 data.zMin = 0;
680 production3Sp(data);
681 }
682 data.i_seed_Pro = data.l_seeds_Pro.begin();
683
684 if (msgLvl(MSG::DEBUG)) {
685 data.nprint=1;
686 dump(data, msg(MSG::DEBUG));
687 }
688}
689
693
694MsgStream& InDet::SiSpacePointsSeedMaker_ATLxk::dump(EventData& data, MsgStream& out) const
695{
696 if (not data.initialized) initializeEventData(data);
697
698 if (data.nprint) return dumpEvent(data, out);
699 return dumpConditions(data, out);
700}
701
705
707{
708 int n = 42-m_spacepointsPixel.key().size();
709 std::string s2;
710 for (int i=0; i<n; ++i) s2.append(" ");
711 s2.append("|");
712 n = 42-m_spacepointsSCT.key().size();
713 std::string s3;
714 for (int i=0; i<n; ++i) s3.append(" ");
715 s3.append("|");
716 n = 42-m_spacepointsOverlap.key().size();
717 std::string s4;
718 for (int i=0; i<n; ++i) s4.append(" ");
719 s4.append("|");
720 n = 42-m_beamSpotKey.key().size();
721 std::string s5;
722 for (int i=0; i<n; ++i) s5.append(" ");
723 s5.append("|");
724
725 out<<"|---------------------------------------------------------------------|"
726 <<endmsg;
727 out<<"| Pixel space points | "<<m_spacepointsPixel.key() <<s2
728 <<endmsg;
729 out<<"| SCT space points | "<<m_spacepointsSCT.key()<<s3
730 <<endmsg;
731 out<<"| Overlap space points | "<<m_spacepointsOverlap.key()<<s4
732 <<endmsg;
733 out<<"| BeamConditionsService | "<<m_beamSpotKey.key()<<s5
734 <<endmsg;
735 out<<"| usePixel | "
736 <<std::setw(12)<<m_pixel
737 <<" |"<<endmsg;
738 out<<"| useSCT | "
739 <<std::setw(12)<<m_sct
740 <<" |"<<endmsg;
741 out<<"| maxSize | "
742 <<std::setw(12)<<m_maxsize
743 <<" |"<<endmsg;
744 out<<"| maxSizeSP | "
745 <<std::setw(12)<<m_maxsizeSP
746 <<" |"<<endmsg;
747 out<<"| pTmin (mev) | "
748 <<std::setw(12)<<std::setprecision(5)<<m_ptmin
749 <<" |"<<endmsg;
750 out<<"| max radius SP | "
751 <<std::setw(12)<<std::setprecision(5)<<m_r_rmax
752 <<" |"<<endmsg;
753 out<<"| radius step | "
754 <<std::setw(12)<<std::setprecision(5)<<m_binSizeR
755 <<" |"<<endmsg;
756 out<<"| min Z-vertex position | "
757 <<std::setw(12)<<std::setprecision(5)<<m_zmin
758 <<" |"<<endmsg;
759 out<<"| max Z-vertex position | "
760 <<std::setw(12)<<std::setprecision(5)<<m_zmax
761 <<" |"<<endmsg;
762 out<<"| min radius first SP(2) | "
763 <<std::setw(12)<<std::setprecision(5)<<m_r1minv
764 <<" |"<<endmsg;
765 out<<"| min radius second SP(2) | "
766 <<std::setw(12)<<std::setprecision(5)<<m_r2minv
767 <<" |"<<endmsg;
768 out<<"| max radius first SP(2) | "
769 <<std::setw(12)<<std::setprecision(5)<<m_r1maxv
770 <<" |"<<endmsg;
771 out<<"| max radius second SP(2) | "
772 <<std::setw(12)<<std::setprecision(5)<<m_r2maxv
773 <<" |"<<endmsg;
774 out<<"| min space points dR | "
775 <<std::setw(12)<<std::setprecision(5)<<m_drmin
776 <<" |"<<endmsg;
777 out<<"| max space points dR | "
778 <<std::setw(12)<<std::setprecision(5)<<m_drmax
779 <<" |"<<endmsg;
780 out<<"| max dZ impact | "
781 <<std::setw(12)<<std::setprecision(5)<<m_dzver
782 <<" |"<<endmsg;
783 out<<"| max dZ/dR impact | "
784 <<std::setw(12)<<std::setprecision(5)<<m_dzdrver
785 <<" |"<<endmsg;
786 out<<"| max impact | "
787 <<std::setw(12)<<std::setprecision(5)<<m_maxdImpact
788 <<" |"<<endmsg;
789 out<<"| max impact sss | "
790 <<std::setw(12)<<std::setprecision(5)<<m_maxdImpactSSS
791 <<" |"<<endmsg;
792 out<<"|---------------------------------------------------------------------|"
793 <<endmsg;
794 out<<"| Beam X center | "
795 <<std::setw(12)<<std::setprecision(5)<<data.xbeam[0]
796 <<" |"<<endmsg;
797 out<<"| Beam Y center | "
798 <<std::setw(12)<<std::setprecision(5)<<data.ybeam[0]
799 <<" |"<<endmsg;
800 out<<"| Beam Z center | "
801 <<std::setw(12)<<std::setprecision(5)<<data.zbeam[0]
802 <<" |"<<endmsg;
803 out<<"| Beam X-axis direction | "
804 <<std::setw(12)<<std::setprecision(5)<<data.xbeam[1]
805 <<std::setw(12)<<std::setprecision(5)<<data.xbeam[2]
806 <<std::setw(12)<<std::setprecision(5)<<data.xbeam[3]
807 <<" |"<<endmsg;
808 out<<"| Beam Y-axis direction | "
809 <<std::setw(12)<<std::setprecision(5)<<data.ybeam[1]
810 <<std::setw(12)<<std::setprecision(5)<<data.ybeam[2]
811 <<std::setw(12)<<std::setprecision(5)<<data.ybeam[3]
812 <<" |"<<endmsg;
813 out<<"| Beam Z-axis direction | "
814 <<std::setw(12)<<std::setprecision(5)<<data.zbeam[1]
815 <<std::setw(12)<<std::setprecision(5)<<data.zbeam[2]
816 <<std::setw(12)<<std::setprecision(5)<<data.zbeam[3]
817 <<" |"<<endmsg;
818 out<<"|---------------------------------------------------------------------|"
819 <<endmsg;
820 return out;
821}
822
826
828{
829 out<<"|---------------------------------------------------------------------|"
830 <<endmsg;
831 out<<"| ns | "
832 <<std::setw(12)<<data.ns
833 <<" |"<<endmsg;
834 out<<"| nsaz | "
835 <<std::setw(12)<<data.nsaz
836 <<" |"<<endmsg;
837 out<<"| nsazv | "
838 <<std::setw(12)<<data.nsazv
839 <<" |"<<endmsg;
840 out<<"| seeds | "
841 <<std::setw(12)<<data.l_seeds_Pro.size()
842 <<" |"<<endmsg;
843 out<<"|---------------------------------------------------------------------|"
844 <<endmsg;
845 return out;
846}
847
851
853{
854 if (data.endlist) return;
855
856 data.i_seede_Pro = data.l_seeds_Pro.begin();
857
858 if (data.mode==0 or data.mode==1) production2Sp(data);
859 else if (data.mode==2 or data.mode==3) production3Sp(data);
860 else if (data.mode==5 or data.mode==6) production3Sp(data);
861
862 data.i_seed_Pro = data.l_seeds_Pro.begin();
863 ++data.nlist;
864}
865
869
870bool InDet::SiSpacePointsSeedMaker_ATLxk::newVertices(EventData& data, const std::list<Trk::Vertex>& lV) const
871{
872 unsigned int s1 = data.l_vertex.size();
873 unsigned int s2 = lV.size();
874
876 data.isvertex = false;
879 if (s1==0 && s2==0) return false;
880
882 data.l_vertex.clear();
884 if (s2 == 0) return false;
885
887 data.isvertex = true;
888 for (const Trk::Vertex& v: lV) {
889 data.l_vertex.insert(static_cast<float>(v.position().z()));
890 }
891
894 data.zminU = (*data.l_vertex. begin())-20.f;
895 if (data.zminU < m_zmin) data.zminU = m_zmin;
896 data.zmaxU = (*data.l_vertex.rbegin())+20.f;
897 if (data.zmaxU > m_zmax) data.zmaxU = m_zmax;
898
901 return false;
902}
903
905// Initiate frame work for seed generator
907
909{
910 // Tell clang to optimize assuming that FP operations may trap.
912
913 m_ptmin = std::max( std::abs(m_ptmin), float(100.*m_fieldScale));
917
919 if (std::abs(m_etamin) < .1) m_etamin = -m_etamax;
921 m_dzdrmax0 = 1.f/std::tan(2.f*std::atan(std::exp(-m_etamax)));
922 m_dzdrmin0 = 1.f/std::tan(2.f*std::atan(std::exp(-m_etamin)));
923
925 m_ipt = 1.f/std::abs(.9f*m_ptmin);
927
936
938 m_nBinsR = static_cast<int>((m_r_rmax+.1f)/m_binSizeR);
939
944
946 constexpr float twoPi = 2.f*M_PI;
947
949 const int nPhiBinsMax = arraySizePhi;
950 const float inverseSizePhiMax = static_cast<float>(nPhiBinsMax)/twoPi;
951
959 constexpr float inverseSizePhiMin = 100./60.;
960
966
969 constexpr float radiusPixelStart = 33.;
970 constexpr float radiusPixelEnd = 150.;
974 const float binSizePhi_PPP = m_pixel ? azimuthalStep(m_ptmin/m_fieldScale,m_maxdImpact,radiusPixelStart,radiusPixelEnd)/3.f : 0.f;
976 constexpr float radiusSctStart = 295.; ;
977 constexpr float radiusSctEnd = 560.;
978 const float binSizePhi_SSS = m_sct ? azimuthalStep(m_ptmin/m_fieldScale,m_maxdImpactSSS,radiusSctStart,radiusSctEnd)/3.f : 0.f;
980 m_inverseBinSizePhi = 1.f/std::max(binSizePhi_PPP, binSizePhi_SSS);
981 }
982 else {
985 float ptm = 400.;
987 if (m_ptmin/m_fieldScale < ptm) ptm = m_ptmin/m_fieldScale;
988 m_inverseBinSizePhi = ptm /60.f;
989 }
990
992 if (m_inverseBinSizePhi > inverseSizePhiMax) m_inverseBinSizePhi = inverseSizePhiMax;
993 else if (m_inverseBinSizePhi < inverseSizePhiMin) m_inverseBinSizePhi = inverseSizePhiMin;
994
996 m_maxPhiBin = static_cast<int>(twoPi*m_inverseBinSizePhi);
998 if (m_maxPhiBin >=nPhiBinsMax) m_maxPhiBin = nPhiBinsMax-1;
999
1002 const int nPhiBinsVertexMax = arraySizePhiV;
1003 const float inverseBinSizePhiVertexMax = static_cast<float>(nPhiBinsVertexMax)/twoPi;
1005 if (m_inverseBinSizePhiVertex > inverseBinSizePhiVertexMax) m_inverseBinSizePhiVertex = inverseBinSizePhiVertexMax;
1006 m_maxBinPhiVertex = static_cast<int>(twoPi*m_inverseBinSizePhiVertex);
1007 if (m_maxBinPhiVertex>=nPhiBinsVertexMax) m_maxBinPhiVertex = nPhiBinsVertexMax-1;
1008
1009
1013
1014 for (int phiBin=0; phiBin<=m_maxPhiBin; ++phiBin) {
1015
1016 int phiBelow = phiBin-1;
1017 if (phiBelow<0) phiBelow=m_maxPhiBin;
1018
1019 int phiAbove = phiBin+1;
1020 if (phiAbove>m_maxPhiBin) phiAbove=0;
1021
1023 for (int z=0; z<arraySizeZ; ++z) {
1024
1027
1028 int twoDbinSamePhi = phiBin * arraySizeZ+z;
1029 int twoDbinLowerPhi = phiBelow*arraySizeZ+z;
1030 int twoDbinHigherPhi = phiAbove* arraySizeZ+z;
1031
1032 m_nNeighbourCellsBottom [twoDbinSamePhi] = 3;
1033 m_nNeighbourCellsTop [twoDbinSamePhi] = 3;
1034
1035 m_neighbourCellsBottom[twoDbinSamePhi][0] = twoDbinSamePhi;
1036 m_neighbourCellsTop[twoDbinSamePhi][0] = twoDbinSamePhi;
1037
1038 m_neighbourCellsBottom[twoDbinSamePhi][1] = twoDbinLowerPhi;
1039 m_neighbourCellsTop[twoDbinSamePhi][1] = twoDbinLowerPhi;
1040
1041 m_neighbourCellsBottom[twoDbinSamePhi][2] = twoDbinHigherPhi;
1042 m_neighbourCellsTop[twoDbinSamePhi][2] = twoDbinHigherPhi;
1043
1052 if (z==5) {
1053 m_nNeighbourCellsTop [twoDbinSamePhi] = 9;
1054 // in the central z region, we include the two neighbouring
1055 // z slices for the top neighbour search
1056
1057 m_neighbourCellsTop[twoDbinSamePhi][3] = twoDbinSamePhi+1;
1058 m_neighbourCellsTop[twoDbinSamePhi][4] = twoDbinLowerPhi+1;
1059 m_neighbourCellsTop[twoDbinSamePhi][5] = twoDbinHigherPhi+1;
1060 m_neighbourCellsTop[twoDbinSamePhi][6] = twoDbinSamePhi-1;
1061 m_neighbourCellsTop[twoDbinSamePhi][7] = twoDbinLowerPhi-1;
1062 m_neighbourCellsTop[twoDbinSamePhi][8] = twoDbinHigherPhi-1;
1063 }
1064 // z > 5: positive z values, |z| > 250mm
1065 else if (z> 5) {
1066 // for the bottom SP search in positive non-central z, we include the
1067 // neighbouring Z region on the left (towards the IP) in the bottom
1068 // neighbour search
1069 m_nNeighbourCellsBottom [twoDbinSamePhi] = 6;
1070 m_neighbourCellsBottom[twoDbinSamePhi][3] = twoDbinSamePhi-1;
1071 m_neighbourCellsBottom[twoDbinSamePhi][4] = twoDbinLowerPhi-1;
1072 m_neighbourCellsBottom[twoDbinSamePhi][5] = twoDbinHigherPhi-1;
1073
1074 if (z<10) {
1080 m_nNeighbourCellsTop [twoDbinSamePhi] = 6;
1081 m_neighbourCellsTop[twoDbinSamePhi][3] = twoDbinSamePhi+1;
1082 m_neighbourCellsTop[twoDbinSamePhi][4] = twoDbinLowerPhi+1;
1083 m_neighbourCellsTop[twoDbinSamePhi][5] = twoDbinHigherPhi+1;
1084 }
1085 }
1086 // z < 5: negative z values, |z| > 250mm
1087 else {
1092 m_nNeighbourCellsBottom [twoDbinSamePhi] = 6;
1093 m_neighbourCellsBottom[twoDbinSamePhi][3] = twoDbinSamePhi+1;
1094 m_neighbourCellsBottom[twoDbinSamePhi][4] = twoDbinLowerPhi+1;
1095 m_neighbourCellsBottom[twoDbinSamePhi][5] = twoDbinHigherPhi+1;
1096
1097 if (z>0) {
1098 // if there is a z region on the left (away from the IP), we include it in the top
1099 // neighbour search
1100 m_nNeighbourCellsTop [twoDbinSamePhi] = 6;
1101 m_neighbourCellsTop[twoDbinSamePhi][3] = twoDbinSamePhi-1;
1102 m_neighbourCellsTop[twoDbinSamePhi][4] = twoDbinLowerPhi-1;
1103 m_neighbourCellsTop[twoDbinSamePhi][5] = twoDbinHigherPhi-1;
1104 }
1105 }
1106
1112 if (z==3) {
1113 m_nNeighbourCellsBottom[twoDbinSamePhi] = 9;
1114 m_neighbourCellsBottom[twoDbinSamePhi][6] = twoDbinSamePhi+2;
1115 m_neighbourCellsBottom[twoDbinSamePhi][7] = twoDbinLowerPhi+2;
1116 m_neighbourCellsBottom[twoDbinSamePhi][8] = twoDbinHigherPhi+2;
1117 } else if (z==7) {
1118 m_nNeighbourCellsBottom[twoDbinSamePhi] = 9;
1119 m_neighbourCellsBottom[twoDbinSamePhi][6] = twoDbinSamePhi-2;
1120 m_neighbourCellsBottom[twoDbinSamePhi][7] = twoDbinLowerPhi-2;
1121 m_neighbourCellsBottom[twoDbinSamePhi][8] = twoDbinHigherPhi-2;
1122 }
1123 }
1124 }
1125
1129 for (int phiBin=0; phiBin<=m_maxBinPhiVertex; ++phiBin) {
1130
1131 int phiBinBelow = phiBin-1;
1132 if (phiBinBelow<0) phiBinBelow=m_maxBinPhiVertex;
1133
1134 int phiBinTop = phiBin+1;
1135 if (phiBinTop>m_maxBinPhiVertex) phiBinTop=0;
1136
1138 for (int zbin=0; zbin<arraySizeZV; ++zbin) {
1139
1140 int twoDbinSamePhi = phiBin *arraySizeZV+zbin;
1141 int twoDbinLowerPhi = phiBinBelow*arraySizeZV+zbin;
1142 int twoDbinHigherPhi = phiBinTop*arraySizeZV+zbin;
1143
1145 m_nNeighboursVertexPhiZ[twoDbinSamePhi] = 3;
1146 m_neighboursVertexPhiZ[twoDbinSamePhi][0] = twoDbinSamePhi;
1147 m_neighboursVertexPhiZ[twoDbinSamePhi][1] = twoDbinLowerPhi;
1148 m_neighboursVertexPhiZ[twoDbinSamePhi][2] = twoDbinHigherPhi;
1149
1151 if (zbin>1) {
1152 m_nNeighboursVertexPhiZ[twoDbinSamePhi] = 6;
1153 m_neighboursVertexPhiZ[twoDbinSamePhi][3] = twoDbinSamePhi-1;
1154 m_neighboursVertexPhiZ[twoDbinSamePhi][4] = twoDbinLowerPhi-1;
1155 m_neighboursVertexPhiZ[twoDbinSamePhi][5] = twoDbinHigherPhi-1;
1156 }
1158 else if (zbin<1) {
1159 m_nNeighboursVertexPhiZ[twoDbinSamePhi] = 6;
1160 m_neighboursVertexPhiZ[twoDbinSamePhi][3] = twoDbinSamePhi+1;
1161 m_neighboursVertexPhiZ[twoDbinSamePhi][4] = twoDbinLowerPhi+1;
1162 m_neighboursVertexPhiZ[twoDbinSamePhi][5] = twoDbinHigherPhi+1;
1163 }
1164 }
1165 }
1166}
1167
1171
1173{
1175
1176 const Amg::Vector3D &bsCentre = beamSpotHandle->beamPos();
1177 double tx = std::tan(beamSpotHandle->beamTilt(0));
1178 double ty = std::tan(beamSpotHandle->beamTilt(1));
1179
1180 double phi = std::atan2(ty,tx);
1181 double theta = std::acos(1./std::sqrt(1.+tx*tx+ty*ty));
1182 double sinTheta = std::sin(theta);
1183 double cosTheta = std::cos(theta);
1184 double sinPhi = std::sin(phi);
1185 double cosPhi = std::cos(phi);
1186
1187 data.xbeam[0] = static_cast<float>(bsCentre.x());
1188 data.xbeam[1] = static_cast<float>(cosTheta*cosPhi*cosPhi+sinPhi*sinPhi);
1189 data.xbeam[2] = static_cast<float>(cosTheta*sinPhi*cosPhi-sinPhi*cosPhi);
1190 data.xbeam[3] =-static_cast<float>(sinTheta*cosPhi );
1191
1192 data.ybeam[0] = static_cast<float>(bsCentre.y());
1193 data.ybeam[1] = static_cast<float>(cosTheta*cosPhi*sinPhi-sinPhi*cosPhi);
1194 data.ybeam[2] = static_cast<float>(cosTheta*sinPhi*sinPhi+cosPhi*cosPhi);
1195 data.ybeam[3] =-static_cast<float>(sinTheta*sinPhi );
1196
1197 data.zbeam[0] = static_cast<float>(bsCentre.z());
1198 data.zbeam[1] = static_cast<float>(sinTheta*cosPhi);
1199 data.zbeam[2] = static_cast<float>(sinTheta*sinPhi);
1200 data.zbeam[3] = static_cast<float>(cosTheta);
1201}
1202
1206
1208(EventData& data, const Trk::SpacePoint*const& sp, std::array<float,3> & r)
1209{
1210 r[0] = static_cast<float>(sp->globalPosition().x())-data.xbeam[0];
1211 r[1] = static_cast<float>(sp->globalPosition().y())-data.ybeam[0];
1212 r[2] = static_cast<float>(sp->globalPosition().z())-data.zbeam[0];
1213}
1214
1218
1220{
1221 constexpr float twoPi = 2.*M_PI;
1222
1223 int firstRadialBin = 0;
1224 bool ibl = false;
1225
1238
1239 const std::map<float, int> ztoBin{
1240 {-2500. ,0},
1241 {-1400. ,1},
1242 {-925. ,2},
1243 {-450. ,3},
1244 {-250 ,4},
1245 { 250 ,5},
1246 { 450 ,6},
1247 { 925 ,7},
1248 { 1400 ,8},
1249 { 2500 ,9},
1250 { 100000 ,10},
1251 };
1252
1253 for (int radialBin=data.r_first; radialBin<m_nBinsR; ++radialBin) {
1254
1256 if (!data.r_map[radialBin]) continue;
1258 if (firstRadialBin == 0) firstRadialBin = radialBin;
1259
1261 if (data.iteration) {
1268 if (!data.r_Sorted[radialBin].front()->spacepoint->clusterList().second) {
1269 if (radialBin * m_binSizeR < m_radiusCutIBL) ibl = true;
1270 }
1271
1273
1275 else if (ibl) {
1276 break;
1277 }
1283 else if (radialBin > 175) {
1284 break;
1285 }
1286 }
1287 // loop over the space points in the r-bin and sort them into the 2d phi-z binning
1288 for (InDet::SiSpacePointForSeed* SP : data.r_Sorted[radialBin]) {
1289
1292 float Phi = SP->phi();
1293 if (Phi<0.) Phi+=twoPi; // phi is defined in [0..2pi] for the binning
1294 int phiBin = static_cast<int>(Phi*m_inverseBinSizePhi);
1296 if (phiBin < 0) {
1297 phiBin = m_maxPhiBin;
1298 } else if (phiBin > m_maxPhiBin) {
1299 phiBin = 0;
1300 }
1301
1302 float Z = SP->z();
1309 int zBin{0};
1310 auto bound = ztoBin.lower_bound(Z);
1312 if (bound == ztoBin.end()){
1313 --bound;
1314 }
1315 zBin=bound->second;
1316
1318 int twoDbin = phiBin*arraySizeZ+zBin;
1321 ++data.nsaz;
1322 // push our space point into the 2D binned array
1323 data.rfz_Sorted[twoDbin].push_back(SP);
1327 if (!data.rfz_map[twoDbin]++) data.rfz_index[data.nrfz++] = twoDbin;
1328 }
1329 }
1332 if (!m_sct && firstRadialBin && static_cast<float>(firstRadialBin)*m_binSizeR < m_radiusCutIBL) {
1334 }
1335
1336 data.state = 0;
1337}
1338
1339
1340
1341float InDet::SiSpacePointsSeedMaker_ATLxk::azimuthalStep(const float pTmin,const float maxd0,const float Rmin,const float Rmax)
1342{
1343 // Tell clang to optimize assuming that FP exceptions can trap.
1344 // Otherwise, it can vectorize the division, which can lead to
1345 // spurious division-by-zero traps from unused vector lanes.
1347
1351 float Rm = pTmin/.6f;
1352
1359 float worstCaseD0 = maxd0;
1360 if (maxd0 > Rmin) worstCaseD0 = Rmin;
1361
1362 float sI = std::abs(std::asin(worstCaseD0/Rmin) - std::asin(worstCaseD0/Rmax));
1363 float sF = std::abs(std::asin(std::min(1.f,Rmax/(2.f*Rm))) -
1364 std::asin(std::min(1.f,Rmin/(2.f*Rm))));
1365 return sI+sF;
1366}
1367
1368
1369
1371// Erase space point information
1373
1375{
1376 for (int i=0; i<data.nrfz; ++i) {
1377 int n = data.rfz_index[i];
1378 data.rfz_map[n] = 0;
1379 data.rfz_Sorted[n].clear();
1380 }
1381
1382 for (int i=0; i<data.nrfzv; ++i) {
1383 int n = data.rfzv_index[i];
1384 data.rfzv_map[n] = 0;
1385 data.rfzv_Sorted[n].clear();
1386 }
1387 data.state = 0;
1388 data.nsaz = 0;
1389 data.nsazv = 0;
1390 data.nrfz = 0;
1391 data.nrfzv = 0;
1392}
1393
1395// Test is space point used
1397
1398
1399
1401// 2 space points seeds production
1403
1405{
1406 if (data.nsazv<2) return;
1407
1408 std::vector<InDet::SiSpacePointForSeed*>::iterator r0,r0e,r,re;
1409 int nseed = 0;
1410
1411 // Loop thorugh all azimuthal regions
1412 //
1413 for (int f=data.fvNmin; f<=m_maxBinPhiVertex; ++f) {
1414
1415 // For each azimuthal region loop through Z regions
1416 //
1417 int z = 0;
1418 if (!data.endlist) z = data.zMin;
1419 for (; z<arraySizeZV; ++z) {
1420
1421 int a = f*arraySizeZV+z;
1422 if (!data.rfzv_map[a]) continue;
1423 r0 = data.rfzv_Sorted[a].begin();
1424 r0e = data.rfzv_Sorted[a].end ();
1425
1426 if (!data.endlist) {
1427 r0 = data.rMin;
1428 data.endlist = true;
1429 }
1430
1431 // Loop through trigger space points
1432 //
1433 for (; r0!=r0e; ++r0) {
1434
1435 float X = (*r0)->x();
1436 float Y = (*r0)->y();
1437 float R = (*r0)->radius();
1438 if (R<m_r2minv) continue;
1439 if (R>m_r2maxv) break;
1440 float Z = (*r0)->z();
1441 float ax = X/R;
1442 float ay = Y/R;
1443
1444 // Bottom links production
1445 //
1446 int numberBottomCells = m_nNeighboursVertexPhiZ[a];
1447 for (int i=0; i<numberBottomCells; ++i) {
1448
1449 int an = m_neighboursVertexPhiZ[a][i];
1450 if (!data.rfzv_map[an]) continue;
1451
1452 r = data.rfzv_Sorted[an].begin();
1453 re = data.rfzv_Sorted[an].end();
1454
1455 for (; r!=re; ++r) {
1456
1457 float Rb =(*r)->radius();
1458 if (Rb<m_r1minv) continue;
1459 if (Rb>m_r1maxv) break;
1460 float dR = R-Rb;
1461 if (dR<m_drminv) break;
1462 if (dR>m_drmax) continue;
1463 float dZ = Z-(*r)->z();
1464 float Tz = dZ/dR;
1465 if (Tz<data.dzdrmin || Tz>data.dzdrmax) continue;
1466 float Zo = Z-R*Tz;
1467
1468 // Comparison with vertices Z coordinates
1469 //
1470 if (!isZCompatible(data, Zo, Rb, Tz)) continue;
1471
1472 // Momentum cut
1473 //
1474 float dx =(*r)->x()-X;
1475 float dy =(*r)->y()-Y;
1476 float x = dx*ax+dy*ay;
1477 float y =-dx*ay+dy*ax;
1478 float xy = x*x+y*y;
1479 if (xy == 0.) continue;
1480 float r2 = 1.f/xy;
1481 float Ut = x*r2;
1482 float Vt = y*r2;
1483 float UR = Ut*R+1.f;
1484 if (UR == 0.) continue;
1485 float A = Vt*R/UR;
1486 float B = Vt-A*Ut;
1487 if (std::abs(B*data.K) > m_ipt*std::sqrt(1.f+A*A)) continue;
1488 ++nseed;
1489 newSeed(data, (*r), (*r0), Zo);
1490 }
1491 }
1492 if (nseed < m_maxsize) continue;
1493 data.endlist=false;
1494 data.rMin = (++r0);
1495 data.fvNmin=f;
1496 data.zMin=z;
1497 return;
1498 }
1499 }
1500 }
1501 data.endlist = true;
1502}
1503
1505// Production 3 space points seeds
1507
1509{
1510
1512 if (data.nsaz<3) return;
1513
1524
1525
1535 const std::array<int,arraySizeZ> zBinIndex {5,6,7,8,9,10,4,3,2,1,0};
1536
1539 std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator,arraySizeNeighbourBins> iter_topCands;
1540 std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator,arraySizeNeighbourBins> iter_endTopCands;
1541 std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator,arraySizeNeighbourBins> iter_bottomCands;
1542 std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator,arraySizeNeighbourBins> iter_endBottomCands;
1543
1545 int nseed = 0;
1547 data.endlist = true;
1548
1550 for (int phiBin=data.fNmin; phiBin<=m_maxPhiBin; ++phiBin) {
1551
1553 int z = 0;
1555 if (!data.endlist) z = data.zMin;
1556
1560 for (; z<arraySizeZ; ++z) {
1561
1562 int phiZbin = phiBin *arraySizeZ+zBinIndex[z];
1563
1565 if (!data.rfz_map[phiZbin]) continue;
1566
1569 int numberBottomCells = 0;
1570 int numberTopCells = 0;
1571
1576 for (int neighbourCellNumber=0; neighbourCellNumber<m_nNeighbourCellsBottom[phiZbin]; ++neighbourCellNumber) {
1577
1578 int theNeighbourCell = m_neighbourCellsBottom[phiZbin][neighbourCellNumber];
1580 if (!data.rfz_map[theNeighbourCell]) continue;
1582 //cppcheck-suppress objectIndex; false positive
1583 iter_bottomCands [numberBottomCells] = data.rfz_Sorted[theNeighbourCell].begin();
1584 //cppcheck-suppress objectIndex; false positive
1585 iter_endBottomCands[numberBottomCells++] = data.rfz_Sorted[theNeighbourCell].end();
1586 }
1587
1592 for (int neighbourCellNumber=0; neighbourCellNumber<m_nNeighbourCellsTop[phiZbin]; ++neighbourCellNumber) {
1593
1594 int theNeighbourCell = m_neighbourCellsTop[phiZbin][neighbourCellNumber];
1596 if (!data.rfz_map[theNeighbourCell]) continue;
1598 //cppcheck-suppress objectIndex; false positive
1599 iter_topCands [numberTopCells] = data.rfz_Sorted[theNeighbourCell].begin();
1600 //cppcheck-suppress objectIndex; false positive
1601 iter_endTopCands[numberTopCells++] = data.rfz_Sorted[theNeighbourCell].end();
1602 }
1603
1605 if (!data.trigger) production3Sp (data, iter_bottomCands, iter_endBottomCands, iter_topCands, iter_endTopCands, numberBottomCells, numberTopCells, nseed,zBinIndex[z]);
1606 else production3SpTrigger(data, iter_bottomCands, iter_endBottomCands, iter_topCands, iter_endTopCands, numberBottomCells, numberTopCells, nseed);
1607 }
1608
1615 if (nseed>=m_maxsize) {
1616 data.endlist=false;
1617 data.fNmin = phiBin+1;
1618 return;
1619 }
1620 }
1622 data.endlist = true;
1623}
1624
1628
1630(EventData& data,
1631 std::array <std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & iter_bottomCands ,
1632 std::array <std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & iter_endBottomCands,
1633 std::array <std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & iter_topCands ,
1634 std::array <std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & iter_endTopCands,
1635 const int numberBottomCells, const int numberTopCells, int& nseed, const int zbin) const
1636{
1642
1644 std::vector<InDet::SiSpacePointForSeed*>::iterator iter_centralSP=iter_bottomCands[0];
1645 std::vector<InDet::SiSpacePointForSeed*>::iterator iter_otherSP;
1646
1652
1654 bool isStrip = ((*iter_centralSP)->spacepoint->clusterList().second);
1655
1658 bool isBarrelRegion = (zbin >=4 && zbin <= 6);
1659 bool isTransitionRegion = (zbin == 3 || zbin == 7); // region in z = 450 - 925mm
1660 float rmin = 40.;
1661 float rmax = 140.;
1662
1663 if (isBarrelRegion){
1664 rmax = 100;
1665 }
1667 if(isStrip) {
1669 rmin = 285.;
1670 rmax = 450.;
1671 if (isBarrelRegion){
1672 rmin = 335;
1673 rmax = 450.;
1674 }
1675 else if (isTransitionRegion){
1676 rmin = 335;
1677 rmax = 550;
1679 }
1680 }
1681
1683 for(; iter_centralSP!=iter_endBottomCands[0]; ++iter_centralSP) {
1684 if((*iter_centralSP)->radius() > rmin) break;
1685 }
1686
1689 iter_topCands[0] = iter_centralSP;
1690 ++iter_topCands[0];
1691
1693 const float ipt2K = data.ipt2K;
1694 const float ipt2C = data.ipt2C;
1695 const float COFK = data.COFK;
1696 const float maxd0cut = m_maxdImpact;
1697 const float maxd0cutstrips = m_maxdImpactDecays;
1698 const float zmin = data.zminU;
1699 const float zmax = data.zmaxU;
1700 const float dzdrmax = data.dzdrmax;
1701 const float dzdrmin = data.dzdrmin;
1702 data.CmSp.clear();
1703
1706 size_t SPcapacity = data.SP.size();
1707
1709 for (; iter_centralSP!=iter_endBottomCands[0]; ++iter_centralSP) {
1710
1711 const float& R = (*iter_centralSP)->radius();
1712 if(R > rmax) break;
1713
1715 const float& X = (*iter_centralSP)->x();
1716 const float& Y = (*iter_centralSP)->y();
1717 const float& Z = (*iter_centralSP)->z();
1718
1721 double absZ = std::abs(Z);
1723 if (isStrip && absZ > 2650. ) continue;
1725 if (!isStrip && absZ > 600.) continue;
1727 if (isStrip && isTransitionRegion && absZ < 750. && R > 450.) continue;
1728
1731 size_t Nb = 0;
1732
1735 for (int cell=0; cell<numberBottomCells; ++cell) {
1737 for (iter_otherSP=iter_bottomCands[cell]; iter_otherSP!=iter_endBottomCands[cell]; ++iter_otherSP) {
1738
1740 const float& Rb =(*iter_otherSP)->radius();
1741 float dR = R-Rb;
1742
1745 if (dR > m_drmax) {
1746 iter_bottomCands[cell]=iter_otherSP;
1747 continue;
1748 }
1751 if (dR < m_drmin || (data.iteration && (*iter_otherSP)->spacepoint->clusterList().second)) break;
1752
1754 const float dZdR = (Z-(*iter_otherSP)->z())/dR;
1756 const float absdZdR = std::abs(dZdR);
1758 if (absdZdR < dzdrmin or absdZdR > dzdrmax) continue;
1759
1762 const float z0 = Z-R*dZdR;
1763 if(z0 > zmax || z0 < zmin) continue;
1765 data.SP[Nb] = (*iter_otherSP);
1766 if(m_writeNtuple) data.SP[Nb]->setDZDR(dZdR);
1770 if (++Nb==SPcapacity){
1771 data.resizeSPCont();
1772 SPcapacity=data.SP.size();
1773 }
1774 }
1775 }
1776
1778 if (!Nb) continue;
1779
1783 size_t Nt = Nb;
1784
1786
1788 for (int cell=0; cell<numberTopCells; ++cell) {
1790 for (iter_otherSP=iter_topCands[cell];iter_otherSP!=iter_endTopCands[cell]; ++iter_otherSP) {
1791
1793 float Rt =(*iter_otherSP)->radius();
1794 float dR = Rt-R;
1795
1797 if (dR<m_drmin) {
1798 iter_topCands[cell]=iter_otherSP;
1799 continue;
1800 }
1802 if (dR>m_drmax) break;
1803
1805 float dZdR = ((*iter_otherSP)->z()-Z)/dR;
1806 float absdZdR = std::abs(dZdR);
1807 if (absdZdR < dzdrmin or absdZdR > dzdrmax) continue;
1808
1811 float z0 = Z-R*dZdR;
1812 if(z0 > zmax || z0 < zmin) continue;
1814 data.SP[Nt] = (*iter_otherSP);
1815 if (m_writeNtuple) data.SP[Nt]->setDZDR(dZdR);
1819 if (++Nt==SPcapacity) {
1820 data.resizeSPCont();
1821 SPcapacity=data.SP.size();
1822 }
1823 }
1824 }
1825
1827 if (!(Nt-Nb)) continue;
1828
1830 float covr0 = (*iter_centralSP)->covr ();
1831 float covz0 = (*iter_centralSP)->covz ();
1832
1834 float ax = X/R;
1835 float ay = Y/R;
1836
1839 for (size_t i=0; i<Nt; ++i) {
1840
1841 InDet::SiSpacePointForSeed* sp = data.SP[i];
1842
1845 float dx = sp->x()-X;
1846 float dy = sp->y()-Y;
1847 float dz = sp->z()-Z;
1848 float x = dx*ax+dy*ay;
1849 float y = dy*ax-dx*ay;
1850
1852 float r2 = 1.f/(x*x+y*y);
1854 float dr = std::sqrt(r2);
1857 float tz = dz*dr;
1858
1861 if (i < Nb) tz = -tz;
1862
1864 data.Tz[i] = tz;
1865 data.Zo[i] = Z-R*tz;
1866 data.R [i] = dr;
1867 data.U [i] = x*r2;
1868 data.V [i] = y*r2;
1869 data.Er[i] = ((covz0+sp->covz())+(tz*tz)*(covr0+sp->covr()))*r2;
1870 }
1871
1872 data.nOneSeeds = 0;
1873 data.mapOneSeeds_Pro.clear();
1874
1877 for (size_t b=0; b<Nb; ++b) {
1878
1880 float Zob = data.Zo[b];
1881 float Tzb = data.Tz[b];
1882 float Erb = data.Er[b];
1883 float Vb = data.V [b];
1884 float Ub = data.U [b];
1885 float Tzb2 = (1.f+Tzb*Tzb);
1886 float sTzb2 = std::sqrt(Tzb2);
1887 float sigmaSquaredScatteringPtDependent = Tzb2*COFK;
1888 float sigmaSquaredScatteringMinPt = Tzb2*ipt2C;
1891 float d0max = maxd0cut;
1893 if (data.SP[b]->spacepoint->clusterList().second) d0max = maxd0cutstrips;
1894
1896 for (size_t t=Nb; t<Nt; ++t) {
1897
1902
1904 float meanOneOverTanTheta = (Tzb+data.Tz[t])/2.f;
1905 float theta = 0.;
1906 if(m_writeNtuple){
1908 theta = std::atan(1.f/meanOneOverTanTheta);
1909 }
1911 float sigmaSquaredSpacePointErrors = Erb+data.Er[t]
1912 + 2.f * covz0 * data.R[t]*data.R[b]
1913 + 2.f * covr0 * data.R[t]*data.R[b] * meanOneOverTanTheta * meanOneOverTanTheta; // mixed term with r-uncertainy on central SP
1915 float remainingSquaredDelta = (Tzb-data.Tz[t])*(Tzb-data.Tz[t]) - sigmaSquaredSpacePointErrors;
1916
1919 if (remainingSquaredDelta - sigmaSquaredScatteringMinPt > 0 ) continue;
1920
1942 float deltaU = data.U[t]-Ub;
1943 if (deltaU == 0.) continue;
1944 float A = (data.V[t]-Vb)/deltaU;
1945 float B = Vb-A*Ub;
1946 float onePlusAsquare = 1.f+A*A;
1947 float BSquare = B*B;
1948
1961 if (BSquare > ipt2K*onePlusAsquare || remainingSquaredDelta*onePlusAsquare > BSquare*sigmaSquaredScatteringPtDependent) continue;
1978 float d0 = 0;
1979 if(std::abs(B) < 1e-10) d0 = std::abs((A-B*R)*R);
1981 else{
1982 float x0 = -A/(2.f*B);
1983 float rTrack = std::sqrt(onePlusAsquare/BSquare)*.5f;
1984 d0 = std::abs(-rTrack + std::sqrt(rTrack*rTrack +2.f*x0*R +R*R));
1985 }
1986
1988 if (d0 <= d0max) {
1990 float dr = data.R[b];
1991 if (data.R[t] < data.R[b]) dr = data.R[t];
1995 data.SP[t]->setScorePenalty(std::abs((Tzb-data.Tz[t])/(dr*sTzb2)));
1996 data.SP[t]->setParam(d0);
1997
1998 if(m_writeNtuple){
2000 data.SP[t]->setEta(-std::log(std::tan(0.5f*theta)));
2001 data.SP[t]->setPt(std::sqrt(onePlusAsquare/BSquare)/(1000.f*data.K));
2002 }
2004 data.CmSp.emplace_back(B/std::sqrt(onePlusAsquare), data.SP[t]);
2006
2007 }
2008
2009 }
2011 if (!data.CmSp.empty()) {
2012 newOneSeedWithCurvaturesComparison(data, data.SP[b], (*iter_centralSP), Zob);
2013 }
2014 }
2016 fillSeeds(data);
2017 nseed += data.fillOneSeeds;
2018
2019 }
2020}
2021
2023// Production 3 space points seeds in ROI
2025
2027(EventData& data,
2028 std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & rb ,
2029 std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & rbe,
2030 std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & rt ,
2031 std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & rte,
2032 const int numberBottomCells, const int numberTopCells, int& nseed) const
2033{
2034 constexpr float twoPi = 2.*M_PI;
2035
2036 std::vector<InDet::SiSpacePointForSeed*>::iterator r0=rb[0],r;
2037
2038 float rmin = 40.;
2039 float rmax = 140.;
2040 if((*r0)->spacepoint->clusterList().second) {
2041 rmin = 280.;
2042 rmax = 540.;
2043 }
2044
2045 for(; r0!=rbe[0]; ++r0) {if((*r0)->radius() > rmin) break;}
2046 rt[0] = r0; ++rt[0];
2047
2048 float ipt2K = data.ipt2K;
2049 float ipt2C = data.ipt2C;
2050 float COFK = data.COFK;
2051 float maxd0cut = m_maxdImpact;
2052 float maxd0cutstrips = m_maxdImpactSSS;
2053
2054 data.CmSp.clear();
2055
2056 // Loop through all trigger space points
2057 //
2058 for (; r0!=rbe[0]; ++r0) {
2059
2060 data.nOneSeeds = 0;
2061 data.mapOneSeeds_Pro.clear();
2062
2063 float R = (*r0)->radius();
2064 if(R>rmax) break;
2065
2066 const Trk::Surface* sur0 = (*r0)->sur();
2067 float X = (*r0)->x();
2068 float Y = (*r0)->y();
2069 float Z = (*r0)->z();
2070 int Nb = 0;
2071
2072 // Bottom links production
2073 //
2074 for (int i=0; i<numberBottomCells; ++i) {
2075
2076 for (r=rb[i]; r!=rbe[i]; ++r) {
2077
2078 float Rb =(*r)->radius();
2079
2080 float dR = R-Rb;
2081 if (dR > m_drmax) {
2082 rb[i]=r;
2083 continue;
2084 }
2085 if ((*r)->sur()==sur0) continue;
2086
2087 if (dR < m_drmin || (data.iteration && (*r)->spacepoint->clusterList().second)) break;
2088
2089 // Comparison with bottom and top Z
2090 //
2091 float Tz = (Z-(*r)->z())/dR;
2092 float Zo = Z-R*Tz;
2093 if (Zo < data.zminB || Zo > data.zmaxB) continue;
2094 float Zu = Z+(550.f-R)*Tz;
2095 if (Zu < data.zminU || Zu > data.zmaxU) continue;
2096 data.SP[Nb] = (*r);
2097 if (++Nb==m_maxsizeSP) goto breakb;
2098 }
2099 }
2100 breakb:
2101 if (!Nb || Nb==m_maxsizeSP) continue;
2102 int Nt = Nb;
2103
2104 // Top links production
2105 //
2106 for (int i=0; i<numberTopCells; ++i) {
2107
2108 for (r=rt[i]; r!=rte[i]; ++r) {
2109
2110 float Rt =(*r)->radius();
2111 float dR = Rt-R;
2112
2113 if (dR<m_drmin) {
2114 rt[i]=r;
2115 continue;
2116 }
2117 if (dR>m_drmax) break;
2118 if ((*r)->sur()==sur0) continue;
2119
2120 // Comparison with bottom and top Z
2121 //
2122 float Tz = ((*r)->z()-Z)/dR;
2123 float Zo = Z-R*Tz;
2124 if (Zo < data.zminB || Zo > data.zmaxB) continue;
2125 float Zu = Z+(550.f-R)*Tz;
2126 if (Zu < data.zminU || Zu > data.zmaxU) continue;
2127 data.SP[Nt] = (*r);
2128 if (++Nt==m_maxsizeSP) goto breakt;
2129 }
2130 }
2131
2132 breakt:
2133 if (!(Nt-Nb)) continue;
2134 float covr0 = (*r0)->covr ();
2135 float covz0 = (*r0)->covz ();
2136
2137 float ax = X/R;
2138 float ay = Y/R;
2139
2140 for (int i=0; i<Nt; ++i) {
2141
2142 InDet::SiSpacePointForSeed* sp = data.SP[i];
2143
2144 float dx = sp->x()-X;
2145 float dy = sp->y()-Y;
2146 float dz = sp->z()-Z;
2147 float x = dx*ax+dy*ay;
2148 float y = dy*ax-dx*ay;
2149 float r2 = 1.f/(x*x+y*y);
2150 float dr = std::sqrt(r2);
2151 float tz = dz*dr;
2152 if (i < Nb) tz = -tz;
2153
2154 data.Tz[i] = tz;
2155 data.Zo[i] = Z-R*tz;
2156 data.R [i] = dr;
2157 data.U [i] = x*r2;
2158 data.V [i] = y*r2;
2159 data.Er[i] = ((covz0+sp->covz())+(tz*tz)*(covr0+sp->covr()))*r2;
2160 }
2161 covr0 *= .5f;
2162 covz0 *= 2.f;
2163
2164 // Three space points comparison
2165 //
2166 for (int b=0; b<Nb; ++b) {
2167
2168 float Zob = data.Zo[b];
2169 float Tzb = data.Tz[b];
2170 float Rb2r = data.R [b]*covr0;
2171 float Rb2z = data.R [b]*covz0;
2172 float Erb = data.Er[b];
2173 float Vb = data.V [b];
2174 float Ub = data.U [b];
2175 float Tzb2 = (1.f+Tzb*Tzb);
2176 float CSA = Tzb2*COFK;
2177 float ICSA = Tzb2*ipt2C;
2178 float d0max = maxd0cut;
2179 if (data.SP[b]->spacepoint->clusterList().second) d0max = maxd0cutstrips;
2180
2181 for (int t=Nb; t<Nt; ++t) {
2182
2183 float dT = ((Tzb-data.Tz[t])*(Tzb-data.Tz[t])-data.R[t]*Rb2z-(Erb+data.Er[t]))-(data.R[t]*Rb2r)*((Tzb+data.Tz[t])*(Tzb+data.Tz[t]));
2184 if ( dT > ICSA) continue;
2185
2186 float deltaU = data.U[t]-Ub;
2187 if (deltaU == 0.) continue;
2188 float A = (data.V[t]-Vb)/deltaU;
2189 float onePlusAsquare = 1.f+A*A;
2190 float B = Vb-A*Ub;
2191 float BSquare = B*B;
2192 if (BSquare > ipt2K*onePlusAsquare || dT*onePlusAsquare > BSquare*CSA) continue;
2193
2194 float Im = std::abs((A-B*R)*R);
2195 if (Im > d0max) continue;
2196
2197 // Azimuthal angle test
2198 //
2199 float y = 1.;
2200 float x = 2.f*B*R-A;
2201 float df = std::abs(std::atan2(ay*y-ax*x,ax*y+ay*x)-data.ftrig);
2202 if (df > M_PI) df = twoPi-df;
2203 if (df > data.ftrigW) continue;
2204 data.CmSp.emplace_back(B/std::sqrt(onePlusAsquare), data.SP[t]);
2205 data.SP[t]->setParam(Im);
2206 }
2207 if (!data.CmSp.empty()) {
2208 newOneSeedWithCurvaturesComparison(data, data.SP[b], (*r0), Zob);
2209 }
2210 }
2211 fillSeeds(data);
2212 nseed += data.fillOneSeeds;
2213 }
2214}
2215
2217// New 3 space points pro seeds
2219
2221(EventData& data,
2223 InDet::SiSpacePointForSeed*& p3, float z, float seedCandidateQuality) const
2224{
2226 float worstQualityInMap = std::numeric_limits<float>::min();
2227 InDet::SiSpacePointsProSeed* worstSeedSoFar = nullptr;
2228 if (!data.mapOneSeeds_Pro.empty()) {
2229 std::multimap<float,InDet::SiSpacePointsProSeed*>::reverse_iterator l = data.mapOneSeeds_Pro.rbegin();
2230 worstQualityInMap = (*l).first;
2231 worstSeedSoFar = (*l).second;
2232 }
2235 if (data.nOneSeeds < data.maxSeedsPerSP
2238 || (data.keepAllConfirmedSeeds && worstQualityInMap <= seedCandidateQuality && isConfirmedSeed(p1,p3,seedCandidateQuality) && data.nOneSeeds < data.seedPerSpCapacity)
2241 || (data.keepAllConfirmedSeeds && worstQualityInMap > seedCandidateQuality && worstSeedSoFar
2242 && isConfirmedSeed(worstSeedSoFar->spacepoint0(),worstSeedSoFar->spacepoint2(),worstQualityInMap) && data.nOneSeeds < data.seedPerSpCapacity)
2243 ){
2244 data.OneSeeds_Pro[data.nOneSeeds].set(p1,p2,p3,z);
2245 data.mapOneSeeds_Pro.insert(std::make_pair(seedCandidateQuality, &data.OneSeeds_Pro[data.nOneSeeds]));
2246 ++data.nOneSeeds;
2247 }
2249 else if (worstSeedSoFar && (worstQualityInMap > seedCandidateQuality)){
2251 worstSeedSoFar->set(p1,p2,p3,z);
2253 std::multimap<float,InDet::SiSpacePointsProSeed*>::iterator
2254 i = data.mapOneSeeds_Pro.insert(std::make_pair(seedCandidateQuality,worstSeedSoFar));
2256 for (++i; i!=data.mapOneSeeds_Pro.end(); ++i) {
2257 if ((*i).second==worstSeedSoFar) {
2258 data.mapOneSeeds_Pro.erase(i);
2259 return;
2260 }
2261 }
2262 }
2263}
2264
2266// New 3 space points pro seeds production
2268namespace {
2269 inline
2270 float computeEta( float r, float z) {
2271 if (r <= 10e-9) return 0;
2272 auto asinh = [] (double x) { return std::log (std::sqrt (x*x+1) + x); };
2273 return asinh (z / r);
2274 }
2275}
2276
2278(EventData& data, SiSpacePointForSeed*& SPb, SiSpacePointForSeed*& SP0, float Zob) const
2279{
2280 constexpr float curvatureInterval = .00003;
2281
2282 bool bottomSPisPixel = !SPb->spacepoint->clusterList().second;
2283 float bottomSPQuality = SPb->quality();
2284 float centralSPQuality = SP0->quality();
2285
2287 if(data.CmSp.size() > 2) std::sort(data.CmSp.begin(), data.CmSp.end(), comCurvature());
2288
2289 float bottomR=SPb->radius();
2290 float bottomZ=SPb->z();
2291
2292 std::vector<std::pair<float,InDet::SiSpacePointForSeed*>>::iterator it_otherSP;
2293 std::vector<std::pair<float,InDet::SiSpacePointForSeed*>>::iterator it_commonTopSP = data.CmSp.begin(), ie = data.CmSp.end();
2294 std::vector<std::pair<float,InDet::SiSpacePointForSeed*>>::iterator it_startInnerLoop=it_commonTopSP;
2295
2297 for (; it_commonTopSP!=ie; ++it_commonTopSP) {
2298
2300 float seedIP = (*it_commonTopSP).second->param();
2301 float seedQuality = seedIP + (*it_commonTopSP).second->scorePenalty();
2302 float originalSeedQuality = seedQuality;
2303
2304 if(m_maxdImpact > 50){ //This only applies to LRT
2305 // Tell clang to optimize assuming that FP exceptions can trap.
2306 // Otherwise, it can vectorize the division, which can lead to
2307 // spurious division-by-zero traps from unused vector lanes.
2309
2310 float topR=(*it_commonTopSP).second->radius();
2311 float topZ=(*it_commonTopSP).second->z();
2312
2313 float Zot = std::abs(topR - bottomR) > 10e-9 ?
2314 bottomZ - (bottomR - originalSeedQuality) * ((topZ - bottomZ) / (topR - bottomR)) : bottomZ;
2315
2316 float eta1 = computeEta(topR - bottomR, topZ - bottomZ);
2317 float eta0 = computeEta(seedIP, Zot);
2318
2319 float deltaEta=std::abs(eta1-eta0); //For LLP daughters, the direction of the track is correlated with the direction of the LLP (which is correlated with the direction of the point of closest approach
2320 //calculate weighted average of d0 and deltaEta, normalized by their maximum values
2321 float f=std::min(0.5f,originalSeedQuality/200.f); //0.5 and 200 are parameters chosen from a grid scan to optimize efficiency
2322 seedQuality*=(1.f-f)/300.f;
2323 seedQuality+=f*deltaEta/2.5f;
2324 }
2325
2326 bool topSPisPixel = !(*it_commonTopSP).second->spacepoint->clusterList().second;
2327
2329 const Trk::Surface* surfaceTopSP = (*it_commonTopSP).second->sur ();
2330 float radiusTopSP = (*it_commonTopSP).second->radius();
2332 float minCurvature =(*it_commonTopSP).first-curvatureInterval;
2333 float maxCurvature =(*it_commonTopSP).first+curvatureInterval;
2334
2341
2343 if (!bottomSPisPixel) seedQuality+=m_seedScoreBonusSSS;
2345 else if ( topSPisPixel) seedQuality+=m_seedScoreBonusPPP;
2346
2347
2353
2354 for (it_otherSP=it_startInnerLoop; it_otherSP!=ie; ++it_otherSP) {
2356 if ( it_otherSP == it_commonTopSP ) continue;
2359 if ( (*it_otherSP).first < minCurvature ) {
2360 it_startInnerLoop=it_otherSP;
2361 ++it_startInnerLoop;
2362 continue;
2363 }
2365 if ( (*it_otherSP).first > maxCurvature ) break;
2367 if ( (*it_otherSP).second->sur()==surfaceTopSP) continue;
2369 float radiusOtherSP = (*it_otherSP).second->radius();
2370 if (std::abs(radiusOtherSP-radiusTopSP) < m_drmin) continue;
2371 // if we have a confirmation seed, we improve the score of the seed.
2372 seedQuality += m_seedScoreBonusConfirmationSeed;
2373 // apply confirmation bonus only once
2374 break;
2375 }
2376
2378 if (seedQuality > data.maxScore) continue;
2379
2382 if (bottomSPisPixel!=topSPisPixel) {
2383 if (seedQuality > 0. ||
2384 (seedQuality > bottomSPQuality && seedQuality > centralSPQuality && seedQuality > (*it_commonTopSP).second->quality())
2385 ) continue;
2386 }
2389 if (!isConfirmedSeed(SPb,it_commonTopSP->second,seedQuality)){
2391 double maxdImpact = m_maxdImpact - (m_dImpactCutSlopeUnconfirmedPPP * (*it_commonTopSP).second->scorePenalty());
2393 if (!bottomSPisPixel) maxdImpact = m_maxdImpactSSS - (m_dImpactCutSlopeUnconfirmedSSS * (*it_commonTopSP).second->scorePenalty());
2394 if (seedIP > maxdImpact) continue;
2395 }
2397 newOneSeed(data, SPb, SP0, (*it_commonTopSP).second, Zob, seedQuality);
2398 }
2399 data.CmSp.clear();
2400}
2401
2405
2407{
2408
2409 data.fillOneSeeds = 0;
2410
2411 std::multimap<float,InDet::SiSpacePointsProSeed*>::iterator it_firstSeedCandidate = data.mapOneSeeds_Pro.begin();
2412 std::multimap<float,InDet::SiSpacePointsProSeed*>::iterator it_seedCandidate = data.mapOneSeeds_Pro.begin();
2413 std::multimap<float,InDet::SiSpacePointsProSeed*>::iterator it_endSeedCandidates = data.mapOneSeeds_Pro.end();
2414
2416 if (it_seedCandidate==it_endSeedCandidates) return;
2417
2418 SiSpacePointsProSeed* theSeed{nullptr};
2419
2421 for (; it_seedCandidate!=it_endSeedCandidates; ++it_seedCandidate) {
2422
2424 float quality = (*it_seedCandidate).first;
2425 theSeed = (*it_seedCandidate).second;
2427 if (it_seedCandidate!=it_firstSeedCandidate && theSeed->spacepoint0()->radius() < m_radiusCutIBL && quality > m_seedScoreThresholdPPPConfirmationSeed) continue;
2429 if (!theSeed->setQuality(quality)) continue;
2430
2432 if (data.i_seede_Pro!=data.l_seeds_Pro.end()) {
2433 theSeed = &(*data.i_seede_Pro++);
2434 *theSeed = *(*it_seedCandidate).second;
2435 } else {
2437 data.l_seeds_Pro.emplace_back(*(*it_seedCandidate).second);
2438 //unused value, keep commented to avoid repetition
2439 //theSeed = &(data.l_seeds_Pro.back());
2440 data.i_seede_Pro = data.l_seeds_Pro.end();
2441 }
2442
2443 ++data.fillOneSeeds;
2444 }
2445}
2446
2448{
2450 if (not data.initialized) initializeEventData(data);
2451
2452 if (data.nspoint==3) {
2453 do {
2455 if (data.i_seed_Pro==data.i_seede_Pro) {
2460 findNext(data);
2462 //cppcheck-suppress identicalInnerCondition
2463 if (data.i_seed_Pro==data.i_seede_Pro) return nullptr;
2464 }
2466 } while (!(*data.i_seed_Pro++).set3(data.seedOutput));
2467
2469 return &data.seedOutput;
2470 } else {
2472 if (data.i_seed_Pro==data.i_seede_Pro) {
2473 findNext(data);
2474 //cppcheck-suppress identicalInnerCondition
2475 if (data.i_seed_Pro==data.i_seede_Pro) return nullptr;
2476 }
2477 (*data.i_seed_Pro++).set2(data.seedOutput);
2478 return &data.seedOutput;
2479 }
2480 return nullptr;
2481}
2482
2486void InDet::SiSpacePointsSeedMaker_ATLxk::writeNtuple(const SiSpacePointsSeed* seed, const Trk::Track* track, int seedType, long eventNumber) const{
2487
2488 if(m_writeNtuple) {
2489 std::lock_guard<std::mutex> lock(m_mutex);
2490
2491 if(track != nullptr) {
2492 m_trackPt = (track->trackParameters()->front()->pT())/1000.f;
2493 m_trackEta = std::abs(track->trackParameters()->front()->eta());
2494 }
2495 else {
2496 m_trackPt = -1.;
2497 m_trackEta = -1.;
2498 }
2499 m_d0 = seed->d0();
2500 m_z0 = seed->zVertex();
2501 m_eta = seed->eta();
2502 m_x1 = seed->x1();
2503 m_x2 = seed->x2();
2504 m_x3 = seed->x3();
2505 m_y1 = seed->y1();
2506 m_y2 = seed->y2();
2507 m_y3 = seed->y3();
2508 m_z1 = seed->z1();
2509 m_z2 = seed->z2();
2510 m_z3 = seed->z3();
2511 m_r1 = seed->r1();
2512 m_r2 = seed->r2();
2513 m_r3 = seed->r3();
2514 m_type = seedType;
2515 m_dzdr_b = seed->dzdr_b();
2516 m_dzdr_t = seed->dzdr_t();
2517 m_pt = seed->pt();
2518 m_givesTrack = !(track == nullptr);
2519 m_eventNumber = eventNumber;
2520
2521 // Ok: protected by mutex.
2522 TTree* outputTree ATLAS_THREAD_SAFE = m_outputTree;
2523 outputTree->Fill();
2524
2525 }
2526
2527}
2528
2530(EventData& data, const float Zv, const float R, const float T) const
2531{
2532 if (Zv < data.zminU or Zv > data.zmaxU) return false;
2533 if (not data.isvertex) return true;
2534 if (data.l_vertex.empty()) return false;
2535
2536 float dZmin = std::numeric_limits<float>::max();
2537 for (const float& v: data.l_vertex) {
2538 float dZ = std::abs(v-Zv);
2539 if (dZ >= dZmin) break;
2540 dZmin = dZ;
2541 }
2542
2543 //return dZmin < (m_dzver+m_dzdrver*R)*sqrt(1.+T*T);
2544 //(Minor) speed-up: Avoid calculation of sqrt, compare squares
2545 return dZmin*dZmin < (m_dzver+m_dzdrver*R)*(m_dzver+m_dzdrver*R)*(1.f+T*T);
2546}
2547
2551
2553(EventData& data, const Trk::SpacePoint*const& sp) const
2554{
2555 InDet::SiSpacePointForSeed* sps = nullptr;
2556
2559 std::array<float,3> r{0,0,0};
2560 convertToBeamFrameWork(data, sp, r);
2561
2563 if (data.checketa) {
2564 float z = (std::abs(r[2])+m_zmax);
2565 float x = r[0]*data.dzdrmin;
2566 float y = r[1]*data.dzdrmin;
2567 if ((z*z )<(x*x+y*y)) return sps;
2568 }
2572 if (data.i_spforseed!=data.l_spforseed.end()) {
2574 sps = &(*data.i_spforseed++);
2577 sps->set(sp,r);
2578 } else {
2580 data.l_spforseed.emplace_back(sp, r);
2582 sps = &(data.l_spforseed.back());
2584 data.i_spforseed = data.l_spforseed.end();
2585 }
2586
2587 return sps;
2588}
2589
2591// New 2 space points seeds
2593
2596{
2597 InDet::SiSpacePointForSeed* p3 = nullptr;
2598
2599 if (data.i_seede_Pro!=data.l_seeds_Pro.end()) {
2600 SiSpacePointsProSeed* s = &(*data.i_seede_Pro++);
2601 s->set(p1, p2, p3, z);
2602 } else {
2603 data.l_seeds_Pro.emplace_back(p1, p2, p3, z);
2604 data.i_seede_Pro = data.l_seeds_Pro.end();
2605 }
2606}
2607
2609 int seedArrayPerSPSize = (m_maxOneSizePPP>m_maxOneSizeSSS ? m_maxOneSizePPP : m_maxOneSizeSSS);
2611 data.initialize(EventData::ToolType::ATLxk,
2613 seedArrayPerSPSize,
2614 0,
2615 m_nBinsR,
2616 0,
2619 m_checketa);
2620}
2621
2623
2625 if (bottomSP->spacepoint->clusterList().second){
2626 return (quality < m_seedScoreThresholdSSSConfirmationSeed);
2627 }
2629 else if (!topSP->spacepoint->clusterList().second){
2630 return (quality < m_seedScoreThresholdPPPConfirmationSeed);
2631 }
2633 else return (quality < 0.);
2634}
2635
const std::regex re(r_e)
#define M_PI
Scalar phi() const
phi method
Scalar theta() const
theta method
#define endmsg
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_ERROR(x,...)
virtual void lock()=0
Interface to allow an object to lock itself when made const in SG.
static Double_t sp
static Double_t a
static Double_t sc
#define X(x)
Definition SiHitPlots.h:102
#define y
#define x
#define z
Define macros for attributes used to control the static checker.
void getInitializedCache(MagField::AtlasFieldCache &cache) const
get B field cache for evaluation as a function of 2-d or 3-d position.
Describes the API of the Region of Ineterest geometry.
virtual double phiPlus() const =0
extreme phi values
virtual double zedPlus() const =0
the zed and eta values at the most forward and most rear ends of the RoI
virtual double phiMinus() const =0
virtual double zedMinus() const =0
virtual double etaMinus() const =0
virtual double etaPlus() const =0
This is a "hash" representation of an Identifier.
Class to hold geometrical description of a silicon detector element.
const Trk::SpacePoint * spacepoint
float quality() const
penalty term in the seed score
void set(const Trk::SpacePoint *, std::span< float const, 3 >)
SiSpacePointForSeed * spacepoint0()
void set(SiSpacePointForSeed *&, SiSpacePointForSeed *&, SiSpacePointForSeed *&, float)
SiSpacePointForSeed * spacepoint2()
FloatProperty m_maxScore
Maximum score to accept.
virtual bool getWriteNtupleBoolProperty() const override
static void newSeed(EventData &data, SiSpacePointForSeed *&p1, SiSpacePointForSeed *&p2, float z)
SG::ReadHandleKey< SpacePointContainer > m_spacepointsSCT
virtual void find3Sp(const EventContext &ctx, EventData &data, const std::list< Trk::Vertex > &lv) const override
with three space points with or without vertex constraint
void newOneSeed(EventData &data, SiSpacePointForSeed *&p1, SiSpacePointForSeed *&p2, SiSpacePointForSeed *&p3, float z, float quality) const
This inserts a seed into the set of saved seeds.
std::array< std::array< int, arraySizeNeighbourBinsVertex >, arraySizePhiZV > m_neighboursVertexPhiZ
virtual void newRegion(const EventContext &ctx, EventData &data, const std::vector< IdentifierHash > &vPixel, const std::vector< IdentifierHash > &vSCT) const override
Initialize tool for new region.
void fillSeeds(EventData &data) const
fills the seeds from the mapOneSeeds_Pro member into the l_seeds_Pro member of the data object,...
virtual void find2Sp(EventData &data, const std::list< Trk::Vertex > &lv) const override
With two space points with or without vertex constraint.
float m_seedScoreThresholdPPPConfirmationSeed
Seed score thresholds defined based on the modifiers defined as configurables above.
SG::ReadHandleKey< SpacePointContainer > m_spacepointsPixel
int m_maxBinPhiVertex
number of bins in phi for vertices
void buildFrameWork()
prepare several data members with cached cut values, conversion factors, binnings,...
BooleanProperty m_optimisePhiBinning
This flag will make the buildFrameWork method determine an optimal phi binning of the search regions ...
bool newVertices(EventData &data, const std::list< Trk::Vertex > &lV) const
This method updates the EventData based on the passed list of vertices.
SG::ReadCondHandleKey< InDet::BeamSpotData > m_beamSpotKey
std::array< int, arraySizePhiZ > m_nNeighbourCellsBottom
arrays associating bins to each other for SP formation
FloatProperty m_dImpactCutSlopeUnconfirmedSSS
these flags allow to dynamically tighten the d0 cut on non-confirmed seeds based on the penalty score...
void fillLists(EventData &data) const
this method populates the data object's "histograms" (implemented as nested vectors).
MsgStream & dumpConditions(EventData &data, MsgStream &out) const
Dumps conditions information into the MsgStream.
void production3SpTrigger(EventData &data, std::array< std::vector< InDet::SiSpacePointForSeed * >::iterator, arraySizeNeighbourBins > &rb, std::array< std::vector< InDet::SiSpacePointForSeed * >::iterator, arraySizeNeighbourBins > &rbe, std::array< std::vector< InDet::SiSpacePointForSeed * >::iterator, arraySizeNeighbourBins > &rt, std::array< std::vector< InDet::SiSpacePointForSeed * >::iterator, arraySizeNeighbourBins > &rte, const int numberBottomCells, const int numberTopCells, int &nseed) const
as above, but for the trigger
void findNext(EventData &data) const
This method is called within next() when we are out of vertices.
static void convertToBeamFrameWork(EventData &data, const Trk::SpacePoint *const &sp, std::array< float, 3 > &r)
This method popualtes the r array with the space point's coordinates relative to the beam spot.
SiSpacePointForSeed * newSpacePoint(EventData &data, const Trk::SpacePoint *const &sp) const
Create a SiSpacePointForSeed from the space point.
static constexpr float m_radiusCutIBL
We detect IBL hits via the seed radial location.
virtual void findVSp(const EventContext &ctx, EventData &data, const std::list< Trk::Vertex > &lv) const override
with variable number space points with or without vertex constraint Variable means (2,...
float m_seedScoreThresholdSSSConfirmationSeed
max (score is assigned negative sign) score for SSS seeds with confirmation seed requirement.
FloatProperty m_seedScoreBonusPPP
Scoring modifiers applied when ranking seeds.
bool isConfirmedSeed(const InDet::SiSpacePointForSeed *bottomSP, const InDet::SiSpacePointForSeed *topSP, float quality) const
Helper method to determine if a seed is 'confirmed' - this means that a second seed exists with compa...
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_fieldCondObjInputKey
Read handle for conditions object to get the field cache.
float m_inverseBinSizePhiVertex
as above but for vertex
std::array< int, arraySizePhiZV > m_nNeighboursVertexPhiZ
SG::ReadHandleKey< SpacePointOverlapCollection > m_spacepointsOverlap
virtual void newEvent(const EventContext &ctx, EventData &data, int iteration) const override
Initialize tool for new event.
float m_inverseBinSizePhi
cache the inverse bin size in phi which we use - needed to evaluate phi bin locations
BooleanProperty m_alwaysKeepConfirmedPixelSeeds
This flag will lead to all confirmed seeds (seeds where a second compatible seed with a different top...
std::array< std::array< int, arraySizeNeighbourBins >, arraySizePhiZ > m_neighbourCellsBottom
mapping of neighbour cells in the 2D phi-z binning to consider for the "bottom SP" search for central...
bool isUsed(const Trk::SpacePoint *sp, const Trk::PRDtoTrackMap &prd_to_track_map) const
int m_nBinsR
number of bins in the radial coordinate
virtual StatusCode initialize() override
Initialisation.
IntegerProperty m_maxOneSizeSSS
maximum number of seeds to keep per central space point.
void buildBeamFrameWork(EventData &data) const
Initiate beam frame work for seed generator.
float m_ipt2
inverse square of 90% of the pt min cut
SG::ReadHandleKey< Trk::PRDtoTrackMap > m_prdToTrackMap
virtual StatusCode finalize() override
Finalize.
static constexpr float m_COF
appears to be an approximated term related to multiple-scattering of particles traversing the ID duri...
void newOneSeedWithCurvaturesComparison(EventData &data, SiSpacePointForSeed *&SPb, SiSpacePointForSeed *&SP0, float Zob) const
This creates all possible seeds with the passed central and bottom SP, using all top SP candidates wh...
virtual const SiSpacePointsSeed * next(const EventContext &ctx, EventData &data) const override
This method will update the data.seedOutput member to be the next seed pointed at by the data....
virtual MsgStream & dump(EventData &data, MsgStream &out) const override
Dumps relevant information into the MsgStream.
std::array< int, arraySizePhiZ > m_nNeighbourCellsTop
number of neighbouring phi-z bins to consider when looking for "top SP" candidates for each phi-z bin
bool isZCompatible(EventData &data, const float Zv, const float R, const float T) const
std::array< std::array< int, arraySizeNeighbourBins >, arraySizePhiZ > m_neighbourCellsTop
mapping of neighbour cells in the 2D phi-z binning to consider for the "top SP" search for central SP...
Gaudi::Property< bool > m_writeNtuple
Flag to write validation ntuples. Turned off by default.
float m_ipt
inverse of 90% of the ptmin cut
@ arraySizePhiZ
capacity for the 2D phi-z arrays
@ arraySizeNeighbourBins
array size to store neighbouring phi-z-regions in the seed finding
@ arraySizePhiZV
array size in phi-Z 2D for the vertexing
@ arraySizePhiV
array size in phi for vertexing
void production3Sp(EventData &data) const
Top-level method for 3-SP seed production.
static MsgStream & dumpEvent(EventData &data, MsgStream &out)
Dumps event information into the MsgStream.
float m_dzdrmin0
conversion factors and cached cut values
static float azimuthalStep(const float pTmin, const float maxd0, const float Rmin, const float Rmax)
Determine the expected azimuthal trajectory displacement in phi in presence of the magnetic field for...
virtual void writeNtuple(const SiSpacePointsSeed *seed, const Trk::Track *track, int seedType, long eventNumber) const override
This method is called by the SiSPSeededTrackFinder algorithm to fill ntuples for seeds seen by the al...
SiSpacePointsSeedMakerEventData EventData
Local cache for magnetic field (based on MagFieldServices/AtlasFieldSvcTLS.h).
bool solenoidOn() const
status of the magnets
void getFieldZR(const double *ATH_RESTRICT xyz, double *ATH_RESTRICT bxyz, double *ATH_RESTRICT deriv=nullptr)
get B field valaue on the z-r plane at given position works only inside the solenoid.
virtual bool isValid() override final
Can the handle be successfully dereferenced?
const_pointer_type cptr()
Dereference the pointer.
const std::pair< const PrepRawData *, const PrepRawData * > & clusterList() const
return the pair of cluster pointers by reference
Abstract Base Class for tracking surfaces.
Definition Surface.h:79
This class is a simplest representation of a vertex candidate.
int ir
counter of the current depth
Definition fastadd.cxx:49
int r
Definition globals.cxx:22
Eigen::Matrix< double, 3, 1 > Vector3D
-event-from-file
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
hold the test vectors and ease the comparison
MsgStream & msg
Definition testRead.cxx:32
Tell the compiler to optimize assuming that FP may trap.
#define CXXUTILS_TRAPPING_FP
Definition trapping_fp.h:24