ATLAS Offline Software
Loading...
Searching...
No Matches
SiSpacePointsSeedMaker_Trigger.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3*/
4
6// Implementation file for class SiSpacePointsSeedMaker_Trigger
8// (c) ATLAS Detector software
10// AlgTool used for TRT_DriftCircleOnTrack object production
12// Version 1.0 21/04/2004 I.Gavrilenko
14
16
17#include <cmath>
18
19#include <iomanip>
20#include <ostream>
21
23// Constructor
25
27(const std::string& t, const std::string& n, const IInterface* p)
28 : base_class(t, n, p)
29{
30}
31
33// Initialisation
35
37{
38 StatusCode sc = AlgTool::initialize();
39
41 ATH_CHECK(m_spacepointsSCT.initialize(m_sct));
43
44 // Get beam geometry
45 //
46 ATH_CHECK(m_beamSpotKey.initialize());
47
49
50 // Build framework
51 //
53
54 // Get output print level
55 //
56 m_outputlevel = msg().level()-MSG::DEBUG;
57 if (m_outputlevel<=0) {
60 data.nprint=0;
61 dump(data, msg(MSG::DEBUG));
62 }
63
64 m_initialized = true;
65
66 return sc;
67}
68
70// Finalize
72
74{
75 return AlgTool::finalize();
76}
77
79// Initialize tool for new event
81
82void InDet::SiSpacePointsSeedMaker_Trigger::newEvent(const EventContext& ctx, EventData& data, int) const
83{
84 if (not data.initialized) initializeEventData(data);
85
86 data.trigger = false;
87 if (!m_pixel && !m_sct) return;
88 erase(data);
90
91 double f[3], gP[3] ={10.,10.,0.};
92
94
95 // Get field cache object
97 const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
98 if (fieldCondObj == nullptr) {
99 ATH_MSG_ERROR("SiSpacePointsSeedMaker_Trigger: Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
100 return;
101 }
102 fieldCondObj->getInitializedCache (fieldCache);
103
104 if (fieldCache.solenoidOn()) {
105 fieldCache.getFieldZR(gP, f);
106
107 data.K = 2./(300.*f[2]);
108 } else {
109 data.K = 2./(300.* 5.);
110 }
111
112 data.i_spforseed = data.l_spforseed.begin();
113
114 float irstep = 1.f/m_r_rstep;
115 int irmax = m_r_size-1;
116
117 // Get pixels space points containers from store gate
118 //
119 if (m_pixel) {
120
122 if (spacepointsPixel.isValid()) {
123
124 for (const SpacePointCollection* spc: *spacepointsPixel) {
125 for (const Trk::SpacePoint* sp: *spc) {
126
127 float r = sp->r();
128 if (r<0. || r>=m_r_rmax) continue;
129
131
132 int ir = static_cast<int>(sps->radius()*irstep);
133 if (ir>irmax) ir = irmax;
134 data.r_Sorted[ir].push_back(sps);
135 ++data.r_map[ir];
136 if (data.r_map[ir]==1) data.r_index[data.nr++] = ir;
137 ++data.ns;
138 }
139 }
140 }
141 }
142
143 // Get sct space points containers from store gate
144 //
145 if (m_sct) {
146
148 if (spacepointsSCT.isValid()) {
149
150 for (const SpacePointCollection* spc: *spacepointsSCT) {
151 for (const Trk::SpacePoint* sp: *spc) {
152
153 float r = sp->r();
154 if (r<0. || r>=m_r_rmax) continue;
155
157
158 int ir = static_cast<int>(sps->radius()*irstep);
159 if (ir>irmax) ir = irmax;
160 data.r_Sorted[ir].push_back(sps);
161 ++data.r_map[ir];
162 if (data.r_map[ir]==1) data.r_index[data.nr++] = ir;
163 ++data.ns;
164 }
165 }
166 }
167
168 // Get sct overlap space points containers from store gate
169 //
170 if (m_useOverlap) {
171
173 if (spacepointsOverlap.isValid()) {
174
175 for (const Trk::SpacePoint* sp: *spacepointsOverlap) {
176
177 float r = sp->r();
178 if (r<0. || r>=m_r_rmax) continue;
179
181
182 int ir = static_cast<int>(sps->radius()*irstep);
183 if (ir>irmax) ir = irmax;
184 data.r_Sorted[ir].push_back(sps);
185 ++data.r_map[ir];
186 if (data.r_map[ir]==1) data.r_index[data.nr++] = ir;
187 ++data.ns;
188 }
189 }
190 }
191 }
193}
194
196// Initialize tool for new region
198
200(const EventContext& ctx, EventData& data, const std::vector<IdentifierHash>& vPixel, const std::vector<IdentifierHash>& vSCT) const
201{
202 if (not data.initialized) initializeEventData(data);
203
204 data.trigger = false;
205 if (!m_pixel && !m_sct) return;
206 erase(data);
207
209
210 double f[3], gP[3] ={10.,10.,0.};
211
212 MagField::AtlasFieldCache fieldCache;
213
214 // Get field cache object
216 const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
217 if (fieldCondObj == nullptr) {
218 ATH_MSG_ERROR("SiSpacePointsSeedMaker_Trigger: Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
219 return;
220 }
221 fieldCondObj->getInitializedCache (fieldCache);
222
223 if (fieldCache.solenoidOn()) {
224 fieldCache.getFieldZR(gP, f);
225
226 data.K = 2./(300.*f[2]);
227 } else {
228 data.K = 2./(300.* 5.);
229 }
230
231 data.i_spforseed = data.l_spforseed.begin();
232
233 float irstep = 1.f/m_r_rstep;
234 int irmax = m_r_size-1;
235
236 // Get pixels space points containers from store gate
237 //
238 if (m_pixel && !vPixel.empty()) {
239
241 if (spacepointsPixel.isValid()) {
242
243 // Loop through all trigger collections
244 //
245 for (const IdentifierHash& l: vPixel) {
246 const auto *w = spacepointsPixel->indexFindPtr(l);
247 if (w==nullptr) continue;
248 for (const Trk::SpacePoint* sp: *w) {
249 float r = sp->r();
250 if (r<0. || r>=m_r_rmax) continue;
252 int ir = static_cast<int>(sps->radius()*irstep);
253 if (ir>irmax) ir = irmax;
254 data.r_Sorted[ir].push_back(sps);
255 ++data.r_map[ir];
256 if (data.r_map[ir]==1) data.r_index[data.nr++] = ir;
257 ++data.ns;
258 }
259 }
260 }
261 }
262
263 // Get sct space points containers from store gate
264 //
265 if (m_sct && !vSCT.empty()) {
266
268 if (spacepointsSCT.isValid()) {
269
270 // Loop through all trigger collections
271 //
272 for (const IdentifierHash& l: vSCT) {
273 const auto *w = spacepointsSCT->indexFindPtr(l);
274 if (w==nullptr) continue;
275 for (const Trk::SpacePoint* sp: *w) {
276 float r = sp->r();
277 if (r<0. || r>=m_r_rmax) continue;
279 int ir = static_cast<int>(sps->radius()*irstep);
280 if (ir>irmax) ir = irmax;
281 data.r_Sorted[ir].push_back(sps);
282 ++data.r_map[ir];
283 if (data.r_map[ir]==1) data.r_index[data.nr++] = ir;
284 ++data.ns;
285 }
286 }
287 }
288 }
290}
291
293// Initialize tool for new region
295
297(const EventContext& ctx, EventData& data,
298 const std::vector<IdentifierHash>& vPixel, const std::vector<IdentifierHash>& vSCT,
299 const IRoiDescriptor& IRD) const
300{
301 if (not data.initialized) initializeEventData(data);
302
303 newRegion(ctx, data, vPixel, vSCT);
304
305 data.trigger = true;
306
307 double dzdrmin = 1./tan(2.*atan(exp(-IRD.etaMinus())));
308 double dzdrmax = 1./tan(2.*atan(exp(-IRD.etaPlus ())));
309
310 data.zminB = IRD.zedMinus()-data.zbeam[0]; // min bottom Z
311 data.zmaxB = IRD.zedPlus ()-data.zbeam[0]; // max bottom Z
312 data.zminU = data.zminB+550.f*dzdrmin;
313 data.zmaxU = data.zmaxB+550.f*dzdrmax;
314 double fmax = IRD.phiPlus ();
315 double fmin = IRD.phiMinus();
316 if (fmin > fmax) fmin-=(2.*M_PI);
317 data.ftrig = (fmin+fmax)*.5;
318 data.ftrigW = (fmax-fmin)*.5;
319}
320
322// Methods to initilize different strategies of seeds production
323// with two space points with or without vertex constraint
325
326void InDet::SiSpacePointsSeedMaker_Trigger::find2Sp(EventData& data, const std::list<Trk::Vertex>& lv) const
327{
328 if (not data.initialized) initializeEventData(data);
329
330 int mode = 0;
331 if (lv.begin()!=lv.end()) mode = 1;
332 bool newv = newVertices(data, lv);
333
334 if (newv || !data.state || data.nspoint!=2 || data.mode!=mode || data.nlist) {
335
336 data.i_seede = data.l_seeds.begin();
337 data.state = 1;
338 data.nspoint = 2;
339 data.nlist = 0;
340 data.mode = mode;
341 data.endlist = true;
342 data.fvNmin = 0;
343 data.fNmin = 0;
344 data.zMin = 0;
346 }
347 data.i_seed = data.l_seeds.begin();
348
349 if (m_outputlevel<=0) {
350 data.nprint=1;
351 dump(data, msg(MSG::DEBUG));
352 }
353}
354
356// Methods to initilize different strategies of seeds production
357// with three space points with or without vertex constraint
359
360void InDet::SiSpacePointsSeedMaker_Trigger::find3Sp(const EventContext&, EventData& data, const std::list<Trk::Vertex>& lv) const
361{
362 if (not data.initialized) initializeEventData(data);
363
364 int mode = 2;
365 if (lv.begin()!=lv.end()) mode = 3;
366 bool newv = newVertices(data, lv);
367
368 if (newv || !data.state || data.nspoint!=3 || data.mode!=mode || data.nlist) {
369
370 data.i_seede = data.l_seeds.begin();
371 data.state = 1;
372 data.nspoint = 3;
373 data.nlist = 0;
374 data.mode = mode;
375 data.endlist = true;
376 data.fvNmin = 0;
377 data.fNmin = 0;
378 data.zMin = 0;
380 }
381 data.i_seed = data.l_seeds.begin();
382 data.seed = data.mapSeeds.begin();
383 data.seede = data.mapSeeds.end ();
384
385 if (m_outputlevel<=0) {
386 data.nprint=1;
387 dump(data, msg(MSG::DEBUG));
388 }
389}
390
391void InDet::SiSpacePointsSeedMaker_Trigger::find3Sp(const EventContext& ctx, EventData& data, const std::list<Trk::Vertex>& lv, const double*) const
392{
393 find3Sp(ctx, data, lv);
394}
395
397// Methods to initilize different strategies of seeds production
398// with variable number space points with or without vertex constraint
399// Variable means (2,3,4,....) any number space points
401
402void InDet::SiSpacePointsSeedMaker_Trigger::findVSp(const EventContext&, EventData& data, const std::list<Trk::Vertex>& lv) const
403{
404 if (not data.initialized) initializeEventData(data);
405
406 int mode = 5;
407 if (lv.begin()!=lv.end()) mode = 6;
408 bool newv = newVertices(data, lv);
409
410 if (newv || !data.state || data.nspoint!=4 || data.mode!=mode || data.nlist) {
411
412 data.i_seede = data.l_seeds.begin();
413 data.state = 1;
414 data.nspoint = 4;
415 data.nlist = 0;
416 data.mode = mode;
417 data.endlist = true;
418 data.fvNmin = 0;
419 data.fNmin = 0;
420 data.zMin = 0;
422 }
423 data.i_seed = data.l_seeds.begin();
424 data.seed = data.mapSeeds.begin();
425 data.seede = data.mapSeeds.end ();
426
427 if (m_outputlevel<=0) {
428 data.nprint=1;
429 dump(data, msg(MSG::DEBUG));
430 }
431}
432
434// Dumps relevant information into the MsgStream
436
438{
439 if (not data.initialized) initializeEventData(data);
440
441 if (data.nprint) return dumpEvent(data, out);
442 return dumpConditions(data, out);
443}
444
446// Dumps conditions information into the MsgStream
448
450{
451 int n = 42-m_spacepointsPixel.key().size();
452 std::string s2; for (int i=0; i<n; ++i) s2.append(" "); s2.append("|");
453 n = 42-m_spacepointsSCT.key().size();
454 std::string s3; for (int i=0; i<n; ++i) s3.append(" "); s3.append("|");
455 n = 42-m_spacepointsOverlap.key().size();
456 std::string s4; for (int i=0; i<n; ++i) s4.append(" "); s4.append("|");
457 n = 42-m_beamSpotKey.key().size();
458 std::string s5; for (int i=0; i<n; ++i) s5.append(" "); s5.append("|");
459
460
461 out<<"|---------------------------------------------------------------------|"
462 <<endmsg;
463 out<<"| Pixel space points | "<<m_spacepointsPixel.key() <<s2
464 <<endmsg;
465 out<<"| SCT space points | "<<m_spacepointsSCT.key() <<s3
466 <<endmsg;
467 out<<"| Overlap space points | "<<m_spacepointsOverlap.key() <<s4
468 <<endmsg;
469 out<<"| BeamConditionsService | "<<m_beamSpotKey.key()<<s5
470 <<endmsg;
471 out<<"| usePixel | "
472 <<std::setw(12)<<m_pixel
473 <<" |"<<endmsg;
474 out<<"| useSCT | "
475 <<std::setw(12)<<m_sct
476 <<" |"<<endmsg;
477 out<<"| maxSize | "
478 <<std::setw(12)<<m_maxsize
479 <<" |"<<endmsg;
480 out<<"| maxSizeSP | "
481 <<std::setw(12)<<m_maxsizeSP
482 <<" |"<<endmsg;
483 out<<"| pTmin (mev) | "
484 <<std::setw(12)<<std::setprecision(5)<<m_ptmin
485 <<" |"<<endmsg;
486 out<<"| |eta| <= | "
487 <<std::setw(12)<<std::setprecision(5)<<m_etamax
488 <<" |"<<endmsg;
489 out<<"| max radius SP | "
490 <<std::setw(12)<<std::setprecision(5)<<m_r_rmax
491 <<" |"<<endmsg;
492 out<<"| radius step | "
493 <<std::setw(12)<<std::setprecision(5)<<m_r_rstep
494 <<" |"<<endmsg;
495 out<<"| min Z-vertex position | "
496 <<std::setw(12)<<std::setprecision(5)<<m_zmin
497 <<" |"<<endmsg;
498 out<<"| max Z-vertex position | "
499 <<std::setw(12)<<std::setprecision(5)<<m_zmax
500 <<" |"<<endmsg;
501 out<<"| min radius first SP(3) | "
502 <<std::setw(12)<<std::setprecision(5)<<m_r1min
503 <<" |"<<endmsg;
504 out<<"| min radius second SP(3) | "
505 <<std::setw(12)<<std::setprecision(5)<<m_r2min
506 <<" |"<<endmsg;
507 out<<"| min radius last SP(3) | "
508 <<std::setw(12)<<std::setprecision(5)<<m_r3min
509 <<" |"<<endmsg;
510 out<<"| max radius first SP(3) | "
511 <<std::setw(12)<<std::setprecision(4)<<m_r1max
512 <<" |"<<endmsg;
513 out<<"| max radius second SP(3) | "
514 <<std::setw(12)<<std::setprecision(5)<<m_r2max
515 <<" |"<<endmsg;
516 out<<"| max radius last SP(3) | "
517 <<std::setw(12)<<std::setprecision(5)<<m_r3max
518 <<" |"<<endmsg;
519 out<<"| min radius first SP(2) | "
520 <<std::setw(12)<<std::setprecision(5)<<m_r1minv
521 <<" |"<<endmsg;
522 out<<"| min radius second SP(2) | "
523 <<std::setw(12)<<std::setprecision(5)<<m_r2minv
524 <<" |"<<endmsg;
525 out<<"| max radius first SP(2) | "
526 <<std::setw(12)<<std::setprecision(5)<<m_r1maxv
527 <<" |"<<endmsg;
528 out<<"| max radius second SP(2) | "
529 <<std::setw(12)<<std::setprecision(5)<<m_r2maxv
530 <<" |"<<endmsg;
531 out<<"| min space points dR | "
532 <<std::setw(12)<<std::setprecision(5)<<m_drmin
533 <<" |"<<endmsg;
534 out<<"| max space points dR | "
535 <<std::setw(12)<<std::setprecision(5)<<m_drmax
536 <<" |"<<endmsg;
537 out<<"| max dZ impact | "
538 <<std::setw(12)<<std::setprecision(5)<<m_dzver
539 <<" |"<<endmsg;
540 out<<"| max dZ/dR impact | "
541 <<std::setw(12)<<std::setprecision(5)<<m_dzdrver
542 <<" |"<<endmsg;
543 out<<"| max impact | "
544 <<std::setw(12)<<std::setprecision(5)<<m_diver
545 <<" |"<<endmsg;
546 out<<"| max impact pps | "
547 <<std::setw(12)<<std::setprecision(5)<<m_diverpps
548 <<" |"<<endmsg;
549 out<<"| max impact sss | "
550 <<std::setw(12)<<std::setprecision(5)<<m_diversss
551 <<" |"<<endmsg;
552 out<<"|---------------------------------------------------------------------|"
553 <<endmsg;
554 out<<"| Beam X center | "
555 <<std::setw(12)<<std::setprecision(5)<<data.xbeam[0]
556 <<" |"<<endmsg;
557 out<<"| Beam Y center | "
558 <<std::setw(12)<<std::setprecision(5)<<data.ybeam[0]
559 <<" |"<<endmsg;
560 out<<"| Beam Z center | "
561 <<std::setw(12)<<std::setprecision(5)<<data.zbeam[0]
562 <<" |"<<endmsg;
563 out<<"| Beam X-axis direction | "
564 <<std::setw(12)<<std::setprecision(5)<<data.xbeam[1]
565 <<std::setw(12)<<std::setprecision(5)<<data.xbeam[2]
566 <<std::setw(12)<<std::setprecision(5)<<data.xbeam[3]
567 <<" |"<<endmsg;
568 out<<"| Beam Y-axis direction | "
569 <<std::setw(12)<<std::setprecision(5)<<data.ybeam[1]
570 <<std::setw(12)<<std::setprecision(5)<<data.ybeam[2]
571 <<std::setw(12)<<std::setprecision(5)<<data.ybeam[3]
572 <<" |"<<endmsg;
573 out<<"| Beam Z-axis direction | "
574 <<std::setw(12)<<std::setprecision(5)<<data.zbeam[1]
575 <<std::setw(12)<<std::setprecision(5)<<data.zbeam[2]
576 <<std::setw(12)<<std::setprecision(5)<<data.zbeam[3]
577 <<" |"<<endmsg;
578 out<<"|---------------------------------------------------------------------|"
579 <<endmsg;
580 return out;
581}
582
584// Dumps event information into the MsgStream
586
588{
589 out<<"|---------------------------------------------------------------------|"
590 <<endmsg;
591 out<<"| data.ns | "
592 <<std::setw(12)<<data.ns
593 <<" |"<<endmsg;
594 out<<"| data.nsaz | "
595 <<std::setw(12)<<data.nsaz
596 <<" |"<<endmsg;
597 out<<"| data.nsazv | "
598 <<std::setw(12)<<data.nsazv
599 <<" |"<<endmsg;
600 out<<"| seeds | "
601 <<std::setw(12)<<data.l_seeds.size()
602 <<" |"<<endmsg;
603 out<<"|---------------------------------------------------------------------|"
604 <<endmsg;
605 return out;
606}
607
609// Find next set space points
611
613{
614 if (data.endlist) return;
615
616 data.i_seede = data.l_seeds.begin();
617
618 if (data.mode==0 || data.mode==1) production2Sp(data);
619 else if (data.mode==2 || data.mode==3) production3Sp(data);
620 else if (data.mode==5 || data.mode==6) production3Sp(data);
621
622 data.i_seed = data.l_seeds.begin();
623 data.seed = data.mapSeeds.begin();
624 data.seede = data.mapSeeds.end ();
625 ++data.nlist;
626}
627
629// New and old list vertices comparison
631
632bool InDet::SiSpacePointsSeedMaker_Trigger::newVertices(EventData& data, const std::list<Trk::Vertex>& lV) const
633{
634 unsigned int s1 = data.l_vertex.size();
635 unsigned int s2 = lV.size();
636
637 if (s1==0 && s2==0) return false;
638
639 data.l_vertex.clear();
640
641 for (const Trk::Vertex& v : lV) {
642 data.l_vertex.insert(static_cast<float>(v.position().z()));
643 if (data.l_vertex.size() >= m_maxNumberVertices) break;
644 }
645 return false;
646}
647
649// Initiate frame work for seed generator
651
653{
654 m_ptmin = std::abs(m_ptmin);
655 if (m_ptmin < 100.) m_ptmin = 100.;
656 m_etamax = std::abs(m_etamax);
657 m_dzdrmax = 1.f/std::tan(2.f*std::atan(exp(-m_etamax)));
660 m_COF = 134*.05*9.;
661 m_ipt = 1.f/std::abs(.9f*m_ptmin);
663
664 // Build radius sorted containers
665 //
666 m_r_size = static_cast<int>((m_r_rmax+.1f)/m_r_rstep);
667
668 // Build radius-azimuthal sorted containers
669 //
670 constexpr float pi2 = 2.*M_PI;
671 const int NFmax = SizeRF;
672 const float sFmax = static_cast<float>(NFmax)/pi2;
673 const float sFmin = 100./60.;
674
675 m_sF = m_ptmin /60.;
676 if (m_sF > sFmax ) m_sF = sFmax;
677 else if (m_sF < sFmin) m_sF = sFmin;
678 m_fNmax = static_cast<int>(pi2*m_sF);
679 if (m_fNmax >=NFmax) m_fNmax = NFmax-1;
680
681 // Build radius-azimuthal-Z sorted containers for Z-vertices
682 //
683 const int NFtmax = 100;
684 const float sFvmax = static_cast<float>(NFtmax)/pi2;
685 m_sFv = m_ptmin/120.f;
686 if (m_sFv>sFvmax) m_sFv = sFvmax;
687 m_fvNmax = static_cast<int>(pi2*m_sFv);
688 if (m_fvNmax>=NFtmax) m_fvNmax = NFtmax-1;
689
690 // Build maps for radius-azimuthal-Z sorted collections
691 //
692 for (int f=0; f<=m_fNmax; ++f) {
693 int fb = f-1;
694 if (fb<0) fb=m_fNmax;
695 int ft = f+1;
696 if (ft>m_fNmax) ft=0;
697
698 // For each azimuthal region loop through all Z regions
699 //
700 for (int z=0; z<SizeZ; ++z) {
701 int a = f *SizeZ+z;
702 int b = fb*SizeZ+z;
703 int c = ft*SizeZ+z;
704 m_rfz_b [a] = 3; m_rfz_t [a] = 3;
705 m_rfz_ib[a][0] = a; m_rfz_it[a][0] = a;
706 m_rfz_ib[a][1] = b; m_rfz_it[a][1] = b;
707 m_rfz_ib[a][2] = c; m_rfz_it[a][2] = c;
708 if (z==5) {
709
710 m_rfz_t [a] = 9;
711 m_rfz_it[a][3] = a+1;
712 m_rfz_it[a][4] = b+1;
713 m_rfz_it[a][5] = c+1;
714 m_rfz_it[a][6] = a-1;
715 m_rfz_it[a][7] = b-1;
716 m_rfz_it[a][8] = c-1;
717 } else if (z> 5) {
718
719 m_rfz_b [a] = 6;
720 m_rfz_ib[a][3] = a-1;
721 m_rfz_ib[a][4] = b-1;
722 m_rfz_ib[a][5] = c-1;
723
724 if (z<10) {
725
726 m_rfz_t [a] = 6;
727 m_rfz_it[a][3] = a+1;
728 m_rfz_it[a][4] = b+1;
729 m_rfz_it[a][5] = c+1;
730 }
731 } else {
732
733 m_rfz_b [a] = 6;
734 m_rfz_ib[a][3] = a+1;
735 m_rfz_ib[a][4] = b+1;
736 m_rfz_ib[a][5] = c+1;
737
738 if (z>0) {
739
740 m_rfz_t [a] = 6;
741 m_rfz_it[a][3] = a-1;
742 m_rfz_it[a][4] = b-1;
743 m_rfz_it[a][5] = c-1;
744 }
745 }
746
747 if (z==3) {
748 m_rfz_b[a] = 9;
749 m_rfz_ib[a][6] = a+2;
750 m_rfz_ib[a][7] = b+2;
751 m_rfz_ib[a][8] = c+2;
752 } else if (z==7) {
753 m_rfz_b[a] = 9;
754 m_rfz_ib[a][6] = a-2;
755 m_rfz_ib[a][7] = b-2;
756 m_rfz_ib[a][8] = c-2;
757 }
758 }
759 }
760
761 // Build maps for radius-azimuthal-Z sorted collections for Z
762 //
763 for (int f=0; f<=m_fvNmax; ++f) {
764
765 int fb = f-1;
766 if (fb<0) fb=m_fvNmax;
767 int ft = f+1;
768 if (ft>m_fvNmax) ft=0;
769
770 // For each azimuthal region loop through central Z regions
771 //
772 for (int z=0; z<SizeZV; ++z) {
773
774 int a = f *SizeZV+z;
775 int b = fb*SizeZV+z;
776 int c = ft*SizeZV+z;
777 m_rfzv_n[a] = 3;
778 m_rfzv_i[a][0] = a;
779 m_rfzv_i[a][1] = b;
780 m_rfzv_i[a][2] = c;
781 if (z>1) {
782 m_rfzv_n[a] = 6;
783 m_rfzv_i[a][3] = a-1;
784 m_rfzv_i[a][4] = b-1;
785 m_rfzv_i[a][5] = c-1;
786 } else if (z<1) {
787 m_rfzv_n[a] = 6;
788 m_rfzv_i[a][3] = a+1;
789 m_rfzv_i[a][4] = b+1;
790 m_rfzv_i[a][5] = c+1;
791 }
792 }
793 }
794}
795
797// Initiate beam frame work for seed generator
799
801{
803
804 const Amg::Vector3D& cb = beamSpotHandle->beamPos();
805 double tx = std::tan(beamSpotHandle->beamTilt(0));
806 double ty = std::tan(beamSpotHandle->beamTilt(1));
807
808 double ph = std::atan2(ty,tx);
809 double th = std::acos(1./std::sqrt(1.+tx*tx+ty*ty));
810 double sint = sin(th);
811 double cost = cos(th);
812 double sinp = sin(ph);
813 double cosp = cos(ph);
814
815 data.xbeam[0] = static_cast<float>(cb.x());
816 data.xbeam[1] = static_cast<float>(cost*cosp*cosp+sinp*sinp);
817 data.xbeam[2] = static_cast<float>(cost*sinp*cosp-sinp*cosp);
818 data.xbeam[3] = -static_cast<float>(sint*cosp );
819
820 data.ybeam[0] = static_cast<float>(cb.y());
821 data.ybeam[1] = static_cast<float>(cost*cosp*sinp-sinp*cosp);
822 data.ybeam[2] = static_cast<float>(cost*sinp*sinp+cosp*cosp);
823 data.ybeam[3] = -static_cast<float>(sint*sinp );
824
825 data.zbeam[0] = static_cast<float>(cb.z());
826 data.zbeam[1] = static_cast<float>(sint*cosp);
827 data.zbeam[2] = static_cast<float>(sint*sinp);
828 data.zbeam[3] = static_cast<float>(cost);
829}
830
832// Initiate beam frame work for seed generator
834
836(EventData& data, const Trk::SpacePoint*const& sp, float* r)
837{
838 r[0] = static_cast<float>(sp->globalPosition().x())-data.xbeam[0];
839 r[1] = static_cast<float>(sp->globalPosition().y())-data.ybeam[0];
840 r[2] = static_cast<float>(sp->globalPosition().z())-data.zbeam[0];
841}
842
844// Initiate space points seed maker
846
848{
849 constexpr float pi2 = 2.*M_PI;
850 std::vector<InDet::SiSpacePointForSeed*>::iterator r;
851
852 for (int i=0; i!= m_r_size; ++i) {
853
854 if (!data.r_map[i]) continue;
855 r = data.r_Sorted[i].begin();
856
857 while (r!=data.r_Sorted[i].end()) {
858
859 // Azimuthal angle sort
860 //
861 float F = (*r)->phi();
862 if (F<0.) F+=pi2;
863
864 int f = static_cast<int>(F*m_sF);
865 if (f < 0) f = m_fNmax;
866 else if (f > m_fNmax) f = 0;
867
868 int z;
869 float Z = (*r)->z();
870
871 // Azimuthal angle and Z-coordinate sort
872 //
873 if (Z>0.) {Z < 250.? z=5 : z=6;}
874 else {Z >-250.? z=5 : z=4;}
875 int n = f*SizeZ+z;
876 ++data.nsaz;
877 data.rfz_Sorted[n].push_back(*r);
878 if (!data.rfz_map[n]++) data.rfz_index[data.nrfz++] = n;
879
880 if ((*r)->spacepoint->clusterList().second == 0 && z>=3 && z<=7) {
881 z<=4 ? z=0 : z>=6 ? z=2 : z=1;
882
883 // Azimuthal angle and Z-coordinate sort for fast vertex search
884 //
885 f = static_cast<int>(F*m_sFv);
886 if (f < 0) f+=m_fvNmax;
887 else if (f > m_fvNmax) f-=m_fvNmax;
888
889 n = f*3+z;
890 ++data.nsazv;
891 data.rfzv_Sorted[n].push_back(*r);
892 if (!data.rfzv_map[n]++) data.rfzv_index[data.nrfzv++] = n;
893 }
894 data.r_Sorted[i].erase(r++);
895 }
896 data.r_map[i] = 0;
897 }
898 data.nr = 0;
899 data.state = 0;
900}
901
903// Erase space point information
905
907{
908 for (int i=0; i!=data.nr; ++i) {
909 int n = data.r_index[i];
910 data.r_map[n] = 0;
911 data.r_Sorted[n].clear();
912 }
913
914 for (int i=0; i!=data.nrfz; ++i) {
915 int n = data.rfz_index[i];
916 data.rfz_map[n] = 0;
917 data.rfz_Sorted[n].clear();
918 }
919
920 for (int i=0; i!=data.nrfzv; ++i) {
921 int n = data.rfzv_index[i];
922 data.rfzv_map[n] = 0;
923 data.rfzv_Sorted[n].clear();
924 }
925 data.state = 0;
926 data.ns = 0;
927 data.nsaz = 0;
928 data.nsazv = 0;
929 data.nr = 0;
930 data.nrfz = 0;
931 data.nrfzv = 0;
932}
933
935// 2 space points seeds production
937
939{
940 if (data.nsazv<2) return;
941
942 std::vector<InDet::SiSpacePointForSeed*>::iterator r0,r0e,r,re;
943 int nseed = 0;
944
945 // Loop thorugh all azimuthal regions
946 //
947 for (int f=data.fvNmin; f<=m_fvNmax; ++f) {
948
949 // For each azimuthal region loop through Z regions
950 //
951 int z = 0;
952 if (!data.endlist) z = data.zMin;
953 for (; z!=3; ++z) {
954
955 int a = f*3+z;
956 if (!data.rfzv_map[a]) continue;
957 r0 = data.rfzv_Sorted[a].begin();
958 r0e = data.rfzv_Sorted[a].end ();
959
960 if (!data.endlist) {
961 r0 = data.rMin;
962 data.endlist = true;
963 }
964
965 // Loop through trigger space points
966 //
967 for (; r0!=r0e; ++r0) {
968
969 float X = (*r0)->x();
970 float Y = (*r0)->y();
971 float R = (*r0)->radius();
972 if (R<m_r2minv) continue;
973 if (R>m_r2maxv) break;
974 float Z = (*r0)->z();
975 float ax = X/R;
976 float ay = Y/R;
977
978 // Bottom links production
979 //
980 int NB = m_rfzv_n[a];
981 for (int i=0; i!=NB; ++i) {
982
983 int an = m_rfzv_i[a][i];
984 if (!data.rfzv_map[an]) continue;
985
986 r = data.rfzv_Sorted[an].begin();
987 re = data.rfzv_Sorted[an].end ();
988
989 for (; r!=re; ++r) {
990
991 float Rb =(*r)->radius();
992 if (Rb<m_r1minv) continue;
993 if (Rb>m_r1maxv) break;
994 float dR = R-Rb;
995 if (dR<m_drminv) break;
996 if (dR>m_drmax) continue;
997 float dZ = Z-(*r)->z();
998 float Tz = dZ/dR;
999 if (Tz<m_dzdrmin || Tz>m_dzdrmax) continue;
1000 float Zo = Z-R*Tz;
1001
1002 // Comparison with vertices Z coordinates
1003 //
1004 if (!isZCompatible(data, Zo, Rb, Tz)) continue;
1005
1006 // Momentum cut
1007 //
1008 float dx =(*r)->x()-X;
1009 float dy =(*r)->y()-Y;
1010 float x = dx*ax+dy*ay;
1011 float y =-dx*ay+dy*ax;
1012 float xy = x*x+y*y; if (xy == 0.) continue;
1013 float r2 = 1./xy;
1014 float Ut = x*r2;
1015 float Vt = y*r2;
1016 float UR = Ut*R+1.; if (UR == 0.) continue;
1017 float A = Vt*R/UR;
1018 float B = Vt-A*Ut;
1019 if (std::abs(B*data.K) > m_ipt*std::sqrt(1.f+A*A)) continue;
1020 ++nseed;
1021 newSeed(data, (*r)->spacepoint, (*r0)->spacepoint, Zo);
1022 }
1023 }
1024 if (nseed < m_maxsize) continue;
1025 data.endlist=false;
1026 data.rMin = (++r0);
1027 data.fvNmin=f;
1028 data.zMin=z;
1029 return;
1030 }
1031 }
1032 }
1033 data.endlist = true;
1034}
1035
1037// Production 3 space points seeds
1039
1041{
1042 data.mapSeeds.clear();
1043 if (data.nsaz<3) return;
1044
1045 const int ZI[SizeZ]= {5,6,7,8,9,10,4,3,2,1,0};
1046 std::vector<InDet::SiSpacePointForSeed*>::iterator rt[9],rte[9],rb[9],rbe[9];
1047 int nseed = 0;
1048
1049 // Loop thorugh all azimuthal regions
1050 //
1051 for (int f=data.fNmin; f<=m_fNmax; ++f) {
1052
1053 // For each azimuthal region loop through all Z regions
1054 //
1055 int z = 0;
1056 if (!data.endlist) z = data.zMin;
1057
1058 for (; z<SizeZ; ++z) {
1059
1060 int a = f *SizeZ+ZI[z];
1061 if (!data.rfz_map[a]) continue;
1062 int NB = 0, NT = 0;
1063 for (int i=0; i!=m_rfz_b[a]; ++i) {
1064
1065 int an = m_rfz_ib[a][i];
1066 if (!data.rfz_map[an]) continue;
1067 rb [NB] = data.rfz_Sorted[an].begin();
1068 rbe[NB++] = data.rfz_Sorted[an].end();
1069 }
1070 for (int i=0; i!=m_rfz_t[a]; ++i) {
1071
1072 int an = m_rfz_it[a][i];
1073 if (!data.rfz_map[an]) continue;
1074 rt [NT] = data.rfz_Sorted[an].begin();
1075 rte[NT++] = data.rfz_Sorted[an].end();
1076 }
1077 if (!data.trigger) production3Sp (data, rb, rbe, rt, rte, NB, NT, nseed);
1078 else production3SpTrigger(data, rb, rbe, rt, rte, NB, NT, nseed);
1079 if (!data.endlist) {
1080 data.fNmin = f;
1081 data.zMin = z;
1082 return;
1083 }
1084 }
1085 }
1086 data.endlist = true;
1087}
1088
1090// Production 3 space points seeds for full scan
1092
1094(EventData& data,
1095 std::vector<InDet::SiSpacePointForSeed*>::iterator* rb ,
1096 std::vector<InDet::SiSpacePointForSeed*>::iterator* rbe,
1097 std::vector<InDet::SiSpacePointForSeed*>::iterator* rt ,
1098 std::vector<InDet::SiSpacePointForSeed*>::iterator* rte,
1099 int NB, int NT, int& nseed) const
1100{
1101 std::vector<InDet::SiSpacePointForSeed*>::iterator r0=rb[0],r;
1102 if (!data.endlist) {
1103 r0 = data.rMin;
1104 data.endlist = true;
1105 }
1106
1107 // Loop through all trigger space points
1108 //
1109 for (; r0!=rbe[0]; ++r0) {
1110
1111 data.nOneSeeds = 0;
1112 data.mapOneSeeds.erase(data.mapOneSeeds.begin(), data.mapOneSeeds.end());
1113
1114 float R = (*r0)->radius();
1115 if (R<m_r2min) continue;
1116 if (R>m_r2max) break;
1117
1118 const Trk::SpacePoint* SP0 = (*r0)->spacepoint;
1119
1120 bool pix = true;
1121 if (SP0->clusterList().second) pix = false;
1122 const Trk::Surface* sur0 = (*r0)->sur();
1123 float X = (*r0)->x();
1124 float Y = (*r0)->y();
1125 float Z = (*r0)->z();
1126 int Nb = 0;
1127
1128 // Bottom links production
1129 //
1130 for (int i=0; i!=NB; ++i) {
1131
1132 for (r=rb[i]; r!=rbe[i]; ++r) {
1133
1134 float Rb =(*r)->radius();
1135 if (Rb<m_r1min) {
1136 rb[i]=r;
1137 continue;
1138 }
1139 if (Rb>m_r1max) break;
1140
1141 float dR = R-Rb;
1142 if (dR<m_drmin) break;
1143
1144 if (dR > m_drmax || (*r)->sur()==sur0) continue;
1145
1146 if ( !pix && !(*r)->spacepoint->clusterList().second) continue;
1147
1148 float Tz = (Z-(*r)->z())/dR;
1149
1150 if (Tz < m_dzdrmin || Tz > m_dzdrmax) continue;
1151
1152 // Comparison with vertices Z coordinates
1153 //
1154 float Zo = Z-R*Tz;
1155 if (!isZCompatible(data, Zo, Rb, Tz)) continue;
1156 data.SP[Nb] = (*r);
1157 if (++Nb==m_maxsizeSP) goto breakb;
1158 }
1159 }
1160 breakb:
1161 if (!Nb || Nb==m_maxsizeSP) continue;
1162 int Nt = Nb;
1163
1164 // Top links production
1165 //
1166 for (int i=0; i!=NT; ++i) {
1167
1168 for (r=rt[i]; r!=rte[i]; ++r) {
1169
1170 float Rt =(*r)->radius();
1171 float dR = Rt-R;
1172 if (dR<m_drmin || Rt<m_r3min) {
1173 rt[i]=r;
1174 continue;
1175 }
1176 if (Rt>m_r3max || dR>m_drmax) break;
1177
1178 if ( (*r)->sur()==sur0) continue;
1179
1180 float Tz = ((*r)->z()-Z)/dR;
1181
1182 if (Tz < m_dzdrmin || Tz > m_dzdrmax) continue;
1183
1184 // Comparison with vertices Z coordinates
1185 //
1186 float Zo = Z-R*Tz;
1187 if (!isZCompatible(data, Zo, Rt, Tz)) continue;
1188 data.SP[Nt] = (*r);
1189 if (++Nt==m_maxsizeSP) goto breakt;
1190 }
1191 }
1192
1193 breakt:
1194 if (!(Nt-Nb)) continue;
1195
1196 float covr0 = (*r0)->covr ();
1197 float covz0 = (*r0)->covz ();
1198
1199 float ax = X/R;
1200 float ay = Y/R;
1201
1202 for (int i=0; i!=Nt; ++i) {
1203
1205
1206 float dx = sp->x()-X;
1207 float dy = sp->y()-Y;
1208 float dz = sp->z()-Z;
1209 float x = dx*ax+dy*ay;
1210 float y =-dx*ay+dy*ax;
1211 float r2 = 1.f/(x*x+y*y);
1212 float dr = std::sqrt(r2);
1213 float tz = dz*dr;
1214 if (i < Nb) tz = -tz;
1215
1216 data.Tz[i] = tz;
1217 data.Zo[i] = Z-R*tz;
1218 data.R [i] = dr;
1219 data.U [i] = x*r2;
1220 data.V [i] = y*r2;
1221 data.Er[i] = (covz0+sp->covz()+tz*tz*(covr0+sp->covr()))*r2;
1222 }
1223
1224 float imc = m_diver;
1225 float imcs = m_diverpps;
1226 float ipt2 = m_ipt2;
1227 float K = data.K;
1228 float K2 = K*K;
1229 float COF = m_COF;
1230 float ipt2K = ipt2/K2;
1231 float ipt2C = ipt2*COF;
1232 float COFK = COF*K2;
1233 covr0 *= 2.f;
1234 covz0 *= 2.f;
1235
1236 // Three space points comparison
1237 //
1238 for (int b=0; b!=Nb; ++b) {
1239
1240 const Trk::SpacePoint* SPb = data.SP[b]->spacepoint;
1241
1242 float Zob = data.Zo[b];
1243 float Tzb = data.Tz[b];
1244 float Rb2r = data.R [b]*covr0;
1245 float Rb2z = data.R [b]*covz0;
1246 float Erb = data.Er[b];
1247 float Vb = data.V [b];
1248 float Ub = data.U [b];
1249 float Tzb2 = (1.f+Tzb*Tzb);
1250 float CSA = Tzb2*COFK;
1251 float ICSA = Tzb2*ipt2C;
1252 float dZ = dZVertexMin(data, Zob);
1253 float Iz = (dZ*dZ)/Tzb2;
1254
1255 for (int t=Nb; t!=Nt; ++t) {
1256
1257 float Ts = .5*(Tzb+data.Tz[t]);
1258 float dt = Tzb-data.Tz[t];
1259 float dT = dt*dt-Erb-data.Er[t]-data.R[t]*(Ts*Ts*Rb2r+Rb2z);
1260 if ( dT > ICSA) continue;
1261 float dU = data.U[t]-Ub;
1262 if (dU == 0.) continue;
1263 float A = (data.V[t]-Vb)/dU;
1264 float S2 = 1.+A*A;
1265 float B = Vb-A*Ub;
1266 float B2 = B*B;
1267 if (B2 > ipt2K*S2 || dT*S2 > B2*CSA) continue;
1268 float Im = std::abs((A-B*R)*R);
1269
1270 if (Im > imc ) continue;
1271
1272 if (pix) {
1273 if (data.SP[t]->spacepoint->clusterList().second && Im > imcs) continue;
1274 }
1275 if (Im > m_diversss) continue;
1276 Im = Im*Im+Iz;
1277 newOneSeed(data, SPb, SP0, data.SP[t]->spacepoint, Zob, Im);
1278 }
1279 }
1280 nseed += data.mapOneSeeds.size();
1281 fillSeeds(data);
1282 if (nseed>=m_maxsize) {
1283 data.endlist=false;
1284 ++r0;
1285 data.rMin = r0;
1286 return;
1287 }
1288 }
1289}
1290
1292// Production 3 space points seeds in ROI
1294
1296(EventData& data,
1297 std::vector<InDet::SiSpacePointForSeed*>::iterator* rb ,
1298 std::vector<InDet::SiSpacePointForSeed*>::iterator* rbe,
1299 std::vector<InDet::SiSpacePointForSeed*>::iterator* rt ,
1300 std::vector<InDet::SiSpacePointForSeed*>::iterator* rte,
1301 int NB, int NT, int& nseed) const
1302{
1303 constexpr float pi2 = 2.*M_PI;
1304
1305 std::vector<InDet::SiSpacePointForSeed*>::iterator r0=rb[0],r;
1306 if (!data.endlist) {
1307 r0 = data.rMin;
1308 data.endlist = true;
1309 }
1310
1311 // Loop through all trigger space points
1312 //
1313 for (; r0!=rbe[0]; ++r0) {
1314
1315 data.nOneSeeds = 0;
1316 data.mapOneSeeds.erase(data.mapOneSeeds.begin(), data.mapOneSeeds.end());
1317
1318 float R = (*r0)->radius();
1319 if (R<m_r2min) continue;
1320 if (R>m_r2max) break;
1321
1322 const Trk::SpacePoint* SP0 = (*r0)->spacepoint;
1323
1324 bool pix = true;
1325 if (SP0->clusterList().second) pix = false;
1326 const Trk::Surface* sur0 = (*r0)->sur();
1327 float X = (*r0)->x();
1328 float Y = (*r0)->y();
1329 float Z = (*r0)->z();
1330 int Nb = 0;
1331
1332 // Bottom links production
1333 //
1334 for (int i=0; i!=NB; ++i) {
1335
1336 for (r=rb[i]; r!=rbe[i]; ++r) {
1337
1338 float Rb =(*r)->radius();
1339 if (Rb<m_r1min) {
1340 rb[i]=r;
1341 continue;
1342 }
1343 if (Rb>m_r1max) break;
1344
1345 float dR = R-Rb;
1346 if (dR<m_drmin) break;
1347
1348 if (dR > m_drmax || (*r)->sur()==sur0) continue;
1349
1350 if ( !pix && !(*r)->spacepoint->clusterList().second) continue;
1351
1352 // Comparison with bottom and top Z
1353 //
1354 float Tz = (Z-(*r)->z())/dR;
1355 float Zo = Z-R*Tz;
1356 if (Zo < data.zminB || Zo > data.zmaxB) continue;
1357 float Zu = Z+(550.f-R)*Tz;
1358 if (Zu < data.zminU || Zu > data.zmaxU) continue;
1359 data.SP[Nb] = (*r);
1360 if (++Nb==m_maxsizeSP) goto breakb;
1361 }
1362 }
1363 breakb:
1364 if (!Nb || Nb==m_maxsizeSP) continue;
1365 int Nt = Nb;
1366
1367 // Top links production
1368 //
1369 for (int i=0; i!=NT; ++i) {
1370
1371 for (r=rt[i]; r!=rte[i]; ++r) {
1372
1373 float Rt =(*r)->radius();
1374 float dR = Rt-R;
1375 if (dR<m_drmin || Rt<m_r3min) {
1376 rt[i]=r;
1377 continue;
1378 }
1379 if (Rt>m_r3max || dR>m_drmax) break;
1380
1381 if ( (*r)->sur()==sur0) continue;
1382
1383 // Comparison with bottom and top Z
1384 //
1385 float Tz = ((*r)->z()-Z)/dR;
1386 float Zo = Z-R*Tz;
1387 if (Zo < data.zminB || Zo > data.zmaxB) continue;
1388 float Zu = Z+(550.f-R)*Tz;
1389 if (Zu < data.zminU || Zu > data.zmaxU) continue;
1390 data.SP[Nt] = (*r);
1391 if (++Nt==m_maxsizeSP) goto breakt;
1392 }
1393 }
1394
1395 breakt:
1396 if (!(Nt-Nb)) continue;
1397
1398 float covr0 = (*r0)->covr ();
1399 float covz0 = (*r0)->covz ();
1400
1401 float ax = X/R;
1402 float ay = Y/R;
1403
1404 for (int i=0; i!=Nt; ++i) {
1405
1407
1408 float dx = sp->x()-X;
1409 float dy = sp->y()-Y;
1410 float dz = sp->z()-Z;
1411 float x = dx*ax+dy*ay;
1412 float y =-dx*ay+dy*ax;
1413 float r2 = 1.f/(x*x+y*y);
1414 float dr = std::sqrt(r2);
1415 float tz = dz*dr;
1416 if (i < Nb) tz = -tz;
1417
1418 data.Tz[i] = tz;
1419 data.Zo[i] = Z-R*tz;
1420 data.R [i] = dr;
1421 data.U [i] = x*r2;
1422 data.V [i] = y*r2;
1423 data.Er[i] = (covz0+sp->covz()+tz*tz*(covr0+sp->covr()))*r2;
1424 }
1425
1426 float imc = m_diver;
1427 float imcs = m_diverpps;
1428 float ipt2 = m_ipt2;
1429 float K = data.K;
1430 float K2 = K*K;
1431 float COF = m_COF;
1432 float ipt2K = ipt2/K2;
1433 float ipt2C = ipt2*COF;
1434 float COFK = COF*K2;
1435 covr0 *= 2.f;
1436 covz0 *= 2.f;
1437
1438 // Three space points comparison
1439 //
1440 for (int b=0; b!=Nb; ++b) {
1441
1442 const Trk::SpacePoint* SPb = data.SP[b]->spacepoint;
1443
1444 float Zob = data.Zo[b];
1445 float Tzb = data.Tz[b];
1446 float Rb2r = data.R [b]*covr0;
1447 float Rb2z = data.R [b]*covz0;
1448 float Erb = data.Er[b];
1449 float Vb = data.V [b];
1450 float Ub = data.U [b];
1451 float Tzb2 = (1.f+Tzb*Tzb);
1452 float CSA = Tzb2*COFK;
1453 float ICSA = Tzb2*ipt2C;
1454 float dZ = dZVertexMin(data, Zob);
1455 float Iz = (dZ*dZ)/Tzb2;
1456
1457 for (int t=Nb; t!=Nt; ++t) {
1458
1459 float Ts = .5f*(Tzb+data.Tz[t]);
1460 float dt = Tzb-data.Tz[t];
1461 float dT = dt*dt-Erb-data.Er[t]-data.R[t]*(Ts*Ts*Rb2r+Rb2z);
1462 if ( dT > ICSA) continue;
1463 float dU = data.U[t]-Ub;
1464 if (dU == 0.) continue;
1465 float A = (data.V[t]-Vb)/dU;
1466 float S2 = 1.f+A*A;
1467 float B = Vb-A*Ub;
1468 float B2 = B*B;
1469 if (B2 > ipt2K*S2 || dT*S2 > B2*CSA) continue;
1470 float Im = std::abs((A-B*R)*R);
1471
1472 if (Im > imc ) continue;
1473 if (pix) {
1474 if (data.SP[t]->spacepoint->clusterList().second && Im > imcs) continue;
1475 }
1476
1477 if (Im > m_diversss) continue;
1478
1479 // Azimuthal angle test
1480 //
1481 float y = 1.;
1482 float x = 2.f*B*R-A;
1483 float df = std::abs(std::atan2(ay*y-ax*x,ax*y+ay*x)-data.ftrig);
1484 if (df > M_PI) df=pi2-df;
1485 if (df > data.ftrigW) continue;
1486 Im = Im*Im+Iz;
1487 newOneSeed(data, SPb, SP0, data.SP[t]->spacepoint, Zob, Im);
1488 }
1489 }
1490 nseed += data.mapOneSeeds.size();
1491 fillSeeds(data);
1492 if (nseed>=m_maxsize) {
1493 data.endlist=false;
1494 ++r0;
1495 data.rMin = r0;
1496 return;
1497 }
1498 }
1499}
1500
1502// New 3 space points seeds from one space points
1504
1506(EventData& data,
1507 const Trk::SpacePoint*& p1,const Trk::SpacePoint*& p2,
1508 const Trk::SpacePoint*& p3,const float& z,const float& q) const
1509{
1510 if (data.nOneSeeds < m_maxOneSize) {
1511
1512 data.OneSeeds[data.nOneSeeds].erase();
1513 data.OneSeeds[data.nOneSeeds].add(p1);
1514 data.OneSeeds[data.nOneSeeds].add(p2);
1515 data.OneSeeds[data.nOneSeeds].add(p3);
1516 data.OneSeeds[data.nOneSeeds].setZVertex(static_cast<double>(z));
1517 data.mapOneSeeds.insert(std::make_pair(q, &(data.OneSeeds[data.nOneSeeds])));
1518 ++data.nOneSeeds;
1519 } else {
1520 std::multimap<float,InDet::SiSpacePointsSeed*>::reverse_iterator
1521 l = data.mapOneSeeds.rbegin();
1522 if ((*l).first <= q) return;
1523
1524 InDet::SiSpacePointsSeed* s = (*l).second;
1525 s->erase ( );
1526 s->add (p1);
1527 s->add (p2);
1528 s->add (p3);
1529 s->setZVertex(static_cast<double>(z));
1530 std::multimap<float,InDet::SiSpacePointsSeed*>::iterator
1531 i = data.mapOneSeeds.insert(std::make_pair(q,s));
1532
1533 for (++i; i!=data.mapOneSeeds.end(); ++i) {
1534 if ((*i).second==s) {
1535 data.mapOneSeeds.erase(i);
1536 return;
1537 }
1538 }
1539 }
1540}
1541
1543{
1544 if (not data.initialized) initializeEventData(data);
1545
1546 if (data.i_seed==data.i_seede) {
1547 findNext(data);
1548 //cppcheck-suppress identicalInnerCondition
1549 if (data.i_seed==data.i_seede) return nullptr;
1550 }
1551 if (data.mode==0 || data.mode==1) return &(*data.i_seed++);
1552 ++data.i_seed;
1553 return (*data.seed++).second;
1554}
1555
1557(EventData& data, float& Zv, float& R, float& T) const
1558{
1559 if (Zv < m_zmin || Zv > m_zmax) return false;
1560
1561 if (data.l_vertex.empty()) return true;
1562
1563 float dZmin = std::numeric_limits<float>::max();
1564 for (const float& v : data.l_vertex) {
1565 float dZ = std::abs(v-Zv);
1566 if (dZ<dZmin) dZmin=dZ;
1567 }
1568 return dZmin < (m_dzver+m_dzdrver*R)*sqrt(1.+T*T);
1569}
1570
1572{
1573 if (data.l_vertex.empty()) return 0.;
1574
1575 float dZmin = std::numeric_limits<float>::max();
1576 for (const float& v : data.l_vertex) {
1577 float dZ = std::abs(v-Z);
1578 if (dZ<dZmin) dZmin = dZ;
1579 }
1580 return dZmin;
1581}
1582
1584// New space point for seeds
1586
1588(EventData& data, const Trk::SpacePoint*const& sp)
1589{
1590 InDet::SiSpacePointForSeed* sps = nullptr;
1591
1592 float r[3];
1594
1595 if (data.i_spforseed!=data.l_spforseed.end()) {
1596 sps = &(*data.i_spforseed++);
1597 sps->set(sp, r);
1598 } else {
1599 data.l_spforseed.emplace_back(sp, r);
1600 sps = &(data.l_spforseed.back());
1601 data.i_spforseed = data.l_spforseed.end();
1602 }
1603
1604 return sps;
1605}
1606
1608// New 2 space points seeds
1610
1612(EventData& data,
1613 const Trk::SpacePoint*& p1,const Trk::SpacePoint*& p2,
1614 const float& z)
1615{
1616 if (data.i_seede!=data.l_seeds.end()) {
1617 InDet::SiSpacePointsSeed* s = &(*data.i_seede++);
1618 s->erase ( );
1619 s->add (p1);
1620 s->add (p2);
1621 s->setZVertex(static_cast<double>(z));
1622 } else {
1623 data.l_seeds.emplace_back(p1, p2, z);
1624 data.i_seede = data.l_seeds.end();
1625 }
1626}
1627
1629// Fill seeds
1631
1633{
1634 std::multimap<float, InDet::SiSpacePointsSeed*>::iterator
1635 l = data.mapOneSeeds.begin(),
1636 le = data.mapOneSeeds.end ();
1637
1638 for (; l!=le; ++l) {
1639 float q = (*l).first;
1640 InDet::SiSpacePointsSeed* s0 = (*l).second;
1641
1642 if ((*s0->spacePoints().rbegin())->clusterList().second) {
1643 (*s0->spacePoints().begin())->clusterList().second ? q+=1000. : q+=10000.;
1644 }
1645
1646 if (data.i_seede!=data.l_seeds.end()) {
1647 InDet::SiSpacePointsSeed* s = &(*data.i_seede++);
1648 *s = *s0;
1649 data.mapSeeds.insert(std::make_pair(q,s));
1650 } else {
1651 data.l_seeds.emplace_back(*s0);
1652 InDet::SiSpacePointsSeed* s = &(data.l_seeds.back());
1653 data.i_seede = data.l_seeds.end();
1654 data.mapSeeds.insert(std::make_pair(q, s));
1655 }
1656 }
1657}
1658
1660 data.initialize(EventData::ToolType::Trigger,
1663 0, // maxsize not used
1664 m_r_size,
1665 0, // sizeRF not used
1666 SizeRFZ,
1667 SizeRFZV,
1668 false); // checkEta not used
1669}
1670
1673
const boost::regex re(r_e)
#define M_PI
#define endmsg
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_ERROR(x)
char data[hepevt_bytes_allocation_ATLAS]
Definition HepEvt.cxx:11
static Double_t sp
static Double_t s0
static Double_t a
static Double_t sc
#define F(x, y, z)
Definition MD5.cxx:112
struct TBPatternUnitContext S2
#define y
#define x
#define z
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.
void set(const Trk::SpacePoint *, std::span< float const, 3 >)
virtual void newEvent(const EventContext &ctx, EventData &data, int iteration) const override
virtual void writeNtuple(const SiSpacePointsSeed *seed, const Trk::Track *track, int seedType, long eventNumber) const override
static void newSeed(EventData &data, const Trk::SpacePoint *&, const Trk::SpacePoint *&, const float &)
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
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_fieldCondObjInputKey
static void convertToBeamFrameWork(EventData &data, const Trk::SpacePoint *const &, float *)
virtual MsgStream & dump(EventData &data, MsgStream &out) const override
void production3SpTrigger(EventData &data, std::vector< InDet::SiSpacePointForSeed * >::iterator *, std::vector< InDet::SiSpacePointForSeed * >::iterator *, std::vector< InDet::SiSpacePointForSeed * >::iterator *, std::vector< InDet::SiSpacePointForSeed * >::iterator *, int, int, int &) const
bool isZCompatible(EventData &data, float &, float &, float &) const
virtual void newRegion(const EventContext &ctx, EventData &data, const std::vector< IdentifierHash > &vPixel, const std::vector< IdentifierHash > &vSCT) const override
SG::ReadCondHandleKey< InDet::BeamSpotData > m_beamSpotKey
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,...
static float dZVertexMin(EventData &data, float &)
virtual const SiSpacePointsSeed * next(const EventContext &ctx, EventData &data) const override
bool newVertices(EventData &data, const std::list< Trk::Vertex > &) const
SG::ReadHandleKey< SpacePointContainer > m_spacepointsSCT
void newOneSeed(EventData &data, const Trk::SpacePoint *&, const Trk::SpacePoint *&, const Trk::SpacePoint *&, const float &, const float &) const
virtual void find2Sp(EventData &data, const std::list< Trk::Vertex > &lv) const override
with two space points with or without vertex constraint
SG::ReadHandleKey< SpacePointContainer > m_spacepointsPixel
static SiSpacePointForSeed * newSpacePoint(EventData &data, const Trk::SpacePoint *const &)
MsgStream & dumpConditions(EventData &data, MsgStream &out) const
SG::ReadHandleKey< SpacePointOverlapCollection > m_spacepointsOverlap
static MsgStream & dumpEvent(EventData &data, MsgStream &out)
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 std::pair< const PrepRawData *, const PrepRawData * > & clusterList() const
return the pair of cluster pointers by reference
Abstract Base Class for tracking surfaces.
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
int cost(std::vector< std::string > &files, node &n, const std::string &directory="", bool deleteref=false, bool relocate=false)
Definition hcg.cxx:922
Eigen::Matrix< double, 3, 1 > Vector3D
-event-from-file
hold the test vectors and ease the comparison
MsgStream & msg
Definition testRead.cxx:32