ATLAS Offline Software
TrigInDetTrackFitter.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 
7 // TrigInDetTrackFitter tool
8 // -------------------------------
9 // ATLAS Collaboration
10 //
11 // 01.09.2005 Package created
12 //
13 // Author: Dmitry Emeliyanov, RAL
14 // e-mail: D.Emeliyanov@rl.ac.uk
15 //
17 
18 
19 #include <cmath>
20 #include <iostream>
21 #include <memory>
22 #include "GaudiKernel/SystemOfUnits.h"
23 
24 #include "TrigInDetEvent/TrigInDetTrack.h"
25 #include "TrigInDetEvent/TrigInDetTrackCollection.h"
26 
27 #include "TrkTrack/Track.h"
28 #include "TrkTrack/TrackInfo.h"
29 #include "TrkTrack/TrackStateOnSurface.h"
30 #include "TrkParameters/TrackParameters.h"
31 
32 #include "TrkDistributedKalmanFilter/TrkBaseNode.h"
33 #include "TrkDistributedKalmanFilter/TrkFilteringNodes.h"
34 #include "TrkDistributedKalmanFilter/TrkTrackState.h"
35 #include "TrkDistributedKalmanFilter/TrkPlanarSurface.h"
36 
37 #include "TrigInDetToolInterfaces/ITrigInDetTrackFitter.h"
38 #include "TrigInDetTrackFitter.h"
39 #include "TrkToolInterfaces/IRIO_OnTrackCreator.h"
40 #include "AthenaBaseComps/AthMsgStreamMacros.h"
41 #include "AthenaBaseComps/AthCheckMacros.h"
42 
43 #include "AtlasDetDescr/AtlasDetectorID.h"
44 #include "InDetIdentifier/SCT_ID.h"
45 #include "InDetIdentifier/PixelID.h"
46 #include "TrigInDetToolInterfaces/TrigL2HitResidual.h"
47 
48 #include "TrkRIO_OnTrack/RIO_OnTrack.h"
49 #include "TrkPrepRawData/PrepRawData.h"
50 #include "TrkSurfaces/Surface.h"
51 
52 #include "TrkEventPrimitives/ParticleHypothesis.h"
53 
54 #include "TrkEventPrimitives/FitQuality.h"
55 #include "TrkEventPrimitives/FitQualityOnSurface.h"
56 #include "MagFieldConditions/AtlasFieldCacheCondObj.h"
57 
58 TrigInDetTrackFitter::TrigInDetTrackFitter(const std::string& t,
59  const std::string& n,
60  const IInterface* p ): AthAlgTool(t,n,p),
61  m_trackMaker("TrigDkfTrackMakerTool")
62 {
63  declareInterface< ITrigInDetTrackFitter >( this );
64 
65  declareProperty( "doMultScattering", m_doMultScatt = true);
66  declareProperty( "doBremmCorrection", m_doBremm=false);
67  declareProperty( "Chi2Cut", m_DChi2 = 1000.0);
68  declareProperty( "correctClusterPos", m_correctClusterPos = false);
69  declareProperty( "ROTcreator", m_ROTcreator, "ROTcreatorTool" );
70  m_nTracksTotal = 0;
71  m_fitErrorsUnresolved=0;
72  m_fitErrorsDivergence=0;
73  m_fitErrorsLowPt=0;
74 }
75 
76 StatusCode TrigInDetTrackFitter::initialize()
77 {
78  ATH_CHECK(m_trackMaker.retrieve());
79  if (m_correctClusterPos) {
80  ATH_CHECK(m_ROTcreator.retrieve());
81  }
82  ATH_CHECK( m_fieldCondObjInputKey.initialize());
83  ATH_CHECK(detStore()->retrieve(m_idHelper, "AtlasID"));
84  ATH_CHECK(detStore()->retrieve(m_pixelId, "PixelID"));
85  ATH_CHECK(detStore()->retrieve(m_sctId, "SCT_ID"));
86 
87  return StatusCode::SUCCESS;
88 
89 }
90 
91 StatusCode TrigInDetTrackFitter::finalize()
92 {
93  ATH_MSG_INFO("==============================================================");
94  ATH_MSG_INFO("TrigInDetTrackFitter::finalize() - LVL2 Track fit Statistics: ");
95  ATH_MSG_INFO(" N tracks = "<<m_nTracksTotal);
96  ATH_MSG_INFO("Problems detected: ");
97  ATH_MSG_INFO("Unresolved spacepoints :"<< m_fitErrorsUnresolved);
98  ATH_MSG_INFO("Extrapolator divergence:"<< m_fitErrorsDivergence);
99  ATH_MSG_INFO("pT falls below 200 MeV :"<< m_fitErrorsLowPt);
100  ATH_MSG_INFO("==============================================================");
101  return StatusCode::SUCCESS;
102 }
103 
104 void TrigInDetTrackFitter::getMagneticField(double r[3],double* B, MagField::AtlasFieldCache& fieldCache) const {
105  B[0]=0.0;B[1]=0.0;B[2]=0.0;
106  double field[3];
107  fieldCache.getField(r,field);//field is returned in kT
108  for(int i=0;i<3;i++) B[i]=field[i]/Gaudi::Units::kilogauss;//convert to kG
109 }
110 
111 void TrigInDetTrackFitter::correctScale(Trk::TrkTrackState* pTS) const {
112 
113  double Rf[5];
114  double Gf[5][5];
115  int i,j;
116 
117  for(i=0;i<4;i++) Rf[i] = pTS->getTrackState(i);
118  Rf[4] = 0.001*pTS->getTrackState(4);
119 
120  for(i=0;i<4;i++)
121  for(j=0;j<4;j++) Gf[i][j] = pTS->getTrackCovariance(i,j);
122 
123  Gf[0][4] = Gf[4][0] = pTS->getTrackCovariance(0,4)/1000.0;
124  Gf[1][4] = Gf[4][1] = pTS->getTrackCovariance(1,4)/1000.0;
125  Gf[2][4] = Gf[4][2] = pTS->getTrackCovariance(2,4)/1000.0;
126  Gf[3][4] = Gf[4][3] = pTS->getTrackCovariance(3,4)/1000.0;
127  Gf[4][4] = pTS->getTrackCovariance(4,4)/1000000.0;
128 
129  pTS->setTrackState(Rf);
130  pTS->setTrackCovariance(Gf);
131 }
132 
133 Trk::TrkTrackState* TrigInDetTrackFitter::extrapolate(Trk::TrkTrackState* pTS,
134  Trk::TrkPlanarSurface* pSB,
135  Trk::TrkPlanarSurface* pSE,
136  MagField::AtlasFieldCache& fieldCache) const
137 {
138  const double C=0.02999975/1000.0;//using GeV internally
139  const double minStep=30.0;
140 
141  double J[5][5],Rf[5],AG[5][5],Gf[5][5],A[5][5];
142  int i,j,m;
143 
144  bool samePlane=false;
145 
146  if(pSB!=nullptr)
147  {
148  double diff=0.0;
149  for(i=0;i<4;i++) diff+=fabs(pSE->getPar(i)-pSB->getPar(i));
150  if(diff<1e-5) {
151  samePlane=true;
152  m_fitErrorsUnresolved++;
153  }
154  }
155 
156  if(!samePlane) {
157 
158  double gP[3],gPi[3],lP[3],gV[3],a,b,c,s,J0[7][5],descr,CQ,Ac,Av,Cc;
159  double V[3],P[3],M[3][3],D[4],Jm[7][7],
160  J1[5][7],gB[3],gBi[3],gBf[3],dBds[3],Buf[5][7],DVx,DVy,DVz;
161  int nStep,nStepMax;
162  double sl,ds;
163 
164  double sint,cost,sinf,cosf;
165  sint=sin(pTS->getTrackState(3));cosf=cos(pTS->getTrackState(2));
166  sinf=sin(pTS->getTrackState(2));cost=cos(pTS->getTrackState(3));
167  gV[0]=sint*cosf;gV[1]=sint*sinf;gV[2]=cost;CQ=C*pTS->getTrackState(4);
168 
169  memset(&J0[0][0],0,sizeof(J0));
170 
171  if(pSB!=nullptr)
172  {
173  double L[3][3];
174  lP[0]=pTS->getTrackState(0);lP[1]=pTS->getTrackState(1);lP[2]=0.0;
175  pSB->transformPointToGlobal(lP,gP);
176  for(i=0;i<3;i++) for(j=0;j<3;j++) L[i][j]=pSB->getInvRotMatrix(i,j);
177 
178  J0[0][0]=L[0][0];J0[0][1]=L[0][1];
179  J0[1][0]=L[1][0];J0[1][1]=L[1][1];
180  J0[2][0]=L[2][0];J0[2][1]=L[2][1];
181  J0[3][2]=-sinf*sint;J0[3][3]=cosf*cost;
182  J0[4][2]= cosf*sint;J0[4][3]=sinf*cost;
183  J0[5][3]=-sint;
184  J0[6][4]=1.0;
185  }
186  else
187  {
188  gP[0]=-pTS->getTrackState(0)*sinf;
189  gP[1]= pTS->getTrackState(0)*cosf;
190  gP[2]= pTS->getTrackState(1);
191  J0[0][0]=-sinf;J0[0][2]=-pTS->getTrackState(0)*cosf;
192  J0[1][0]= cosf;J0[1][2]=-pTS->getTrackState(0)*sinf;
193  J0[2][1]=1.0;
194  J0[3][2]=-sinf*sint;J0[3][3]=cosf*cost;
195  J0[4][2]= cosf*sint;J0[4][3]=sinf*cost;
196  J0[5][3]=-sint;
197  J0[6][4]=1.0;
198  }
199  for(i=0;i<4;i++) D[i]=pSE->getPar(i);
200  for(i=0;i<3;i++) gPi[i]=gP[i];
201 
202  getMagneticField(gP,gB, fieldCache);
203 
204  for(i=0;i<3;i++) gBi[i]=gB[i];
205 
206  c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
207  b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
208  a=0.5*CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
209  gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
210  gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
211 
212  descr=b*b-4.0*a*c;
213 
214  if(descr<0.0)
215  {
216  // printf("D<0 - extrapolation failed\n");
217  return nullptr;
218  }
219 
220  bool useExpansion=true;
221  double ratio = 4*a*c/(b*b);
222 
223  if(fabs(ratio)>0.1)
224  useExpansion = false;
225 
226  if(useExpansion) {
227  sl=-c/b;
228  sl=sl*(1-a*sl/b);
229  }
230  else {
231  int signb = (b<0.0)?-1:1;
232  sl = (-b+signb*sqrt(descr))/(2*a);
233  }
234 
235  if(fabs(sl)<minStep) nStepMax=1;
236  else
237  {
238  nStepMax=(int)(fabs(sl)/minStep)+1;
239  }
240  if((nStepMax<0)||(nStepMax>1000))
241  {
242  return nullptr;
243  }
244  Av=sl*CQ;
245  Ac=0.5*sl*Av;
246  DVx=gV[1]*gB[2]-gV[2]*gB[1];
247  DVy=gV[2]*gB[0]-gV[0]*gB[2];
248  DVz=gV[0]*gB[1]-gV[1]*gB[0];
249 
250  P[0]=gP[0]+gV[0]*sl+Ac*DVx;
251  P[1]=gP[1]+gV[1]*sl+Ac*DVy;
252  P[2]=gP[2]+gV[2]*sl+Ac*DVz;
253  V[0]=gV[0]+Av*DVx;
254  V[1]=gV[1]+Av*DVy;
255  V[2]=gV[2]+Av*DVz;
256 
257  getMagneticField(P,gB,fieldCache);
258 
259  for(i=0;i<3;i++) gBf[i]=gB[i];
260  for(i=0;i<3;i++)
261  {
262  dBds[i]=(gBf[i]-gBi[i])/sl;
263  gB[i]=gBi[i];
264  }
265  nStep=nStepMax;
266  while(nStep>0)
267  {
268  c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
269  b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
270  a=0.5*CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
271  gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
272  gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
273 
274  ratio = 4*a*c/(b*b);
275  if(fabs(ratio)>0.1)
276  useExpansion = false;
277  else useExpansion = true;
278 
279  if(useExpansion) {
280  sl=-c/b;
281  sl=sl*(1-a*sl/b);
282  }
283  else {
284  descr=b*b-4.0*a*c;
285  if(descr<0.0)
286  {
287  // printf("D<0 - extrapolation failed\n");
288  return nullptr;
289  }
290  int signb = (b<0.0)?-1:1;
291  sl = (-b+signb*sqrt(descr))/(2*a);
292  }
293 
294  ds=sl/nStep;
295  Av=ds*CQ;
296  Ac=0.5*ds*Av;
297  DVx=gV[1]*gB[2]-gV[2]*gB[1];
298  DVy=gV[2]*gB[0]-gV[0]*gB[2];
299  DVz=gV[0]*gB[1]-gV[1]*gB[0];
300 
301  P[0]=gP[0]+gV[0]*ds+Ac*DVx;
302  P[1]=gP[1]+gV[1]*ds+Ac*DVy;
303  P[2]=gP[2]+gV[2]*ds+Ac*DVz;
304  V[0]=gV[0]+Av*DVx;
305  V[1]=gV[1]+Av*DVy;
306  V[2]=gV[2]+Av*DVz;
307  for(i=0;i<3;i++)
308  {
309  gV[i]=V[i];gP[i]=P[i];
310  }
311  for(i=0;i<3;i++) gB[i]+=dBds[i]*ds;
312  nStep--;
313  }
314  pSE->transformPointToLocal(gP,lP);
315  Rf[0]=lP[0];Rf[1]=lP[1];
316  Rf[2]=atan2(V[1],V[0]);
317 
318  if(fabs(V[2])>1.0)
319  {
320  return nullptr;
321  }
322 
323  Rf[3]=acos(V[2]);
324  Rf[4]=pTS->getTrackState(4);
325 
326  gV[0]=sint*cosf;gV[1]=sint*sinf;gV[2]=cost;
327 
328  for(i=0;i<4;i++) D[i]=pSE->getPar(i);
329  for(i=0;i<3;i++) gP[i]=gPi[i];
330 
331  for(i=0;i<3;i++)
332  {
333  gB[i]=0.5*(gBi[i]+gBf[i]);
334  }
335 
336  c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
337  b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
338  a=CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
339  gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
340  gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
341 
342  ratio = 4*a*c/(b*b);
343  if(fabs(ratio)>0.1)
344  useExpansion = false;
345  else useExpansion = true;
346 
347  if(useExpansion) {
348  s=-c/b;
349  s=s*(1-a*s/b);
350  }
351  else {
352  descr=b*b-4.0*a*c;
353  if(descr<0.0)
354  {
355  // printf("D<0 - extrapolation failed\n");
356  return nullptr;
357  }
358  int signb = (b<0.0)?-1:1;
359  s = (-b+signb*sqrt(descr))/(2*a);
360  }
361 
362  Av=s*CQ;
363  Ac=0.5*s*Av;
364  Cc=0.5*s*s*C;
365 
366  DVx=gV[1]*gB[2]-gV[2]*gB[1];
367  DVy=gV[2]*gB[0]-gV[0]*gB[2];
368  DVz=gV[0]*gB[1]-gV[1]*gB[0];
369 
370  P[0]=gP[0]+gV[0]*s+Ac*DVx;
371  P[1]=gP[1]+gV[1]*s+Ac*DVy;
372  P[2]=gP[2]+gV[2]*s+Ac*DVz;
373 
374  V[0]=gV[0]+Av*DVx;V[1]=gV[1]+Av*DVy;V[2]=gV[2]+Av*DVz;
375  if (std::abs(V[2]) > 1.0) {
376  return nullptr;
377  }
378 
379  pSE->transformPointToLocal(P,lP);
380 
381  memset(&Jm[0][0],0,sizeof(Jm));
382 
383  for(i=0;i<3;i++) for(j=0;j<3;j++) M[i][j]=pSE->getRotMatrix(i,j);
384 
385  double coeff[3], dadVx,dadVy,dadVz,dadQ,dsdx,dsdy,dsdz,dsdVx,dsdVy,dsdVz,dsdQ;
386  coeff[0]=-c*c/(b*b*b);
387  coeff[1]=c*(1.0+3.0*c*a/(b*b))/(b*b);
388  coeff[2]=-(1.0+2.0*c*a/(b*b))/b;
389 
390  dadVx=0.5*CQ*(-D[1]*gB[2]+D[2]*gB[1]);
391  dadVy=0.5*CQ*( D[0]*gB[2]-D[2]*gB[0]);
392  dadVz=0.5*CQ*(-D[0]*gB[1]+D[1]*gB[0]);
393  dadQ=0.5*C*(D[0]*DVx+D[1]*DVy+D[2]*DVz);
394 
395  dsdx=coeff[2]*D[0];
396  dsdy=coeff[2]*D[1];
397  dsdz=coeff[2]*D[2];
398  dsdVx=coeff[0]*dadVx+coeff[1]*D[0];
399  dsdVy=coeff[0]*dadVy+coeff[1]*D[1];
400  dsdVz=coeff[0]*dadVz+coeff[1]*D[2];
401  dsdQ=coeff[0]*dadQ;
402 
403  Jm[0][0]=1.0+V[0]*dsdx;
404  Jm[0][1]= V[0]*dsdy;
405  Jm[0][2]= V[0]*dsdz;
406 
407  Jm[0][3]= s+V[0]*dsdVx;
408  Jm[0][4]= V[0]*dsdVy+Ac*gB[2];
409  Jm[0][5]= V[0]*dsdVz-Ac*gB[1];
410  Jm[0][6]= V[0]*dsdQ+Cc*DVx;
411 
412  Jm[1][0]= V[1]*dsdx;
413  Jm[1][1]=1.0+V[1]*dsdy;
414  Jm[1][2]= V[1]*dsdz;
415 
416  Jm[1][3]= V[1]*dsdVx-Ac*gB[2];
417  Jm[1][4]= s+V[1]*dsdVy;
418  Jm[1][5]= V[1]*dsdVz+Ac*gB[0];
419  Jm[1][6]= V[1]*dsdQ+Cc*DVy;
420 
421  Jm[2][0]= V[2]*dsdx;
422  Jm[2][1]= V[2]*dsdy;
423  Jm[2][2]=1.0+V[2]*dsdz;
424  Jm[2][3]= V[2]*dsdVx+Ac*gB[1];
425  Jm[2][4]= V[2]*dsdVy-Ac*gB[0];
426  Jm[2][5]= s+V[2]*dsdVz;
427  Jm[2][6]= V[2]*dsdQ+Cc*DVz;
428 
429  Jm[3][0]=dsdx*CQ*DVx;
430  Jm[3][1]=dsdy*CQ*DVx;
431  Jm[3][2]=dsdz*CQ*DVx;
432 
433  Jm[3][3]=1.0+dsdVx*CQ*DVx;
434  Jm[3][4]=CQ*(dsdVy*DVx+s*gB[2]);
435  Jm[3][5]=CQ*(dsdVz*DVx-s*gB[1]);
436 
437  Jm[3][6]=(CQ*dsdQ+C*s)*DVx;
438 
439  Jm[4][0]=dsdx*CQ*DVy;
440  Jm[4][1]=dsdy*CQ*DVy;
441  Jm[4][2]=dsdz*CQ*DVy;
442 
443  Jm[4][3]=CQ*(dsdVx*DVy-s*gB[2]);
444  Jm[4][4]=1.0+dsdVy*CQ*DVy;
445  Jm[4][5]=CQ*(dsdVz*DVy+s*gB[0]);
446 
447  Jm[4][6]=(CQ*dsdQ+C*s)*DVy;
448 
449  Jm[5][0]=dsdx*CQ*DVz;
450  Jm[5][1]=dsdy*CQ*DVz;
451  Jm[5][2]=dsdz*CQ*DVz;
452  Jm[5][3]=CQ*(dsdVx*DVz+s*gB[1]);
453  Jm[5][4]=CQ*(dsdVy*DVz-s*gB[0]);
454  Jm[5][5]=1.0+dsdVz*CQ*DVz;
455  Jm[5][6]=(CQ*dsdQ+C*s)*DVz;
456 
457  Jm[6][6]=1.0;
458 
459  memset(&J1[0][0],0,sizeof(J1));
460 
461  J1[0][0]=M[0][0];J1[0][1]=M[0][1];J1[0][2]=M[0][2];
462  J1[1][0]=M[1][0];J1[1][1]=M[1][1];J1[1][2]=M[1][2];
463  J1[2][3]=-V[1]/(V[0]*V[0]+V[1]*V[1]);
464  J1[2][4]= V[0]/(V[0]*V[0]+V[1]*V[1]);
465  J1[3][5]=-1.0/sqrt(1-V[2]*V[2]);
466  J1[4][6]=1.0;
467 
468  for(i=0;i<7;i++)
469  {
470  for(j=0;j<2;j++)
471  Buf[j][i]=J1[j][0]*Jm[0][i]+J1[j][1]*Jm[1][i]+J1[j][2]*Jm[2][i];
472  Buf[2][i]=J1[2][3]*Jm[3][i]+J1[2][4]*Jm[4][i];
473  Buf[3][i]=J1[3][5]*Jm[5][i];
474  Buf[4][i]=Jm[6][i];
475  }
476 
477  if(pSB!=nullptr)
478  {
479  for(i=0;i<5;i++)
480  {
481  J[i][0]=Buf[i][0]*J0[0][0]+Buf[i][1]*J0[1][0]+Buf[i][2]*J0[2][0];
482  J[i][1]=Buf[i][0]*J0[0][1]+Buf[i][1]*J0[1][1]+Buf[i][2]*J0[2][1];
483  J[i][2]=Buf[i][3]*J0[3][2]+Buf[i][4]*J0[4][2];
484  J[i][3]=Buf[i][3]*J0[3][3]+Buf[i][4]*J0[4][3]+Buf[i][5]*J0[5][3];
485  J[i][4]=Buf[i][6];
486  }
487  }
488  else
489  {
490  for(i=0;i<5;i++)
491  {
492  J[i][0]=Buf[i][0]*J0[0][0]+Buf[i][1]*J0[1][0];
493  J[i][1]=Buf[i][2];
494  J[i][2]=Buf[i][0]*J0[0][2]+Buf[i][1]*J0[1][2]+Buf[i][3]*J0[3][2]+Buf[i][4]*J0[4][2];
495  J[i][3]=Buf[i][3]*J0[3][3]+Buf[i][4]*J0[4][3]+Buf[i][5]*J0[5][3];
496  J[i][4]=Buf[i][6];
497  }
498  }
499  }
500  else {
501  Rf[0]=pTS->getTrackState(0);
502  Rf[1]=pTS->getTrackState(1);
503  Rf[2]=pTS->getTrackState(2);
504  Rf[3]=pTS->getTrackState(3);
505  Rf[4]=pTS->getTrackState(4);
506  memset(&J[0][0],0,sizeof(J));
507  for(i=0;i<5;i++) J[i][i]=1.0;
508  }
509 
510  for(i=0;i<5;i++) for(j=0;j<5;j++)
511  {
512  AG[i][j]=0.0;for(m=0;m<5;m++) AG[i][j]+=J[i][m]*pTS->getTrackCovariance(m,j);
513  }
514  for(i=0;i<5;i++) for(j=i;j<5;j++)
515  {
516  Gf[i][j]=0.0;
517  for(m=0;m<5;m++) Gf[i][j]+=AG[i][m]*J[j][m];
518  Gf[j][i]=Gf[i][j];
519  }
520 
521  Trk::TrkTrackState* pTE=new Trk::TrkTrackState(pTS);
522 
523  //workaround to keep using existing TrkTrackState code
524  double Rtmp[5];
525 
526  for(i=0;i<4;i++) Rtmp[i] = Rf[i];
527  Rtmp[4] = 0.001*Rf[4];//GeV->MeV
528 
529  pTE->setTrackState(Rtmp);
530  pTE->setTrackCovariance(Gf);
531  pTE->attachToSurface(pSE);
532 
533  if(m_doMultScatt)
534  pTE->applyMultipleScattering();
535 
536  pTE->setTrackState(Rf);//restore
537 
538  if(m_doBremm)
539  pTE->applyEnergyLoss(1);
540 
541  //quick check for covariance sanity
542 
543  for(int idx=0;idx<5;idx++) {
544  if(pTE->getTrackCovariance(idx,idx) < 0) {
545  ATH_MSG_DEBUG("REGTEST: cov(" << idx << "," << idx << ") =" << pTE->getTrackCovariance(idx,idx) << " < 0, reject track");
546  delete pTE;
547  return nullptr;
548  }
549  }
550 
551  AmgSymMatrix(5) Gi;
552  for(i=0;i<5;i++) for(j=i;j<5;j++)
553  {
554  Gi.fillSymmetric(i, j, pTE->getTrackCovariance(i,j));
555  }
556  Gi = Gi.inverse();
557 
558  for(i=0;i<5;i++) for(j=0;j<5;j++)
559  {
560  A[i][j]=0.0;
561  for(m=0;m<5;m++) A[i][j]+=AG[m][i]*Gi(m,j);
562  }
563  pTE->setPreviousState(pTS);
564  pTE->setSmootherGain(A);
565 
566  return pTE;
567 }
568 
569 void TrigInDetTrackFitter::fit(const TrackCollection& inputTracks, TrackCollection& fittedTracks, const EventContext& ctx, const Trk::ParticleHypothesis& matEffects) const
570 {
571  TrackCollection tmp;
572  fit(inputTracks,fittedTracks,tmp,ctx,matEffects,false);
573 }
574 
575 void TrigInDetTrackFitter::fit(const TrackCollection& inputTracks, TrackCollection& fittedTracks, TrackCollection& fittedTrackswTP, const EventContext& ctx, const Trk::ParticleHypothesis& matEffects, const bool addTPtoTSoS) const
576 {
577  MagField::AtlasFieldCache fieldCache;
578 
579  SG::ReadCondHandle<AtlasFieldCacheCondObj> fieldCondObj{m_fieldCondObjInputKey, ctx};
580  if (!fieldCondObj.isValid()) {
581  ATH_MSG_ERROR("Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
582  return;
583  }
584 
585  fieldCondObj->getInitializedCache (fieldCache);
586  fittedTracks.reserve(inputTracks.size());
587  if( addTPtoTSoS ) fittedTrackswTP.reserve(inputTracks.size());
588  for(auto trIt = inputTracks.begin(); trIt != inputTracks.end(); ++trIt) {
589  Trk::Track* fittedTrack = nullptr;
590  Trk::Track* fittedTrackwTP = nullptr;
591  std::tie(fittedTrack,fittedTrackwTP) = fitTrack(**trIt, fieldCache, ctx, matEffects, addTPtoTSoS);
592  if (fittedTrack!=nullptr) {
593  fittedTracks.push_back(fittedTrack);
594  }
595  if (addTPtoTSoS && fittedTrackwTP!=nullptr) {
596  fittedTrackswTP.push_back(fittedTrackwTP);
597  }
598  }
599 }
600 
601 std::pair<Trk::Track*,Trk::Track*> TrigInDetTrackFitter::fitTrack(const Trk::Track& recoTrack, MagField::AtlasFieldCache& fieldCache, const EventContext& ctx, const Trk::ParticleHypothesis& matEffects, const bool addTPtoTSoS) const {
602 
603  const Trk::TrackParameters* trackPars = recoTrack.perigeeParameters();
604  if(trackPars==nullptr) {
605  ATH_MSG_WARNING("Fit Failed -- TrkTrack has no parameters");
606  return std::make_pair(nullptr,nullptr);
607  }
608 
609  // 1. Create initial track state:
610  double Rk[5];
611  Rk[0] = trackPars->parameters()[Trk::d0];
612  Rk[1] = trackPars->parameters()[Trk::z0];
613  Rk[2] = trackPars->parameters()[Trk::phi0];
614  if(Rk[2]>M_PI) Rk[2]-=2*M_PI;
615  if(Rk[2]<-M_PI) Rk[2]+=2*M_PI;
616  double trk_theta = trackPars->parameters()[Trk::theta];
617  Rk[3] = trk_theta;
618  double trk_qOverP = trackPars->parameters()[Trk::qOverP];
619  Rk[4] = 1000.0*trk_qOverP;//MeV->GeV
620  double trk_Pt = sin(trk_theta)/trk_qOverP;
621 
622  if(fabs(trk_Pt)<100.0)
623  {
624  ATH_MSG_DEBUG("Estimated Pt is too low "<<trk_Pt<<" - skipping fit");
625  return std::make_pair(nullptr,nullptr);
626  }
627 
628  // 2. Create filtering nodes
629 
630  std::vector<Trk::TrkBaseNode*> vpTrkNodes;
631  std::vector<Trk::TrkTrackState*> vpTrackStates;
632  vpTrackStates.reserve(vpTrkNodes.size() + 1);
633  bool trackResult = m_trackMaker->createDkfTrack(recoTrack,vpTrkNodes, m_DChi2);
634  int nHits=vpTrkNodes.size();
635  ATH_MSG_VERBOSE(nHits<<" filtering nodes created");
636 
637  if(!trackResult) return std::make_pair(nullptr,nullptr);
638 
639  // 3. Main algorithm: filter and smoother (Rauch-Tung-Striebel)
640  m_nTracksTotal++;
641  Trk::TrkTrackState* pTS = new Trk::TrkTrackState(Rk);
642  double Gk[5][5] = {{100.0, 0, 0, 0, 0},
643  {0, 100.0, 0, 0, 0},
644  {0, 0, 0.01, 0, 0},
645  {0, 0, 0, 0.01, 0},
646  {0, 0, 0, 0, 0.1}};
647  pTS->setTrackCovariance(Gk);
648  if(m_doMultScatt)
649  pTS->setScatteringMode(1);
650  if(m_doBremm)
651  pTS->setScatteringMode(2);
652  vpTrackStates.push_back(pTS);
653 
654  //ATH_MSG_DEBUG("Initial chi2: "<<recoTrack.chi2()<<" track authorId: "<<recoTrack.algorithmId());
655  ATH_MSG_VERBOSE("Initial params: locT="<<Rk[0]<<" locL="<<Rk[1]<<" phi="<<Rk[2]
656  <<" theta="<<Rk[3]<<" Q="<<Rk[4]<<" pT="<<sin(Rk[3])/Rk[4]<<" GeV");
657 
658  bool OK=true;
659 
660  double chi2tot=0.0;
661  int ndoftot=-5;
662 
663  Trk::TrkPlanarSurface* pSB=nullptr;
664  Trk::TrkPlanarSurface* pSE=nullptr;
665  for(auto pnIt = vpTrkNodes.begin(); pnIt!=vpTrkNodes.end(); ++pnIt) {
666  pSE=(*pnIt)->getSurface();
667  Trk::TrkTrackState* pNS=extrapolate(pTS,pSB,pSE,fieldCache);
668 
669  pSB=pSE;
670  if(pNS!=nullptr) {
671  vpTrackStates.push_back(pNS);
672 
673  (*pnIt)->validateMeasurement(pNS);
674  ATH_MSG_VERBOSE("dChi2="<<(*pnIt)->getChi2());
675  if((*pnIt)->isValidated())
676  {
677  chi2tot+=(*pnIt)->getChi2();
678  ndoftot+=(*pnIt)->getNdof();
679  }
680  (*pnIt)->updateTrackState(pNS);
681  pTS=pNS;
682  double est_Pt = 1000.0*sin(pTS->getTrackState(3))/pTS->getTrackState(4);
683  if(fabs(est_Pt)<200.0)
684  {
685  ATH_MSG_VERBOSE("Estimated Pt is too low "<<est_Pt<<" - skipping fit");
686  m_fitErrorsLowPt++;
687  OK=false;break;
688  }
689  }
690  else
691  {
692  m_fitErrorsDivergence++;
693  OK=false;break;
694  }
695  }
696  Trk::Track* fittedTrack = nullptr;
697  Trk::Track* fittedTrackwTP = nullptr;
698  if(OK)
699  {
700  for(auto ptsIt = vpTrackStates.rbegin();ptsIt!=vpTrackStates.rend();++ptsIt)
701  {
702  (*ptsIt)->runSmoother();
703  }
704  pTS=(*vpTrackStates.begin());
705  //correct GeV->MeV
706 
707  correctScale(pTS);
708 
709  double qOverP = pTS->getTrackState(4);
710  double pt=sin(pTS->getTrackState(3))/pTS->getTrackState(4);
711  double phi0 = pTS->getTrackState(2);
712  if(phi0>M_PI) phi0-=2*M_PI;
713  if(phi0<-M_PI) phi0+=2*M_PI;
714  double theta = pTS->getTrackState(3);
715  double z0 = pTS->getTrackState(1);
716  double d0 = pTS->getTrackState(0);
717  bool bad_cov = false;
718  auto cov = AmgSymMatrix(5){};
719  for(int i=0;i<5;i++) {
720  double cov_diag = pTS->getTrackCovariance(i,i);
721  if (cov_diag < 0) {
722  bad_cov = true;//Diagonal elements must be positive
723  ATH_MSG_DEBUG("REGTEST: cov(" << i << "," << i << ") =" << cov_diag << " < 0, reject track");
724  break;
725  }
726  (cov)(i, i) = pTS->getTrackCovariance(i,i);
727  for(int j=i+1;j<5;j++) {
728  cov.fillSymmetric(i, j, pTS->getTrackCovariance(i,j));
729  }
730  }
731 
732  if((ndoftot<0) || (fabs(pt)<100.0) || (std::isnan(pt)) || bad_cov)
733  {
734  ATH_MSG_DEBUG("Fit failed - possibly floating point problem");
735  }
736  else
737  {
738  Trk::PerigeeSurface perigeeSurface;
739  auto perigee = std::make_unique<Trk::Perigee>(d0, z0, phi0, theta, qOverP, perigeeSurface, cov);
740  ATH_MSG_VERBOSE("perigee: " << *perigee);
741  std::unique_ptr<Trk::Perigee> perigeewTP = (addTPtoTSoS) ? std::make_unique<Trk::Perigee>(d0, z0, phi0, theta, qOverP, perigeeSurface, cov) : nullptr;
742 
743  std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern;
744  typePattern.set(Trk::TrackStateOnSurface::Perigee);
745  auto pParVec = std::make_unique<Trk::TrackStates>();
746  auto pParVecwTP = std::make_unique<Trk::TrackStates>();
747  if (m_correctClusterPos) {
748  pParVec->reserve(vpTrkNodes.size()+1);
749  pParVec->push_back(new Trk::TrackStateOnSurface(nullptr, std::move(perigee),nullptr, typePattern));
750  if (addTPtoTSoS) {
751  std::bitset<
752  Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes>
753  typePatternwTP;
754  typePatternwTP.set(Trk::TrackStateOnSurface::Perigee);
755  pParVecwTP->reserve(vpTrkNodes.size() + 1);
756  pParVecwTP->push_back(
757  new Trk::TrackStateOnSurface(nullptr, std::move(perigeewTP), nullptr, typePatternwTP));
758  }
759  for (auto pnIt = vpTrkNodes.begin(); pnIt != vpTrkNodes.end(); ++pnIt) {
760  if((*pnIt)->isValidated()) {
761  Trk::TrackStateOnSurface* pTSS = createTrackStateOnSurface(*pnIt,false, ctx);
762  Trk::TrackStateOnSurface* pTSSwTP = nullptr;
763  if( addTPtoTSoS ) pTSSwTP = createTrackStateOnSurface(*pnIt,true,ctx);
764  if(pTSS!=nullptr) {
765  pParVec->push_back(pTSS);
766  }
767  if(pTSSwTP!=nullptr) {
768  pParVecwTP->push_back(pTSSwTP);
769  }
770  }
771  }
772  }
773  else {
774  pParVec->reserve(recoTrack.trackStateOnSurfaces()->size());
775  pParVec->push_back(new Trk::TrackStateOnSurface(nullptr, std::move(perigee),nullptr, typePattern));
776 
777  for (auto tSOS = recoTrack.trackStateOnSurfaces()->begin(); tSOS != recoTrack.trackStateOnSurfaces()->end(); ++tSOS) {
778  //Don't store perigee - new perigee created above
779  if ((*tSOS)->type(Trk::TrackStateOnSurface::Perigee) == false) {
780  pParVec->push_back((*tSOS)->clone());
781  }
782  }
783  }
784  ATH_MSG_VERBOSE("Total chi2 ="<<chi2tot<<" NDOF="<<ndoftot);
785  if(msgLvl(MSG::VERBOSE)) {
786  double eta = -log(tan(0.5*theta));
787  ATH_MSG_VERBOSE("Fitted parameters: d0="<<d0<<" phi0="<<phi0<<" z0="<<z0
788  <<" eta0="<<eta<<" pt="<<pt);
789  }
790  auto pFQ= std::make_unique<Trk::FitQuality>(chi2tot,ndoftot);
791  Trk::TrackInfo info(recoTrack.info());
792  info.setParticleHypothesis(matEffects);
793  fittedTrack = new Trk::Track(info, std::move(pParVec), std::move(pFQ));//fittedTrack now owns pParVec and pFQ
794  if( addTPtoTSoS ) {
795  auto pFQwTP=std::make_unique<Trk::FitQuality>(chi2tot,ndoftot);
796  Trk::TrackInfo infowTP(recoTrack.info());
797  infowTP.setParticleHypothesis(matEffects);
798  fittedTrackwTP = new Trk::Track(infowTP, std::move(pParVecwTP), std::move(pFQwTP));//fittedTrack now owns pParVecwTP and pFQwTP
799  }
800  }
801  }
802  else
803  {
804  ATH_MSG_DEBUG("Forward Kalman filter: extrapolation failure ");
805  }
806 
807  for(auto pnIt = vpTrkNodes.begin(); pnIt!=vpTrkNodes.end(); ++pnIt) {
808  delete((*pnIt)->getSurface());
809  delete (*pnIt);
810  }
811  vpTrkNodes.clear();
812  for(auto ptsIt = vpTrackStates.begin();ptsIt!=vpTrackStates.end();++ptsIt) {
813  delete (*ptsIt);
814  }
815  vpTrackStates.clear();
816 
817  return std::make_pair(fittedTrack,fittedTrackwTP);
818 }
819 
820 Trk::TrackStateOnSurface* TrigInDetTrackFitter::createTrackStateOnSurface(Trk::TrkBaseNode* pN, const bool addTPtoTSoS, const EventContext& ctx) const
821 {
822  Trk::TrackStateOnSurface* pTSS=nullptr;
823  char type=pN->getNodeType();
824  std::unique_ptr<Trk::TrackParameters> pTP{};
825  if(type==0) return pTSS;
826 
827 
828  Trk::TrkTrackState* pTS=pN->getTrackState();
829  auto pM = AmgSymMatrix(5){};
830  for(int i=0;i<5;i++) {
831  for(int j=0;j<5;j++) {
832  (pM)(i,j)=pTS->getTrackCovariance(i,j);
833  }
834  }
835  const Trk::PrepRawData* pPRD=pN->getPrepRawData();
836 
837  if((type==1)||(type==2))
838  {
839  const Trk::Surface& rS = pPRD->detectorElement()->surface();
840  const Trk::PlaneSurface* pPS = dynamic_cast<const Trk::PlaneSurface*>(&rS);
841  if(pPS==nullptr) return pTSS;
842 
843  pTP=std::make_unique<Trk::AtaPlane>(pTS->getTrackState(0),
844  pTS->getTrackState(1),
845  pTS->getTrackState(2),
846  pTS->getTrackState(3),
847  pTS->getTrackState(4),*pPS,
848  std::move(pM));
849  }
850  else if(type==3)
851  {
852  const Trk::Surface& rS = pPRD->detectorElement()->surface(pPRD->identify());
853  const Trk::StraightLineSurface* pLS=dynamic_cast<const Trk::StraightLineSurface*>(&rS);
854  if(pLS==nullptr) return pTSS;
855 
856  if((pTS->getTrackState(2)<-M_PI) ||(pTS->getTrackState(2)>M_PI)) {
857  ATH_MSG_WARNING("Phi out of range when correcting Trk::TrackStateOnSurface");
858  }
859 
860 
861  pTP=std::make_unique<Trk::AtaStraightLine>(pTS->getTrackState(0),
862  pTS->getTrackState(1),
863  pTS->getTrackState(2),
864  pTS->getTrackState(3),
865  pTS->getTrackState(4),
866  *pLS,
867  std::move(pM));
868  }
869  if(pTP==nullptr) return nullptr;
870  std::unique_ptr<Trk::RIO_OnTrack> pRIO{m_ROTcreator->correct(*pPRD,*pTP,ctx)};
871  if(pRIO==nullptr) {
872  return nullptr;
873  }
874  std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern;
875  typePattern.set(Trk::TrackStateOnSurface::Measurement);
876  typePattern.set(Trk::TrackStateOnSurface::Scatterer);
877  auto pFQ=Trk::FitQualityOnSurface(pN->getChi2(),pN->getNdof());
878  if( addTPtoTSoS ) {
879  std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePatternwTP;
880  typePatternwTP.set(Trk::TrackStateOnSurface::Measurement);
881  typePatternwTP.set(Trk::TrackStateOnSurface::Scatterer);
882  auto pFQwTP=Trk::FitQualityOnSurface(pN->getChi2(),pN->getNdof());
883  pTSS = new Trk::TrackStateOnSurface(pFQwTP, pRIO->uniqueClone(), std::move(pTP), nullptr, typePatternwTP);
884  }
885  else {
886  pTSS = new Trk::TrackStateOnSurface(pFQ, std::move(pRIO), nullptr, nullptr, typePattern);
887  }
888  return pTSS;
889 }
890 
891 StatusCode TrigInDetTrackFitter::getUnbiasedResiduals(const Trk::Track& pT,
892  std::vector<TrigL2HitResidual>& vResid, const EventContext& ctx) const {
893 
894  const Trk::TrackParameters* trackPars = pT.perigeeParameters();
895  if(trackPars==nullptr) {
896  ATH_MSG_WARNING("Fit Failed -- TrkTrack has no parameters");
897  return StatusCode::FAILURE;
898  }
899 
900  MagField::AtlasFieldCache fieldCache;
901 
902  SG::ReadCondHandle<AtlasFieldCacheCondObj> fieldCondObj{m_fieldCondObjInputKey, ctx};
903  if (!fieldCondObj.isValid()) {
904  ATH_MSG_ERROR("Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
905  return StatusCode::FAILURE;
906  }
907 
908  fieldCondObj->getInitializedCache (fieldCache);
909  std::vector<Trk::TrkBaseNode*> vpTrkNodes;
910  std::vector<Trk::TrkTrackState*> vpTrackStates;
911  vResid.clear();
912  double trk_theta = trackPars->parameters()[Trk::theta];
913  double trk_qOverP = trackPars->parameters()[Trk::qOverP];
914  double Pt = sin(trk_theta)/trk_qOverP;
915  if(fabs(Pt)<100.0)
916  {
917  ATH_MSG_DEBUG("TrigL2ResidualCalculator failed -- Estimated Pt is too low "<<Pt);
918  return StatusCode::FAILURE;
919  }
920 
921  // 1. Create filtering nodes
922 
923  bool trackResult = m_trackMaker->createDkfTrack(pT,vpTrkNodes, m_DChi2);
924  if(!trackResult)
925  {
926  ATH_MSG_DEBUG("TrigDkfTrackMaker failed");
927  return StatusCode::FAILURE;
928  }
929 
930  bool OK=true;
931  std::vector<Trk::TrkBaseNode*>::iterator pnIt,pnEnd(vpTrkNodes.end());
932 
933  for(std::vector<Trk::TrkBaseNode*>::iterator pNodeIt=vpTrkNodes.begin();pNodeIt!=vpTrkNodes.end();
934  ++pNodeIt)
935  {
936  Trk::TrkBaseNode* pMaskedNode=(*pNodeIt);
937  Trk::TrkTrackState* pMaskedState=nullptr;
938 
939  // 2. Create initial track state:
940 
941  double Rk[5];
942  Rk[0] = trackPars->parameters()[Trk::d0];
943  Rk[1] = trackPars->parameters()[Trk::z0];
944  Rk[2] = trackPars->parameters()[Trk::phi0];
945  if(Rk[2]>M_PI) Rk[2]-=2*M_PI;
946  if(Rk[2]<-M_PI) Rk[2]+=2*M_PI;
947  trk_theta = trackPars->parameters()[Trk::theta];
948  Rk[3] = trk_theta;
949  Rk[4] = 1000.0*trackPars->parameters()[Trk::qOverP];//MeV->GeV
950  //No need to correct scale back - not returning track
951 
952  // 3. Main algorithm: filter and smoother (Rauch-Tung-Striebel)
953 
954  Trk::TrkTrackState* pTS = new Trk::TrkTrackState(Rk);
955  double Gk[5][5] = {{100.0, 0, 0, 0, 0},
956  {0, 100.0, 0, 0, 0},
957  {0, 0, 0.01, 0, 0},
958  {0, 0, 0, 0.01, 0},
959  {0, 0, 0, 0, 0.1}};
960  pTS->setTrackCovariance(Gk);
961  if(m_doMultScatt)
962  pTS->setScatteringMode(1);
963  if(m_doBremm)
964  pTS->setScatteringMode(2);
965  vpTrackStates.push_back(pTS);
966 
967  ATH_MSG_DEBUG("Initial params: locT="<<Rk[0]<<" locL="<<Rk[1]<<" phi="<<Rk[2]
968  <<" theta="<<Rk[3]<<" Q="<<Rk[4]<<" pT="<<sin(Rk[3])/Rk[4]);
969 
970  OK=true;
971  Trk::TrkPlanarSurface *pSB=nullptr,*pSE;
972 
973  for(pnIt=vpTrkNodes.begin();pnIt!=pnEnd;++pnIt)
974  {
975  pSE=(*pnIt)->getSurface();
976  Trk::TrkTrackState* pNS=extrapolate(pTS,pSB,pSE,fieldCache);
977 
978  pSB=pSE;
979  if(pNS!=nullptr)
980  {
981  vpTrackStates.push_back(pNS);
982 
983  (*pnIt)->validateMeasurement(pNS);
984  ATH_MSG_DEBUG("dChi2="<<(*pnIt)->getChi2());
985  if((*pnIt)!=pMaskedNode)
986  {
987  (*pnIt)->updateTrackState(pNS);
988  }
989  else
990  {
991  pMaskedState=pNS;
992  }
993  pTS=pNS;
994  Pt=sin(pTS->getTrackState(3))/pTS->getTrackState(4);
995  if(fabs(Pt)<0.2)
996  {
997  ATH_MSG_DEBUG("Estimated Pt is too low "<<Pt<<" - skipping fit");
998  OK=false;break;
999  }
1000  }
1001  else
1002  {
1003  OK=false;break;
1004  }
1005  }
1006  if(OK)
1007  {
1008  std::vector<Trk::TrkTrackState*>::reverse_iterator ptsrIt(vpTrackStates.rbegin()),
1009  ptsrEnd(vpTrackStates.rend());
1010 
1011  for(;ptsrIt!=ptsrEnd;++ptsrIt)
1012  {
1013  (*ptsrIt)->runSmoother();
1014  }
1015 
1016  pMaskedNode->validateMeasurement(pMaskedState);
1017 
1018  double r[2],V[2][2];
1019 
1020  int nSize=pMaskedNode->getResiduals(r);
1021  nSize=pMaskedNode->getInverseResidualVariance(V);
1022  const Trk::PrepRawData* pPRD = pMaskedNode->getPrepRawData();
1023 
1024  Identifier id = pPRD->identify();
1025 
1026  Region region = Region::Undefined;
1027  if(m_idHelper->is_pixel(id))
1028  {
1029  region=(m_pixelId->is_barrel(id)) ? Region::PixBarrel: Region::PixEndcap;
1030  if (m_pixelId->is_blayer(id)) {
1031  region = Region::IBL;
1032  }
1033  }
1034  if(m_idHelper->is_sct(id))
1035  {
1036  region=(m_sctId->is_barrel(id)) ? Region::SctBarrel : Region::SctEndcap;
1037  }
1038  if(nSize==1) {
1039  if(V[0][0]>0.0) {
1040  vResid.push_back(TrigL2HitResidual(id,region,r[0],r[0]*sqrt(V[0][0])));
1041  }
1042  else {
1043  OK=false;
1044  break;
1045  }
1046  }
1047  else {
1048  if((V[0][0]>0.0) && (V[1][1]>0.0)) {
1049  vResid.push_back(TrigL2HitResidual(id,region,r[0],r[0]*sqrt(V[0][0]),
1050  r[1],r[1]*sqrt(V[1][1])));
1051  }
1052  else {
1053  OK=false;
1054  break;
1055  }
1056  }
1057  }
1058  else
1059  {
1060  ATH_MSG_DEBUG("Forward Kalman filter: extrapolation failure ");
1061  vResid.clear();
1062  }
1063  for(std::vector<Trk::TrkTrackState*>::iterator ptsIt=vpTrackStates.begin();
1064  ptsIt!=vpTrackStates.end();++ptsIt) delete (*ptsIt);
1065  vpTrackStates.clear();
1066  if(!OK) break;
1067  }
1068  pnIt=vpTrkNodes.begin();pnEnd=vpTrkNodes.end();
1069  for(;pnIt!=pnEnd;++pnIt)
1070  {
1071  delete((*pnIt)->getSurface());
1072  delete (*pnIt);
1073  }
1074  vpTrkNodes.clear();
1075 
1076  if(OK) return StatusCode::SUCCESS;
1077  else return StatusCode::FAILURE;
1078 
1079 }
TrigInDetTrackFitter::m_idHelper
const AtlasDetectorID * m_idHelper
Definition: TrigInDetTrackFitter.h:71
DataVector::reserve
void reserve(size_type n)
Attempt to preallocate enough memory for a specified number of elements.
grepfile.info
info
Definition: grepfile.py:38
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
PixelID.h
This is an Identifier helper class for the Pixel subdetector. This class is a factory for creating co...
AtlasDetectorID::is_pixel
bool is_pixel(Identifier id) const
Definition: AtlasDetectorID.h:760
ITrigInDetTrackFitter.h
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
beamspotman.r
def r
Definition: beamspotman.py:676
cost
int cost(std::vector< std::string > &files, node &n, const std::string &directory="", bool deleteref=false, bool relocate=false)
Definition: hcg.cxx:921
Region::SctBarrel
@ SctBarrel
TrigInDetTrackFitter::fitTrack
std::pair< Trk::Track *, Trk::Track * > fitTrack(const Trk::Track &, MagField::AtlasFieldCache &, const EventContext &ctx, const Trk::ParticleHypothesis &matEffects=Trk::pion, const bool addTPtoTSoS=false) const
Definition: TrigInDetTrackFitter.cxx:601
CalculateHighPtTerm.pT
pT
Definition: ICHEP2016/CalculateHighPtTerm.py:57
Trk::TrackInfo
Contains information about the 'fitter' of this track.
Definition: Tracking/TrkEvent/TrkTrack/TrkTrack/TrackInfo.h:32
Trk::TrackStateOnSurface::Perigee
@ Perigee
This represents a perigee, and so will contain a Perigee object only.
Definition: TrackStateOnSurface.h:117
checkxAOD.ds
ds
Definition: Tools/PyUtils/bin/checkxAOD.py:260
Trk::TrkTrackState::applyMultipleScattering
void applyMultipleScattering()
Definition: TrkTrackState.cxx:158
Trk::TrkPlanarSurface::transformPointToLocal
void transformPointToLocal(const double *, double *)
Definition: TrkPlanarSurface.cxx:113
SCT_ID.h
This is an Identifier helper class for the SCT subdetector. This class is a factory for creating comp...
Region::PixBarrel
@ PixBarrel
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
TrackParameters.h
AthCheckMacros.h
SG::ReadCondHandle
Definition: ReadCondHandle.h:44
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
AthMsgStreamMacros.h
Surface.h
Trk::TrkPlanarSurface::transformPointToGlobal
void transformPointToGlobal(const double *, double *)
Definition: TrkPlanarSurface.cxx:122
Trk::Track
The ATLAS Track class.
Definition: Tracking/TrkEvent/TrkTrack/TrkTrack/Track.h:73
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
TrigInDetTrackFitter::m_DChi2
double m_DChi2
Definition: TrigInDetTrackFitter.h:62
Trk::PerigeeSurface
Definition: PerigeeSurface.h:43
AtlasDetectorID::is_sct
bool is_sct(Identifier id) const
Definition: AtlasDetectorID.h:770
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:83
AtlasFieldCacheCondObj.h
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
Trk::Track::trackStateOnSurfaces
const Trk::TrackStates * trackStateOnSurfaces() const
return a pointer to a const DataVector of const TrackStateOnSurfaces.
TrigInDetTrackFitter::m_nTracksTotal
std::atomic< size_t > m_nTracksTotal
Definition: TrigInDetTrackFitter.h:57
Trk::Track::info
const TrackInfo & info() const
Returns a const ref to info of a const tracks.
DMTest::P
P_v1 P
Definition: P.h:23
Region::SctEndcap
@ SctEndcap
IRIO_OnTrackCreator.h
InDetAccessor::phi0
@ phi0
Definition: InDetAccessor.h:33
TrigInDetTrackFitter::m_fieldCondObjInputKey
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_fieldCondObjInputKey
Definition: TrigInDetTrackFitter.h:68
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:75
Trk::TrkTrackState
Definition: Tracking/TrkFitter/TrkDistributedKalmanFilter/TrkDistributedKalmanFilter/TrkTrackState.h:24
TrigInDetTrackFitter::initialize
virtual StatusCode initialize()
Definition: TrigInDetTrackFitter.cxx:76
DMTest::C
C_v1 C
Definition: C.h:26
xAOD::JetInput::Track
@ Track
Definition: JetContainerInfo.h:61
Trk::TrkPlanarSurface::getRotMatrix
double getRotMatrix(int, int)
Definition: TrkPlanarSurface.cxx:84
Trk::TrkPlanarSurface
Definition: Tracking/TrkFitter/TrkDistributedKalmanFilter/TrkDistributedKalmanFilter/TrkPlanarSurface.h:25
python.SystemOfUnits.kilogauss
int kilogauss
Definition: SystemOfUnits.py:231
plotBeamSpotVxVal.cov
cov
Definition: plotBeamSpotVxVal.py:201
test_pyathena.pt
pt
Definition: test_pyathena.py:11
M_PI
#define M_PI
Definition: ActiveFraction.h:11
TrigInDetTrackFitter::fit
void fit(const TrackCollection &, TrackCollection &, const EventContext &, const Trk::ParticleHypothesis &matEffects=Trk::pion) const
Definition: TrigInDetTrackFitter.cxx:569
mc.diff
diff
Definition: mc.SFGenPy8_MuMu_DD.py:14
TrigInDetTrackFitter::m_sctId
const SCT_ID * m_sctId
Definition: TrigInDetTrackFitter.h:70
TrigInDetTrackFitter::m_doBremm
bool m_doBremm
Definition: TrigInDetTrackFitter.h:64
Trk::z0
@ z0
Definition: ParamDefs.h:64
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
Trk::TrkBaseNode::validateMeasurement
virtual void validateMeasurement(TrkTrackState *)=0
FitQualityOnSurface.h
TrigInDetTrackCollection.h
ReadOfcFromCool.field
field
Definition: ReadOfcFromCool.py:48
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
Trk::TrkBaseNode::getNodeType
virtual char getNodeType()
Definition: TrkBaseNode.cxx:61
SG::VarHandleKey::key
const std::string & key() const
Return the StoreGate ID for the referenced object.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:141
InDetAccessor::qOverP
@ qOverP
perigee
Definition: InDetAccessor.h:35
PixelID::is_blayer
bool is_blayer(const Identifier &id) const
Test for b-layer - WARNING: id MUST be pixel id, otherwise answer is not accurate....
Definition: PixelID.h:633
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition: EventPrimitives.h:50
PrepRawData.h
AthCommonDataStore< AthCommonMsg< AlgTool > >::detStore
const ServiceHandle< StoreGateSvc > & detStore() const
The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.
Definition: AthCommonDataStore.h:95
Trk::TrkBaseNode::getNdof
int getNdof() const
Definition: TrkBaseNode.cxx:65
Region::PixEndcap
@ PixEndcap
Track.h
TrigInDetTrackFitter::m_ROTcreator
ToolHandle< Trk::IRIO_OnTrackCreator > m_ROTcreator
Definition: TrigInDetTrackFitter.h:67
Trk::FitQualityOnSurface
Definition: FitQualityOnSurface.h:19
Trk::ParticleHypothesis
ParticleHypothesis
Definition: ParticleHypothesis.h:25
TrigInDetTrack.h
AtlasDetectorID.h
This class provides an interface to generate or decode an identifier for the upper levels of the dete...
TrigInDetTrackFitter::finalize
virtual StatusCode finalize()
Definition: TrigInDetTrackFitter.cxx:91
PixelID::is_barrel
bool is_barrel(const Identifier &id) const
Test for barrel - WARNING: id MUST be pixel id, otherwise answer is not accurate. Use SiliconID for g...
Definition: PixelID.h:605
A
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
TrkTrackState.h
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
python.TrigInDetConfig.inputTracks
inputTracks
Definition: TrigInDetConfig.py:183
lumiFormat.i
int i
Definition: lumiFormat.py:85
beamspotman.n
n
Definition: beamspotman.py:731
TrkBaseNode.h
Trk::theta
@ theta
Definition: ParamDefs.h:66
Trk::TrkBaseNode
Definition: Tracking/TrkFitter/TrkDistributedKalmanFilter/TrkDistributedKalmanFilter/TrkBaseNode.h:23
Trk::TrkTrackState::setScatteringMode
void setScatteringMode(int)
Definition: TrkTrackState.cxx:130
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
TrigL2HitResidual.h
TRT::Track::d0
@ d0
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:62
Trk::TrkDetElementBase::surface
virtual const Surface & surface() const =0
Return surface associated with this detector element.
TrigInDetTrackFitter::correctScale
void correctScale(Trk::TrkTrackState *) const
Definition: TrigInDetTrackFitter.cxx:111
Region
Region
Definition: TrigL2HitResidual.h:14
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Region::IBL
@ IBL
Trk::TrkPlanarSurface::getPar
double getPar(int)
Definition: TrkPlanarSurface.cxx:80
python.TransformConfig.descr
descr
print "%s.properties()" % self.__name__
Definition: TransformConfig.py:360
TrigInDetTrackFitter::m_pixelId
const PixelID * m_pixelId
Definition: TrigInDetTrackFitter.h:69
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
Trk::ParametersBase
Definition: ParametersBase.h:55
Trk::TrkPlanarSurface::getInvRotMatrix
double getInvRotMatrix(int, int)
Definition: TrkPlanarSurface.cxx:88
TRT::Track::z0
@ z0
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:63
DeMoUpdate.tmp
string tmp
Definition: DeMoUpdate.py:1167
Trk::TrkTrackState::attachToSurface
void attachToSurface(TrkPlanarSurface *)
Definition: TrkTrackState.cxx:101
DataVector< Trk::Track >
Trk::TrkBaseNode::getTrackState
TrkTrackState * getTrackState()
Definition: TrkBaseNode.cxx:43
ParticleHypothesis.h
Trk::TrkTrackState::applyEnergyLoss
void applyEnergyLoss(int)
Definition: TrkTrackState.cxx:184
Trk::TrkTrackState::setPreviousState
void setPreviousState(TrkTrackState *)
Definition: TrkTrackState.cxx:148
Trk::PrepRawData
Definition: PrepRawData.h:62
Trk::TrkBaseNode::getResiduals
virtual int getResiduals(double[2])=0
Trk::PrepRawData::identify
Identifier identify() const
return the identifier
TrigInDetTrackFitter::m_correctClusterPos
bool m_correctClusterPos
Definition: TrigInDetTrackFitter.h:65
Trk::Track::perigeeParameters
const Perigee * perigeeParameters() const
return Perigee.
Definition: Tracking/TrkEvent/TrkTrack/src/Track.cxx:163
TrigInDetTrackFitter::m_fitErrorsLowPt
std::atomic< size_t > m_fitErrorsLowPt
Definition: TrigInDetTrackFitter.h:60
Trk::TrackStateOnSurface
represents the track state (measurement, material, fit parameters and quality) at a surface.
Definition: TrackStateOnSurface.h:71
TrigInDetTrackFitter::extrapolate
Trk::TrkTrackState * extrapolate(Trk::TrkTrackState *, Trk::TrkPlanarSurface *, Trk::TrkPlanarSurface *, MagField::AtlasFieldCache &) const
Definition: TrigInDetTrackFitter.cxx:133
Trk::d0
@ d0
Definition: ParamDefs.h:63
TrackInfo.h
RIO_OnTrack.h
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
Region::Undefined
@ Undefined
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
SG::CondHandleKey::initialize
StatusCode initialize(bool used=true)
dqt_zlumi_alleff_HIST.B
B
Definition: dqt_zlumi_alleff_HIST.py:110
TrkFilteringNodes.h
TrigInDetTrackFitter::TrigInDetTrackFitter
TrigInDetTrackFitter(const std::string &, const std::string &, const IInterface *)
Definition: TrigInDetTrackFitter.cxx:58
MuonCalib::Legendre::coeff
constexpr double coeff(unsigned int l, unsigned int k)
Calculates the n-th coefficient of the legendre polynomial series.
Definition: LegendrePoly.h:73
TrigL2HitResidual
Definition: TrigL2HitResidual.h:24
TrigInDetTrackFitter::m_trackMaker
ToolHandle< ITrigDkfTrackMakerTool > m_trackMaker
Definition: TrigInDetTrackFitter.h:66
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
Trk::TrkTrackState::setSmootherGain
void setSmootherGain(double A[5][5])
Definition: TrkTrackState.cxx:152
TrigInDetTrackFitter::m_fitErrorsUnresolved
std::atomic< size_t > m_fitErrorsUnresolved
Definition: TrigInDetTrackFitter.h:58
Prompt::Def::Pt
@ Pt
Definition: VarHolder.h:76
a
TList * a
Definition: liststreamerinfos.cxx:10
python.compareTCTs.ratio
ratio
Definition: compareTCTs.py:295
TrigInDetTrackFitter::getMagneticField
void getMagneticField(double[3], double *, MagField::AtlasFieldCache &) const
Definition: TrigInDetTrackFitter.cxx:104
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Trk::PlaneSurface
Definition: PlaneSurface.h:64
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
MagField::AtlasFieldCache
Local cache for magnetic field (based on MagFieldServices/AtlasFieldSvcTLS.h)
Definition: AtlasFieldCache.h:43
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:67
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
LArNewCalib_DelayDump_OFC_Cali.idx
idx
Definition: LArNewCalib_DelayDump_OFC_Cali.py:69
Trk::TrkTrackState::setTrackCovariance
void setTrackCovariance(double A[5][5])
Definition: TrkTrackState.cxx:142
Trk::TrackStateOnSurface::Scatterer
@ Scatterer
This represents a scattering point on the track, and so will contain TrackParameters and MaterialEffe...
Definition: TrackStateOnSurface.h:113
TrigInDetTrackFitter::m_doMultScatt
bool m_doMultScatt
Definition: TrigInDetTrackFitter.h:63
MagField::AtlasFieldCache::getField
void getField(const double *ATH_RESTRICT xyz, double *ATH_RESTRICT bxyz, double *ATH_RESTRICT deriv=nullptr)
get B field value at given position xyz[3] is in mm, bxyz[3] is in kT if deriv[9] is given,...
Definition: AtlasFieldCache.cxx:42
Trk::TrackInfo::setParticleHypothesis
void setParticleHypothesis(const ParticleHypothesis &hypothesis)
Method re-setting the ParticleHypothesis.
Trk::TrkTrackState::setTrackState
void setTrackState(const double A[5])
Definition: TrkTrackState.cxx:138
Trk::TrkBaseNode::getChi2
double getChi2() const
Definition: TrkBaseNode.cxx:69
Trk::TrkTrackState::getTrackState
double getTrackState(int i)
Definition: Tracking/TrkFitter/TrkDistributedKalmanFilter/TrkDistributedKalmanFilter/TrkTrackState.h:47
Trk::TrkBaseNode::getPrepRawData
virtual const PrepRawData * getPrepRawData()
Definition: TrkBaseNode.cxx:39
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:14
TrigInDetTrackFitter::getUnbiasedResiduals
StatusCode getUnbiasedResiduals(const Trk::Track &, std::vector< TrigL2HitResidual > &, const EventContext &) const
Definition: TrigInDetTrackFitter.cxx:891
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
AthAlgTool
Definition: AthAlgTool.h:26
TrigInDetTrackFitter::m_fitErrorsDivergence
std::atomic< size_t > m_fitErrorsDivergence
Definition: TrigInDetTrackFitter.h:59
TrigInDetTrackFitter::createTrackStateOnSurface
Trk::TrackStateOnSurface * createTrackStateOnSurface(Trk::TrkBaseNode *pN, const bool, const EventContext &ctx) const
Definition: TrigInDetTrackFitter.cxx:820
FitQuality.h
TrigInDetTrackFitter.h
Trk::Surface
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:75
python.compressB64.c
def c
Definition: compressB64.py:93
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
TrkPlanarSurface.h
Trk::phi0
@ phi0
Definition: ParamDefs.h:65
Trk::TrkBaseNode::getInverseResidualVariance
virtual int getInverseResidualVariance(double[2][2])=0
Trk::TrackStateOnSurface::Measurement
@ Measurement
This is a measurement, and will at least contain a Trk::MeasurementBase.
Definition: TrackStateOnSurface.h:101
Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes
@ NumberOfTrackStateOnSurfaceTypes
Definition: TrackStateOnSurface.h:152
TrackStateOnSurface.h
Trk::TrkTrackState::getTrackCovariance
double getTrackCovariance(int i, int j)
Definition: Tracking/TrkFitter/TrkDistributedKalmanFilter/TrkDistributedKalmanFilter/TrkTrackState.h:51
Trk::StraightLineSurface
Definition: StraightLineSurface.h:51
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
SCT_ID::is_barrel
bool is_barrel(const Identifier &id) const
Test for barrel - WARNING: id MUST be sct id, otherwise answer is not accurate. Use SiliconID for gen...
Definition: SCT_ID.h:721
Trk::PrepRawData::detectorElement
virtual const TrkDetElementBase * detectorElement() const =0
return the detector element corresponding to this PRD The pointer will be zero if the det el is not d...
Identifier
Definition: IdentifierFieldParser.cxx:14
Generated on Mon Dec 23 2024 21:21:20 for ATLAS Offline Software by   doxygen 1.8.18