ATLAS Offline Software
SiTrajectoryElement_xk.icc
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include <type_traits>
6 
7 ///////////////////////////////////////////////////////////////////
8 // Search clusters compatible with track
9 ///////////////////////////////////////////////////////////////////
10 
11 template<typename T>
12 int InDet::SiTrajectoryElement_xk::searchClustersSub
13 (Trk::PatternTrackParameters& Tp, SiClusterLink_xk* L) {
14  /// Case a): No PRD association in use
15  if (not m_useassoTool) {
16  /// a1): No stereo
17  if (not m_stereo) {
18  if (m_ndf == 2) {
19  return searchClustersWithoutStereoPIX<T>(Tp, L);
20  } else {
21  return searchClustersWithoutStereoSCT<T>(Tp, L);
22  }
23  }
24  /// a2): With stereo
25  else {
26  return searchClustersWithStereo<T>(Tp, L);
27  }
28  }
29  /// Case b) Using PRD association tool
30  else {
31  /// b1): no stereo
32  if (not m_stereo) {
33  if (m_ndf == 2) {
34  return searchClustersWithoutStereoAssPIX<T>(Tp, L,*m_prdToTrackMap);
35  } else {
36  return searchClustersWithoutStereoAssSCT<T>(Tp, L,*m_prdToTrackMap);
37  }
38  }
39  /// b2): with stereo
40  else {
41  return searchClustersWithStereoAss<T>(Tp, L, *m_prdToTrackMap);
42  }
43  }
44 }
45 
46 ///////////////////////////////////////////////////////////////////
47 // Search closest cluster with stereo angle
48 ///////////////////////////////////////////////////////////////////
49 
50 template <typename T>
51 int InDet::SiTrajectoryElement_xk::searchClustersWithStereo
52 (Trk::PatternTrackParameters& Tp,InDet::SiClusterLink_xk* L)
53 {
54  if (m_detstatus<=0) return 0;
55 
56  const AmgSymMatrix(5) & Tp_cov = *Tp.covariance();
57 
58  int nLinksFound = 0;
59  /// predicted local position
60  double P0 = Tp.parameters()[0];
61  double P1 = Tp.parameters()[1];
62  /// covariance for prediction
63  double PV0 = Tp_cov(0, 0); /// X,X
64  double PV1 = Tp_cov(0, 1); /// X,Y
65  double PV2 = Tp_cov(1, 1); /// Y,Y
66  /// chi² cut values
67  double Xc = m_xi2maxlink;
68  double bestX2 = m_xi2maxlink;
69  double Xm = m_xi2max ;
70  int nMaxClusters = m_tools->isITkGeometry() ? m_tools->maxclusters() : 9;
71  int nMaxClustersForBestChi2 = m_tools->isITkGeometry() ? nMaxClusters-1 : nMaxClusters;
72 
73  /// holder for best cluster seen not pasing cuts
74  const InDet::SiCluster* cl = nullptr;
75 
76  T* sibegin = std::any_cast<T>(&m_sibegin);
77  T* siend = std::any_cast<T>(&m_siend);
78  if (sibegin==nullptr or siend==nullptr) return 0;
79 
80  /// loop over all clusters on this element
81  for (T p=*sibegin; p!=*siend; ++p) {
82  const InDet::SiCluster* c = static_cast<const InDet::SiCluster*>(*p);
83  const Amg::Vector2D& M = c->localPosition();
84  const Amg::MatrixX& V = c->localCovariance();
85 
86  double MV0 = V(0,0);
87  double MV1 = V(1,0);
88  double MV2 = V(1,1);
89  /// add the two covariances
90  double v0 = MV0+PV0;
91  double v1 = MV1+PV1;
92  double v2 = MV2+PV2;
93  /// get the distance from the cluster to the predicted location
94  double r0 = M[0]-P0;
95  double r1 = M[1]-P1;
96  /// determinant of the cov matrix
97  double d = v0*v2-v1*v1;
98  if(d<=0.) continue;
99  d=1./d;
100 
101  /// chi² calculation
102  double x = (r0*(r0*v2-r1*v1)+r1*(r1*v0-r0*v1))*d;
103  if(x > Xc) continue; /// if we don't satisfy the minimum, don't bother further
104 
105  /// refine y estimate
106  r1 = fabs(r1+d*((PV1*v2-PV2*v1)*r0+(PV2*v0-PV1*v1)*r1));
107  x -= (r1*r1)/MV2 ;
108  r1 -= m_halflength ;
109 
110  /// update the chi² estimate, and check again the cut
111  if(r1 > 0.){
112  x+=(r1*r1)/PV2;
113  if (x > Xc) continue;
114  }
115  /// if we satisfy the chi2 criterion:
116  if(x < Xm) {
117  /// create a cluster link to this cluster
118  InDet::SiClusterLink_xk l(c,x);
119 
120  /// This misleadingly named method actually is a sorting
121  /// Mechanism: Our current link will take the place of the one
122  /// it replaces in ascending chi², and we will be
123  /// left with a link to the worst chi² one
124  for(int i=0; i!=nLinksFound; ++i) L[i].Comparison(l);
125  /// Now our link points not the the cluster we just processed but the worst
126  /// so far
127 
128  /// If we have less than nm saved, save it in the next slot
129  if(nLinksFound<=nMaxClusters) {
130  L[nLinksFound++]=l;
131  } else {
132  /// otherwise, don't add it, and update the chi2 cut to
133  /// the worst one in our list, to get rid of further
134  /// candidates early on.
135  Xm=L[nMaxClustersForBestChi2].xi2();
136  }
137 
138  /// finally, update the chi2 cut - don't bother with anything worse
139  /// than the current candidate by 6 units in chi2
140  Xc = Xm+6.;
141  } /// done with branch if we satisfy the chi2 criterion
142  /// if we have no links yet and this is the best chi2 so far
143  else if(!nLinksFound && x < bestX2) {
144  /// update best chi2
145  bestX2 = x;
146  /// update max cut
147  Xc = x+6.;
148  /// update best cluster pointer
149  cl = c;
150  }
151  } /// end cluster loop
152  /// If we didn't find any cluster satisfying the chi2 cut,
153  /// but have at least one cluster,
154  /// we write one link to the best (of the bad) candidate
155  if(cl && !nLinksFound) {L[nLinksFound++].Set(cl,bestX2);}
156  /// return number of found links
157  return nLinksFound;
158 }
159 
160 ///////////////////////////////////////////////////////////////////
161 // Search closest cluster without stereo angle for pixels
162 ///////////////////////////////////////////////////////////////////
163 
164 template <typename T>
165 int InDet::SiTrajectoryElement_xk::searchClustersWithoutStereoPIX
166 (Trk::PatternTrackParameters& Tp,InDet::SiClusterLink_xk* L)
167 {
168  if (m_detstatus<=0) return 0;
169 
170  const AmgSymMatrix(5) & Tp_cov = *Tp.covariance();
171  int nLinksFound = 0;
172  /// predicted local position
173  double P0 = Tp.parameters()[0];
174  double P1 = Tp.parameters()[1];
175  /// covariance of prediction
176  double PV0 = Tp_cov(0, 0);
177  double PV1 = Tp_cov(0, 1);
178  double PV2 = Tp_cov(1, 1);
179  /// chi2 cuts
180  double Xc = m_xi2maxlink;
181  double Xm = m_xi2max ;
182  int nMaxClusters = m_tools->isITkGeometry() ? m_tools->maxclusters() : 9;
183  int nMaxClustersForBestChi2 = m_tools->isITkGeometry() ? nMaxClusters-1 : nMaxClusters;
184 
185  /// best cluster seen
186  const InDet::SiCluster* cl = nullptr;
187 
188  T* sibegin = std::any_cast<T>(&m_sibegin);
189  T* siend = std::any_cast<T>(&m_siend);
190  if (sibegin==nullptr or siend==nullptr) return 0;
191 
192  /// loop over all hits
193  for (T p=*sibegin; p!=*siend; ++p) {
194  const InDet::SiCluster* c = static_cast<const InDet::SiCluster*>(*p);
195  const Amg::Vector2D& M = c->localPosition();
196 
197  /// estimate the covariance by the pitch
198  /// Factors 1/sqrt(12) will be added later
199  double MV0 = c->width().phiR();
200  double MV2 = c->width().z ();
201 
202  /// squared distance between prediction and hit location,
203  double r0 = M[0]-P0,
204  r02 = r0*r0;
205  double r1 = M[1]-P1,
206  r12 = r1*r1;
207 
208  /// sum cov on X - using 1/sqrt(12) update of hit cov
209  double v0 = s_oneOverTwelve*(MV0*MV0)+PV0;
210  /// bail out if X alone sufficient to fail chi2
211  if(r02 >(Xc*v0)) continue;
212 
213  /// sum cov on Y - using 1/sqrt(12) update of hit cov
214  double v2 = s_oneOverTwelve*(MV2*MV2)+PV2;
215  /// bail out if Y alone sufficient to fail chi2
216  if(r12 >(Xc*v2)) continue;
217  double v1 = PV1;
218  /// determinant of cov matrix
219  double d = v0*v2-v1*v1;
220  if( d<=0. ) continue;
221  /// chi2
222  double x = (r02*v2+r12*v0-(r0*r1)*(2.*v1))/d;
223 
224  /// skip clusters failing the cut
225  if(x>Xc) continue;
226 
227  /// if we satisfy the tighter requirement:
228  if(x < Xm) {
229  /// set up a cluster link
230  InDet::SiClusterLink_xk l(c,x);
231  /// This misleadingly named method actually is a sorting
232  /// Mechanism: Our current link will take the place of the one
233  /// it replaces in ascending chi², and we will be
234  /// left with a link to the worst chi² one
235  for(int i=0; i!=nLinksFound; ++i) L[i].Comparison(l);
236  /// and place the worst link (remember that l is now a different link, see above)!
237  /// If within capacity, increment counter and add
238  if(nLinksFound<=nMaxClusters)
239  L[nLinksFound++]=l;
240  /// otherwise increment cut to avoid bothering with lower qualities in the future
241  else Xm=L[nMaxClustersForBestChi2].xi2();
242  /// update looser chi2 cut to match the tighter one - we found at least one good cluster
243  Xc = Xm;
244  }
245  /// if we found no clusters yet satisfying the chi2 cut, use this one
246  else if(!nLinksFound) {
247  Xc = x; cl = c;
248  }
249  } /// end of cluster loop
250  /// if we found no cluster satisfying the chi2, we return one link to the best (least bad) we saw
251  if(cl && !nLinksFound) {
252  L[nLinksFound++].Set(cl,Xc);
253  }
254  return nLinksFound;
255 }
256 
257 ///////////////////////////////////////////////////////////////////
258 // Search closest cluster without stereo angle for SCT
259 ///////////////////////////////////////////////////////////////////
260 
261 template <typename T>
262 int InDet::SiTrajectoryElement_xk::searchClustersWithoutStereoSCT
263 (Trk::PatternTrackParameters& Tp,InDet::SiClusterLink_xk* L)
264 {
265  if (m_detstatus<=0) return 0;
266 
267  const AmgSymMatrix(5) & Tp_cov = *Tp.covariance();
268  int nLinksFound = 0;
269  /// predicted local position
270  double P0 = Tp.parameters()[0];
271  double P1 = Tp.parameters()[1];
272  /// and cov of prediction
273  double PV0 = Tp_cov(0, 0);
274  double PV1 = Tp_cov(0, 1);
275  double PV2 = Tp_cov(1, 1);
276  /// chi2 cuts
277  double Xc = m_xi2maxlink;
278  double Xm = m_xi2max ;
279  int nMaxClusters = m_tools->isITkGeometry() ? m_tools->maxclusters() : 9;
280  int nMaxClustersForBestChi2 = m_tools->isITkGeometry() ? nMaxClusters-1 : nMaxClusters;
281  /// best cluster seen
282  const InDet::SiCluster* cl = nullptr;
283 
284  T* sibegin = std::any_cast<T>(&m_sibegin);
285  T* siend = std::any_cast<T>(&m_siend);
286  if (sibegin==nullptr or siend==nullptr) return 0;
287 
288  /// loop over all hits on this element
289  for (T p=*sibegin; p!=*siend; ++p) {
290  const InDet::SiCluster* c = static_cast<const InDet::SiCluster*>(*p);
291  const Amg::Vector2D& M = c->localPosition();
292 
293  /// estimate covariance by pitch
294  double MV0 = c->width().phiR() ;
295  /// include 1/sqrt(12) factor and add to contribution from pred
296  double v0 = s_oneOverTwelve*(MV0*MV0)+PV0;
297  /// get residual in X
298  double r0 = M[0]-P0;
299  /// one over cov
300  double d = 1./v0;
301  /// chi2
302  double x = (r0*r0)*d;
303 
304  /// stop if chi2 too large
305  if(x>Xc) continue;
306  /// estimate of Y residual
307  double dP1 = (P1-M[1])+PV1*d*r0;
308 
309  /// if residual beyond half-length: assume hit is on very end of strip and
310  /// re-evaluate chi2 with a penalty term on the Y coordinate
311  if(std::abs(dP1) > m_halflength) {
312  double r1 = M[1]-P1;
313  if(m_tools->isITkGeometry()){
314  if(dP1 > m_halflength) r1 += m_halflength;
315  else r1 -= m_halflength;
316  }
317  else{
318  r1 = (dP1 > m_halflength ?m_halflength-P1 : -(m_halflength+P1));
319  }
320  double v1 = PV1;
321  double v2 = PV2;
322  /// determinant of 2D covariance (only contribution from hit on X,X)
323  d = v0*v2-v1*v1 ;
324  if(d<=0.) continue;
325  /// updated chi2 estimate, with penalty for being at end of strip
326  x = (r0*(r0*v2-r1*v1)+r1*(r1*v0-r0*v1))/d;
327  /// and re-apply cut
328  if(x>Xc) continue;
329  } /// end of branch for Y residual beyond half-length
330 
331  /// if we satisfy the minimum criterion:
332  if(x < Xm) {
333  /// create a cluster link for our firend
334  InDet::SiClusterLink_xk l(c,x);
335  /// This misleadingly named method actually is a sorting
336  /// Mechanism: Our current link will take the place of the one
337  /// it replaces in ascending chi², and we will be
338  /// left with a link to the worst chi² one
339  for(int i=0; i!=nLinksFound; ++i) L[i].Comparison(l);
340  /// and place the worst link so far at the end if space.
341  if(nLinksFound<=nMaxClusters)
342  L[nLinksFound++]=l;
343  /// otherwise, update chi2 cut to not bother with worse clusters anymore
344  else Xm=L[nMaxClustersForBestChi2].xi2();
345  /// update looser chi2 to match the tighter one - we have at least one good link now
346  Xc = Xm;
347  }
348  /// if we didn't find any good links yet, keep this one as best so far
349  /// update looser chi2 cut to only keep better candidates in the future
350  else if(!nLinksFound) {Xc = x; cl = c;}
351  }
352  /// if we found no clusters satisfying the chi2 cut, we return
353  /// a single link to the best (least bad) candidate
354  if(cl && !nLinksFound) {L[nLinksFound++].Set(cl,Xc);}
355  return nLinksFound;
356 }
357 
358 ///////////////////////////////////////////////////////////////////
359 // Search closest cluster with stereo angle
360 ///////////////////////////////////////////////////////////////////
361 
362 template <typename T>
363 int InDet::SiTrajectoryElement_xk::searchClustersWithStereoAss
364 (Trk::PatternTrackParameters& Tp,InDet::SiClusterLink_xk* L, const Trk::PRDtoTrackMap &prd_to_track_map)
365 {
366 
367  /// for detailed documentation, please see searchClustersWithStereo
368 
369  if (m_detstatus<=0) return 0;
370 
371  const AmgSymMatrix(5) & Tp_cov = *Tp.covariance();
372  int nLinksFound = 0;
373  double P0 = Tp.parameters()[0];
374  double P1 = Tp.parameters()[1];
375  double PV0 = Tp_cov(0, 0);
376  double PV1 = Tp_cov(0, 1);
377  double PV2 = Tp_cov(1, 1);
378  double Xc = m_xi2maxlink;
379  double Xl = m_xi2maxlink;
380  double Xm = m_xi2max ;
381  int nMaxClusters = m_tools->isITkGeometry() ? m_tools->maxclusters() : 9;
382  int nMaxClustersForBestChi2 = m_tools->isITkGeometry() ? nMaxClusters-1 : nMaxClusters;
383 
384  const InDet::SiCluster* cl = nullptr;
385 
386  T* sibegin = std::any_cast<T>(&m_sibegin);
387  T* siend = std::any_cast<T>(&m_siend);
388  if (sibegin==nullptr or siend==nullptr) return 0;
389  for (T p=*sibegin; p!=*siend; ++p) {
390  const InDet::SiCluster* c = static_cast<const InDet::SiCluster*>(*p);
391 
392  /// check the PRD association map, skip clusters already in use
393  if (prd_to_track_map.isUsed(*c)) continue;
394 
395  const Amg::Vector2D& M = c->localPosition();
396  const Amg::MatrixX& V = c->localCovariance();
397 
398  double MV0 = V(0,0);
399  double MV1 = V(1,0);
400  double MV2 = V(1,1);
401  double v0 = MV0+PV0;
402  double v1 = MV1+PV1;
403  double v2 = MV2+PV2;
404  double r0 = M[0]-P0;
405  double r1 = M[1]-P1;
406  double d = v0*v2-v1*v1;
407  if(d<=0.) continue;
408  d=1./d;
409  double x = (r0*(r0*v2-r1*v1)+r1*(r1*v0-r0*v1))*d;
410  if(x > Xc) continue;
411 
412  r1 = fabs(r1+d*((PV1*v2-PV2*v1)*r0+(PV2*v0-PV1*v1)*r1));
413  x -= (r1*r1)/MV2 ;
414  r1 -= m_halflength ;
415 
416  if(r1 > 0. && (x+=((r1*r1)/PV2)) > Xc) continue;
417 
418  if(x < Xm) {
419  InDet::SiClusterLink_xk l(c,x);
420  /// This misleadingly named method actually is a sorting
421  /// Mechanism: Our current link will take the place of the one
422  /// it replaces in ascending chi², and we will be
423  /// left with a link to the worst chi² one
424  for(int i=0; i!=nLinksFound; ++i) L[i].Comparison(l);
425  if(nLinksFound<=nMaxClusters)
426  L[nLinksFound++]=l;
427  else Xm=L[nMaxClustersForBestChi2].xi2();
428  Xc = Xm+6.;
429  }
430  else if(!nLinksFound && x < Xl) {Xl = x; Xc = x+6.; cl = c;}
431  }
432  if(cl && !nLinksFound) {L[nLinksFound++].Set(cl,Xl);}
433  return nLinksFound;
434 }
435 
436 ///////////////////////////////////////////////////////////////////
437 // Search closest cluster without stereo angle for pixels
438 ///////////////////////////////////////////////////////////////////
439 
440 template <typename T>
441 int InDet::SiTrajectoryElement_xk::searchClustersWithoutStereoAssPIX
442 (Trk::PatternTrackParameters& Tp,InDet::SiClusterLink_xk* L, const Trk::PRDtoTrackMap &prd_to_track_map)
443 {
444  /// for detailed documentation, please see searchClustersWithoutStereoPIX
445  if (m_detstatus<=0) return 0;
446 
447  const AmgSymMatrix(5) & Tp_cov = *Tp.covariance();
448  int nLinksFound = 0;
449  double P0 = Tp.parameters()[0];
450  double P1 = Tp.parameters()[1];
451  double PV0 = Tp_cov(0, 0);
452  double PV1 = Tp_cov(0, 1);
453  double PV2 = Tp_cov(1, 1);
454  double Xc = m_xi2maxlink;
455  double Xm = m_xi2max ;
456  int nMaxClusters = m_tools->isITkGeometry() ? m_tools->maxclusters() : 9;
457  int nMaxClustersForBestChi2 = m_tools->isITkGeometry() ? nMaxClusters-1 : nMaxClusters;
458 
459  const InDet::SiCluster* cl = nullptr;
460 
461  T* sibegin = std::any_cast<T>(&m_sibegin);
462  T* siend = std::any_cast<T>(&m_siend);
463  if (sibegin==nullptr or siend==nullptr) return 0;
464 
465  for (T p=*sibegin; p!=*siend; ++p) {
466  const InDet::SiCluster* c = static_cast<const InDet::SiCluster*>(*p);
467  if (prd_to_track_map.isUsed(*c)) continue;
468  const Amg::Vector2D& M = c->localPosition();
469 
470  double MV0 = c->width().phiR();
471  double MV2 = c->width().z ();
472 
473  double r0 = M[0]-P0, r02 = r0*r0;
474  double r1 = M[1]-P1, r12 = r1*r1;
475 
476  double v0 = s_oneOverTwelve*(MV0*MV0)+PV0; if(r02 >(Xc*v0)) continue;
477  double v2 = s_oneOverTwelve*(MV2*MV2)+PV2; if(r12 >(Xc*v2)) continue;
478  double v1 = PV1;
479  double d = v0*v2-v1*v1; if( d<=0. ) continue;
480  double x = (r02*v2+r12*v0-(r0*r1)*(2.*v1))/d;
481 
482  if(x>Xc) continue;
483 
484  if(x < Xm) {
485  InDet::SiClusterLink_xk l(c,x);
486  /// This misleadingly named method actually is a sorting
487  /// Mechanism: Our current link will take the place of the one
488  /// it replaces in ascending chi², and we will be
489  /// left with a link to the worst chi² one
490  for(int i=0; i!=nLinksFound; ++i) L[i].Comparison(l);
491  if(nLinksFound<=nMaxClusters)
492  L[nLinksFound++]=l;
493  else Xm=L[nMaxClustersForBestChi2].xi2();
494  Xc = Xm;
495  }
496  else if(!nLinksFound) {Xc = x; cl = c;}
497  }
498  if(cl && !nLinksFound) {L[nLinksFound++].Set(cl,Xc);}
499  return nLinksFound;
500 }
501 
502 ///////////////////////////////////////////////////////////////////
503 // Search closest cluster without stereo angle for SCT
504 ///////////////////////////////////////////////////////////////////
505 
506 template <typename T>
507 int InDet::SiTrajectoryElement_xk::searchClustersWithoutStereoAssSCT
508 (Trk::PatternTrackParameters& Tp,InDet::SiClusterLink_xk* L, const Trk::PRDtoTrackMap &prd_to_track_map)
509 {
510  /// for detailed documentation, please see searchClustersWithoutStereoSCT
511  if (m_detstatus<=0) return 0;
512 
513  const AmgSymMatrix(5) & Tp_cov = *Tp.covariance();
514  int nLinksFound = 0;
515  double P0 = Tp.parameters()[0];
516  double P1 = Tp.parameters()[1];
517  double PV0 = Tp_cov(0, 0);
518  double PV1 = Tp_cov(0, 1);
519  double PV2 = Tp_cov(1, 1);
520  double Xc = m_xi2maxlink;
521  double Xm = m_xi2max ;
522  int nMaxClusters = m_tools->isITkGeometry() ? m_tools->maxclusters() : 9;
523  int nMaxClustersForBestChi2 = m_tools->isITkGeometry() ? nMaxClusters-1 : nMaxClusters;
524 
525  const InDet::SiCluster* cl = nullptr;
526 
527  T* sibegin = std::any_cast<T>(&m_sibegin);
528  T* siend = std::any_cast<T>(&m_siend);
529  if (sibegin==nullptr or siend==nullptr) return 0;
530 
531  for (T p=*sibegin; p!=*siend; ++p) {
532  const InDet::SiCluster* c = static_cast<const InDet::SiCluster*>(*p);
533  if (prd_to_track_map.isUsed(*c)) continue;
534  const Amg::Vector2D& M = c->localPosition();
535 
536  double MV0 = c->width().phiR() ;
537  double v0 = s_oneOverTwelve*(MV0*MV0)+PV0;
538  double r0 = M[0]-P0;
539  double d = 1./v0;
540  double x = (r0*r0)*d;
541 
542  if(x>Xc) continue;
543 
544  double dP1 = (P1-M[1])+PV1*d*r0;
545 
546  if(fabs(dP1) > m_halflength) {
547 
548  double r1 = M[1]-P1;
549  if(m_tools->isITkGeometry()){
550  if(dP1 > m_halflength) r1 += m_halflength;
551  else r1 -= m_halflength;
552  }
553  else{
554  r1 = (dP1 > m_halflength ?m_halflength-P1 : -(m_halflength+P1));
555  }
556 
557  double v1 = PV1;
558  double v2 = PV2;
559  d = v0*v2-v1*v1 ; if(d<=0.) continue;
560  x = (r0*(r0*v2-r1*v1)+r1*(r1*v0-r0*v1))/d;
561  if(x>Xc) continue;
562  }
563 
564  if(x < Xm) {
565  InDet::SiClusterLink_xk l(c,x);
566  /// This misleadingly named method actually is a sorting
567  /// Mechanism: Our current link will take the place of the one
568  /// it replaces in ascending chi², and we will be
569  /// left with a link to the worst chi² one
570  for(int i=0; i!=nLinksFound; ++i) L[i].Comparison(l);
571  if(nLinksFound<=nMaxClusters)
572  L[nLinksFound++]=l;
573  else Xm=L[nMaxClustersForBestChi2].xi2();
574  Xc = Xm;
575  }
576  else if(!nLinksFound) {Xc = x; cl = c;}
577  }
578  if(cl && !nLinksFound) {L[nLinksFound++].Set(cl,Xc);}
579  return nLinksFound;
580 }
581 
582 ///////////////////////////////////////////////////////////////////
583 // Set trajectory element
584 ///////////////////////////////////////////////////////////////////
585 
586 // T = InDet::SiClusterCollection::const_iterator or
587 // InDet::PixelClusterCollection::const_iterator or
588 // InDet::SCT_ClusterCollection::const_iterator
589 template <typename T>
590 bool InDet::SiTrajectoryElement_xk::set
591 (int st,
592  const InDet::SiDetElementBoundaryLink_xk*& dl,
593  const T& sb,
594  const T& se,
595  const InDet::SiCluster* si,
596  [[maybe_unused]] const EventContext& ctx)
597 {
598  if (std::is_same<T, InDet::SiClusterCollection::const_iterator>::value) m_itType = SiClusterColl;
599  else if (std::is_same<T, InDet::PixelClusterCollection::const_iterator>::value) m_itType = PixelClusterColl;
600  else if (std::is_same<T, InDet::SCT_ClusterCollection::const_iterator>::value) m_itType = SCT_ClusterColl;
601 
602  m_fieldMode = false;
603  if(m_tools->fieldTool().magneticFieldMode()!=0) m_fieldMode = true;
604  m_status = 0 ;
605  m_detstatus = st ;
606  m_nMissing = 0 ;
607  m_nlinksForward = 0 ;
608  m_nlinksBackward = 0 ;
609  m_nholesForward = 0 ;
610  m_nholesBackward = 0 ;
611  m_dholesForward = 0 ;
612  m_dholesBackward = 0 ;
613  m_nclustersForward = 0 ;
614  m_nclustersBackward = 0 ;
615  m_npixelsBackward = 0 ;
616  m_ndfForward = 0 ;
617  m_ndfBackward = 0 ;
618  m_ntsos = 0 ;
619  m_detelement = dl->detElement() ;
620  m_detlink = dl ;
621  m_surface = &m_detelement->surface(); if(!m_surface && m_tools->isITkGeometry()) return false;
622  m_sibegin = sb ;
623  m_siend = se ;
624  m_cluster = si ;
625  m_clusterOld = si ;
626  m_clusterNoAdd = 0 ;
627  m_stereo = false ;
628  m_xi2Forward = 10000. ;
629  m_xi2Backward = 10000. ;
630  m_xi2totalForward = 0. ;
631  m_xi2totalBackward = 0. ;
632  m_tools->electron() ? m_xi2max = m_tools->xi2maxBrem() : m_xi2max = m_tools->xi2max();
633  m_halflength = 0. ;
634  m_detelement->isSCT() ? m_ndf=1 : m_ndf=2;
635 
636  if(m_tools->heavyion()) {
637  if(m_ndf==2) {m_xi2max = 13.; m_xi2maxNoAdd = 13.;}
638  else {m_xi2max = 4.; m_xi2maxNoAdd = 8.;}
639  }
640 
641  noiseInitiate() ;
642  (m_detelement->isSCT() && (m_detelement->design().shape()==InDetDD::Trapezoid || m_detelement->design().shape()==InDetDD::Annulus) ) ?
643  m_stereo = true : m_stereo = false;
644 
645  if(m_detstatus && m_ndf == 1){
646  m_halflength = si ? si->width().z()*.5 : (*sb)->width().z()*.5;
647  }
648 
649  if(!m_detstatus) {
650 
651  IdentifierHash idHash = m_detelement->identifyHash();
652  if (m_ndf==2) {
653  if (m_tools->pixelStatus() ? !m_tools->pixelStatus()->at(idHash) : !m_tools->pixcond()->isGood(idHash, ctx) ) {
654  m_detstatus = -1;
655  }
656  }
657  else {
658  if (m_tools->sctStatus() ? !m_tools->sctStatus()->at(idHash) : !m_tools->sctcond()->isGood(idHash, ctx) ) {
659  m_detstatus = -1;
660  }
661  }
662  VALIDATE_STATUS_ARRAY(m_tools->pixelStatus() && m_tools->sctStatus(), ((m_ndf==2) ? m_tools->pixelStatus() : m_tools->sctStatus())->at(idHash), ( (m_ndf==2) ? m_tools->pixcond() : m_tools->sctcond())->isGood(idHash, ctx));
663 
664  }
665 
666  const Amg::Transform3D& tr = m_surface->transform();
667  //Initialise material corrections
668  m_radlength = 0.03;
669  if(m_tools->isITkGeometry()) m_radlength = 0.04;
670  /// add material for the forward pixels
671  if(m_ndf == 2 && std::abs(tr(2,2)) >= 1.) m_radlength = .07;
672  if(m_tools->isITkGeometry()) m_radlengthN=m_radlength;
673 
674  m_localTransform[ 0] = tr(0,0); m_localTransform[ 1]=tr(1,0); m_localTransform[ 2]=tr(2,0);
675  m_localTransform[ 3] = tr(0,1); m_localTransform[ 4]=tr(1,1); m_localTransform[ 5]=tr(2,1);
676  m_localTransform[ 6] = tr(0,2); m_localTransform[ 7]=tr(1,2); m_localTransform[ 8]=tr(2,2);
677  m_localTransform[ 9] = tr(0,3); m_localTransform[10]=tr(1,3); m_localTransform[11]=tr(2,3);
678  m_localTransform[12] = m_localTransform[ 9]*m_localTransform[ 6]+m_localTransform[10]*m_localTransform[ 7]+m_localTransform[11]*m_localTransform[ 8];
679  m_localDir[0] = 1.;
680  m_localDir[1] = 0.;
681  m_localDir[2] = 0.;
682  return true;
683 }
684 
685 ///////////////////////////////////////////////////////////////////
686 // Forward propagation for search closest cluster
687 ///////////////////////////////////////////////////////////////////
688 
689 // T = InDet::SiClusterCollection::const_iterator or
690 // InDet::PixelClusterCollection::const_iterator or
691 // InDet::SCT_ClusterCollection::const_iterator
692 template <typename T>
693 bool InDet::SiTrajectoryElement_xk::ForwardPropagationForClusterSeach
694 (int n,
695  const Trk::TrackParameters& Tpa,
696  const InDet::SiDetElementBoundaryLink_xk*& dl,
697  const T& sb,
698  const T& se,
699  const EventContext& ctx)
700 {
701  // remove case if you have trajectory element without actual detector element
702  // this happens if you have added a dead cylinder
703  if(!m_detelement) {
704  return false;
705  }
706 
707  if (std::is_same<T, InDet::SiClusterCollection::const_iterator>::value) m_itType = SiClusterColl;
708  else if (std::is_same<T, InDet::PixelClusterCollection::const_iterator>::value) m_itType = PixelClusterColl;
709  else if (std::is_same<T, InDet::SCT_ClusterCollection::const_iterator>::value) m_itType = SCT_ClusterColl;
710 
711  m_detstatus = 1 ;
712  m_sibegin = sb ;
713  m_siend = se ;
714  m_detelement = dl->detElement() ;
715  m_detlink = dl ;
716  m_surface = &m_detelement->surface();
717  m_detelement->isSCT() ? m_ndf=1 : m_ndf=2;
718  m_halflength = 0. ;
719  m_stereo = false ;
720 
721  (m_detelement->isSCT() && m_detelement->design().shape()==InDetDD::Trapezoid) ?
722  m_stereo = true : m_stereo = false;
723 
724  if(m_detstatus && m_ndf == 1) m_halflength = (*sb)->width().z()*.5;
725 
726  if(!n) {
727  Trk::PatternTrackParameters Tp; if(!Tp.production(&Tpa)) return false;
728  if(!propagateParameters(Tp,m_parametersPredForward,m_step)) return false;
729 
730  if(!m_parametersPredForward.iscovariance()) {
731 
732  double cv[15]={ .02 ,
733  0., .02,
734  0., 0.,.000001,
735  0., 0., 0.,.000001,
736  0., 0., 0., 0.,.000000001};
737 
738  m_parametersPredForward.setCovariance(cv);
739  }
740  }
741  else {
742  if(!propagate(m_parametersPredForward,m_parametersPredForward,m_step,ctx)) return false;
743  }
744  m_nlinksForward=searchClusters(m_parametersPredForward,m_linkForward);
745  return true;
746 }
747 
748 ///////////////////////////////////////////////////////////////////
749 // Forward propagation for search closest cluster
750 ///////////////////////////////////////////////////////////////////
751 
752 // T = InDet::SiClusterCollection::const_iterator or
753 // InDet::PixelClusterCollection::const_iterator or
754 // InDet::SCT_ClusterCollection::const_iterator
755 template <typename T>
756 void InDet::SiTrajectoryElement_xk::CloseClusterSeach
757 (Trk::PatternTrackParameters& Tpa,
758  const InDet::SiDetElementBoundaryLink_xk*& dl,
759  const T& sb,
760  const T& se)
761 {
762  if (std::is_same<T, InDet::SiClusterCollection::const_iterator>::value) m_itType = SiClusterColl;
763  else if (std::is_same<T, InDet::PixelClusterCollection::const_iterator>::value) m_itType = PixelClusterColl;
764  else if (std::is_same<T, InDet::SCT_ClusterCollection::const_iterator>::value) m_itType = SCT_ClusterColl;
765 
766  m_detstatus = 1 ;
767  m_cluster = 0 ;
768  m_sibegin = sb ;
769  m_siend = se ;
770  m_detelement = dl->detElement() ;
771  m_detlink = dl ;
772  m_surface = &m_detelement->surface();
773  m_detelement->isSCT() ? m_ndf=1 : m_ndf=2;
774  m_halflength = 0. ;
775  m_stereo = false ;
776 
777  (m_detelement->isSCT() && m_detelement->design().shape()==InDetDD::Trapezoid) ?
778  m_stereo = true : m_stereo = false;
779  if(m_detstatus && m_ndf == 1) m_halflength = (*sb)->width().z()*.5;
780  /// in this case, we do not wish to overwrite the m_inside member.
781  /// So store the old value and restore after the check.
782  int inside_originalVal = m_inside;
783  checkBoundaries(Tpa);
784  if(m_inside > 0 || !searchClusters(Tpa,m_linkForward)) {
785  m_inside = inside_originalVal;
786  return;
787  }
788  m_inside = inside_originalVal;
789  m_cluster = m_linkForward[0].cluster();
790  m_xi2Forward = m_linkForward[0].xi2 ();
791 }
792 
793 
794 inline const Trk::PatternTrackParameters *
795 InDet::SiTrajectoryElement_xk::parameters() const
796 {
797  // logic from SiTrajectoryElement_xk::trackParameters(bool, int Q) with Q==1
798 
799  if (m_status == 1) {
800  if (m_cluster) {
801  return &m_parametersUpdatedForward;
802  } else {
803  return &m_parametersPredForward;
804  }
805  } else if (m_status == 2) {
806  if (m_cluster) {
807  return &m_parametersUpdatedBackward;
808  } else {
809  return &m_parametersPredBackward;
810  }
811  } else if (m_status == 3) {
812  if (m_cluster) {
813  return &m_parametersUpdatedBackward;
814  }
815  }
816  return nullptr;
817 }
818 
Generated on Sun Dec 22 2024 21:18:20 for ATLAS Offline Software by   doxygen 1.8.18