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  // pTE->applyMaterialEffects();
534  if(m_doMultScatt)
535  pTE->applyMultipleScattering();
536 
537  pTE->setTrackState(Rf);//restore
538 
539  if(m_doBremm)
540  pTE->applyEnergyLoss(1);
541 
542  AmgSymMatrix(5) Gi;
543  for(i=0;i<5;i++) for(j=i;j<5;j++)
544  {
545  Gi.fillSymmetric(i, j, pTE->getTrackCovariance(i,j));
546  }
547  Gi = Gi.inverse();
548 
549  for(i=0;i<5;i++) for(j=0;j<5;j++)
550  {
551  A[i][j]=0.0;
552  for(m=0;m<5;m++) A[i][j]+=AG[m][i]*Gi(m,j);
553  }
554  pTE->setPreviousState(pTS);
555  pTE->setSmootherGain(A);
556 
557  return pTE;
558 }
559 
560 void TrigInDetTrackFitter::fit(const TrackCollection& inputTracks, TrackCollection& fittedTracks, const EventContext& ctx, const Trk::ParticleHypothesis& matEffects) const
561 {
562  TrackCollection tmp;
563  fit(inputTracks,fittedTracks,tmp,ctx,matEffects,false);
564 }
565 
566 void TrigInDetTrackFitter::fit(const TrackCollection& inputTracks, TrackCollection& fittedTracks, TrackCollection& fittedTrackswTP, const EventContext& ctx, const Trk::ParticleHypothesis& matEffects, const bool addTPtoTSoS) const
567 {
568  MagField::AtlasFieldCache fieldCache;
569 
570  SG::ReadCondHandle<AtlasFieldCacheCondObj> fieldCondObj{m_fieldCondObjInputKey, ctx};
571  if (!fieldCondObj.isValid()) {
572  ATH_MSG_ERROR("Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
573  return;
574  }
575 
576  fieldCondObj->getInitializedCache (fieldCache);
577  fittedTracks.reserve(inputTracks.size());
578  if( addTPtoTSoS ) fittedTrackswTP.reserve(inputTracks.size());
579  for(auto trIt = inputTracks.begin(); trIt != inputTracks.end(); ++trIt) {
580  Trk::Track* fittedTrack = nullptr;
581  Trk::Track* fittedTrackwTP = nullptr;
582  std::tie(fittedTrack,fittedTrackwTP) = fitTrack(**trIt, fieldCache, matEffects, addTPtoTSoS);
583  if (fittedTrack!=nullptr) {
584  fittedTracks.push_back(fittedTrack);
585  }
586  if (addTPtoTSoS && fittedTrackwTP!=nullptr) {
587  fittedTrackswTP.push_back(fittedTrackwTP);
588  }
589  }
590 }
591 
592 std::pair<Trk::Track*,Trk::Track*> TrigInDetTrackFitter::fitTrack(const Trk::Track& recoTrack, MagField::AtlasFieldCache& fieldCache, const Trk::ParticleHypothesis& matEffects, const bool addTPtoTSoS) const {
593 
594  const Trk::TrackParameters* trackPars = recoTrack.perigeeParameters();
595  if(trackPars==nullptr) {
596  ATH_MSG_WARNING("Fit Failed -- TrkTrack has no parameters");
597  return std::make_pair(nullptr,nullptr);
598  }
599 
600  // 1. Create initial track state:
601  double Rk[5];
602  Rk[0] = trackPars->parameters()[Trk::d0];
603  Rk[1] = trackPars->parameters()[Trk::z0];
604  Rk[2] = trackPars->parameters()[Trk::phi0];
605  if(Rk[2]>M_PI) Rk[2]-=2*M_PI;
606  if(Rk[2]<-M_PI) Rk[2]+=2*M_PI;
607  double trk_theta = trackPars->parameters()[Trk::theta];
608  Rk[3] = trk_theta;
609  double trk_qOverP = trackPars->parameters()[Trk::qOverP];
610  Rk[4] = 1000.0*trk_qOverP;//MeV->GeV
611  double trk_Pt = sin(trk_theta)/trk_qOverP;
612 
613  if(fabs(trk_Pt)<100.0)
614  {
615  ATH_MSG_DEBUG("Estimated Pt is too low "<<trk_Pt<<" - skipping fit");
616  return std::make_pair(nullptr,nullptr);
617  }
618 
619  // 2. Create filtering nodes
620 
621  std::vector<Trk::TrkBaseNode*> vpTrkNodes;
622  std::vector<Trk::TrkTrackState*> vpTrackStates;
623  vpTrackStates.reserve(vpTrkNodes.size() + 1);
624  bool trackResult = m_trackMaker->createDkfTrack(recoTrack,vpTrkNodes, m_DChi2);
625  int nHits=vpTrkNodes.size();
626  ATH_MSG_VERBOSE(nHits<<" filtering nodes created");
627 
628  if(!trackResult) return std::make_pair(nullptr,nullptr);
629 
630  // 3. Main algorithm: filter and smoother (Rauch-Tung-Striebel)
631  m_nTracksTotal++;
632  Trk::TrkTrackState* pTS = new Trk::TrkTrackState(Rk);
633  double Gk[5][5] = {{100.0, 0, 0, 0, 0},
634  {0, 100.0, 0, 0, 0},
635  {0, 0, 0.01, 0, 0},
636  {0, 0, 0, 0.01, 0},
637  {0, 0, 0, 0, 0.1}};
638  pTS->setTrackCovariance(Gk);
639  if(m_doMultScatt)
640  pTS->setScatteringMode(1);
641  if(m_doBremm)
642  pTS->setScatteringMode(2);
643  vpTrackStates.push_back(pTS);
644 
645  //ATH_MSG_DEBUG("Initial chi2: "<<recoTrack.chi2()<<" track authorId: "<<recoTrack.algorithmId());
646  ATH_MSG_VERBOSE("Initial params: locT="<<Rk[0]<<" locL="<<Rk[1]<<" phi="<<Rk[2]
647  <<" theta="<<Rk[3]<<" Q="<<Rk[4]<<" pT="<<sin(Rk[3])/Rk[4]<<" GeV");
648 
649  bool OK=true;
650 
651  double chi2tot=0.0;
652  int ndoftot=-5;
653 
654  Trk::TrkPlanarSurface* pSB=nullptr;
655  Trk::TrkPlanarSurface* pSE=nullptr;
656  for(auto pnIt = vpTrkNodes.begin(); pnIt!=vpTrkNodes.end(); ++pnIt) {
657  pSE=(*pnIt)->getSurface();
658  Trk::TrkTrackState* pNS=extrapolate(pTS,pSB,pSE,fieldCache);
659 
660  pSB=pSE;
661  if(pNS!=nullptr) {
662  vpTrackStates.push_back(pNS);
663 
664  (*pnIt)->validateMeasurement(pNS);
665  ATH_MSG_VERBOSE("dChi2="<<(*pnIt)->getChi2());
666  if((*pnIt)->isValidated())
667  {
668  chi2tot+=(*pnIt)->getChi2();
669  ndoftot+=(*pnIt)->getNdof();
670  }
671  (*pnIt)->updateTrackState(pNS);
672  pTS=pNS;
673  double est_Pt = 1000.0*sin(pTS->getTrackState(3))/pTS->getTrackState(4);
674  if(fabs(est_Pt)<200.0)
675  {
676  ATH_MSG_VERBOSE("Estimated Pt is too low "<<est_Pt<<" - skipping fit");
677  m_fitErrorsLowPt++;
678  OK=false;break;
679  }
680  }
681  else
682  {
683  m_fitErrorsDivergence++;
684  OK=false;break;
685  }
686  }
687  Trk::Track* fittedTrack = nullptr;
688  Trk::Track* fittedTrackwTP = nullptr;
689  if(OK)
690  {
691  for(auto ptsIt = vpTrackStates.rbegin();ptsIt!=vpTrackStates.rend();++ptsIt)
692  {
693  (*ptsIt)->runSmoother();
694  }
695  pTS=(*vpTrackStates.begin());
696  //correct GeV->MeV
697 
698  correctScale(pTS);
699 
700  double qOverP = pTS->getTrackState(4);
701  double pt=sin(pTS->getTrackState(3))/pTS->getTrackState(4);
702  double phi0 = pTS->getTrackState(2);
703  if(phi0>M_PI) phi0-=2*M_PI;
704  if(phi0<-M_PI) phi0+=2*M_PI;
705  double theta = pTS->getTrackState(3);
706  double z0 = pTS->getTrackState(1);
707  double d0 = pTS->getTrackState(0);
708  bool bad_cov = false;
709  auto cov = AmgSymMatrix(5){};
710  for(int i=0;i<5;i++) {
711  double cov_diag = pTS->getTrackCovariance(i,i);
712  if (cov_diag < 0) {
713  bad_cov = true;//Diagonal elements must be positive
714  break;
715  ATH_MSG_DEBUG("REGTEST: cov(" << i << "," << i << ") =" << cov_diag << " < 0, reject track");
716  }
717  (cov)(i, i) = pTS->getTrackCovariance(i,i);
718  for(int j=i+1;j<5;j++) {
719  cov.fillSymmetric(i, j, pTS->getTrackCovariance(i,j));
720  }
721  }
722 
723  if((ndoftot<0) || (fabs(pt)<100.0) || (std::isnan(pt)) || bad_cov)
724  {
725  ATH_MSG_DEBUG("Fit failed - possibly floating point problem");
726  }
727  else
728  {
729  Trk::PerigeeSurface perigeeSurface;
730  auto perigee = std::make_unique<Trk::Perigee>(d0, z0, phi0, theta, qOverP, perigeeSurface, cov);
731  ATH_MSG_VERBOSE("perigee: " << *perigee);
732  std::unique_ptr<Trk::Perigee> perigeewTP = (addTPtoTSoS) ? std::make_unique<Trk::Perigee>(d0, z0, phi0, theta, qOverP, perigeeSurface, cov) : nullptr;
733 
734  std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern;
735  typePattern.set(Trk::TrackStateOnSurface::Perigee);
736  auto pParVec = std::make_unique<Trk::TrackStates>();
737  auto pParVecwTP = std::make_unique<Trk::TrackStates>();
738  if (m_correctClusterPos) {
739  pParVec->reserve(vpTrkNodes.size()+1);
740  pParVec->push_back(new Trk::TrackStateOnSurface(nullptr, std::move(perigee),nullptr, typePattern));
741  if (addTPtoTSoS) {
742  std::bitset<
743  Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes>
744  typePatternwTP;
745  typePatternwTP.set(Trk::TrackStateOnSurface::Perigee);
746  pParVecwTP->reserve(vpTrkNodes.size() + 1);
747  pParVecwTP->push_back(
748  new Trk::TrackStateOnSurface(nullptr, std::move(perigeewTP), nullptr, typePatternwTP));
749  }
750  for (auto pnIt = vpTrkNodes.begin(); pnIt != vpTrkNodes.end(); ++pnIt) {
751  if((*pnIt)->isValidated()) {
752  Trk::TrackStateOnSurface* pTSS = createTrackStateOnSurface(*pnIt,false);
753  Trk::TrackStateOnSurface* pTSSwTP = nullptr;
754  if( addTPtoTSoS ) pTSSwTP = createTrackStateOnSurface(*pnIt,true);
755  if(pTSS!=nullptr) {
756  pParVec->push_back(pTSS);
757  }
758  if(pTSSwTP!=nullptr) {
759  pParVecwTP->push_back(pTSSwTP);
760  }
761  }
762  }
763  }
764  else {
765  pParVec->reserve(recoTrack.trackStateOnSurfaces()->size());
766  pParVec->push_back(new Trk::TrackStateOnSurface(nullptr, std::move(perigee),nullptr, typePattern));
767 
768  for (auto tSOS = recoTrack.trackStateOnSurfaces()->begin(); tSOS != recoTrack.trackStateOnSurfaces()->end(); ++tSOS) {
769  //Don't store perigee - new perigee created above
770  if ((*tSOS)->type(Trk::TrackStateOnSurface::Perigee) == false) {
771  pParVec->push_back((*tSOS)->clone());
772  }
773  }
774  }
775  ATH_MSG_VERBOSE("Total chi2 ="<<chi2tot<<" NDOF="<<ndoftot);
776  if(msgLvl(MSG::VERBOSE)) {
777  double eta = -log(tan(0.5*theta));
778  ATH_MSG_VERBOSE("Fitted parameters: d0="<<d0<<" phi0="<<phi0<<" z0="<<z0
779  <<" eta0="<<eta<<" pt="<<pt);
780  }
781  auto pFQ= std::make_unique<Trk::FitQuality>(chi2tot,ndoftot);
782  Trk::TrackInfo info(recoTrack.info());
783  info.setParticleHypothesis(matEffects);
784  fittedTrack = new Trk::Track(info, std::move(pParVec), std::move(pFQ));//fittedTrack now owns pParVec and pFQ
785  if( addTPtoTSoS ) {
786  auto pFQwTP=std::make_unique<Trk::FitQuality>(chi2tot,ndoftot);
787  Trk::TrackInfo infowTP(recoTrack.info());
788  infowTP.setParticleHypothesis(matEffects);
789  fittedTrackwTP = new Trk::Track(infowTP, std::move(pParVecwTP), std::move(pFQwTP));//fittedTrack now owns pParVecwTP and pFQwTP
790  }
791  }
792  }
793  else
794  {
795  ATH_MSG_DEBUG("Forward Kalman filter: extrapolation failure ");
796  }
797 
798  for(auto pnIt = vpTrkNodes.begin(); pnIt!=vpTrkNodes.end(); ++pnIt) {
799  delete((*pnIt)->getSurface());
800  delete (*pnIt);
801  }
802  vpTrkNodes.clear();
803  for(auto ptsIt = vpTrackStates.begin();ptsIt!=vpTrackStates.end();++ptsIt) {
804  delete (*ptsIt);
805  }
806  vpTrackStates.clear();
807 
808  return std::make_pair(fittedTrack,fittedTrackwTP);
809 }
810 
811 Trk::TrackStateOnSurface* TrigInDetTrackFitter::createTrackStateOnSurface(Trk::TrkBaseNode* pN, const bool addTPtoTSoS) const
812 {
813  Trk::TrackStateOnSurface* pTSS=nullptr;
814  char type=pN->getNodeType();
815  std::unique_ptr<Trk::TrackParameters> pTP{};
816  if(type==0) return pTSS;
817 
818 
819  Trk::TrkTrackState* pTS=pN->getTrackState();
820  auto pM = AmgSymMatrix(5){};
821  for(int i=0;i<5;i++) {
822  for(int j=0;j<5;j++) {
823  (pM)(i,j)=pTS->getTrackCovariance(i,j);
824  }
825  }
826  const Trk::PrepRawData* pPRD=pN->getPrepRawData();
827 
828  if((type==1)||(type==2))
829  {
830  const Trk::Surface& rS = pPRD->detectorElement()->surface();
831  const Trk::PlaneSurface* pPS = dynamic_cast<const Trk::PlaneSurface*>(&rS);
832  if(pPS==nullptr) return pTSS;
833 
834  pTP=std::make_unique<Trk::AtaPlane>(pTS->getTrackState(0),
835  pTS->getTrackState(1),
836  pTS->getTrackState(2),
837  pTS->getTrackState(3),
838  pTS->getTrackState(4),*pPS,
839  std::move(pM));
840  }
841  else if(type==3)
842  {
843  const Trk::Surface& rS = pPRD->detectorElement()->surface(pPRD->identify());
844  const Trk::StraightLineSurface* pLS=dynamic_cast<const Trk::StraightLineSurface*>(&rS);
845  if(pLS==nullptr) return pTSS;
846 
847  if((pTS->getTrackState(2)<-M_PI) ||(pTS->getTrackState(2)>M_PI)) {
848  ATH_MSG_WARNING("Phi out of range when correcting Trk::TrackStateOnSurface");
849  }
850 
851 
852  pTP=std::make_unique<Trk::AtaStraightLine>(pTS->getTrackState(0),
853  pTS->getTrackState(1),
854  pTS->getTrackState(2),
855  pTS->getTrackState(3),
856  pTS->getTrackState(4),
857  *pLS,
858  std::move(pM));
859  }
860  if(pTP==nullptr) return nullptr;
861  std::unique_ptr<Trk::RIO_OnTrack> pRIO{m_ROTcreator->correct(*pPRD,*pTP)};
862  if(pRIO==nullptr) {
863  return nullptr;
864  }
865  std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern;
866  typePattern.set(Trk::TrackStateOnSurface::Measurement);
867  typePattern.set(Trk::TrackStateOnSurface::Scatterer);
868  auto pFQ=Trk::FitQualityOnSurface(pN->getChi2(),pN->getNdof());
869  if( addTPtoTSoS ) {
870  std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePatternwTP;
871  typePatternwTP.set(Trk::TrackStateOnSurface::Measurement);
872  typePatternwTP.set(Trk::TrackStateOnSurface::Scatterer);
873  auto pFQwTP=Trk::FitQualityOnSurface(pN->getChi2(),pN->getNdof());
874  pTSS = new Trk::TrackStateOnSurface(pFQwTP, pRIO->uniqueClone(), std::move(pTP), nullptr, typePatternwTP);
875  }
876  else {
877  pTSS = new Trk::TrackStateOnSurface(pFQ, std::move(pRIO), nullptr, nullptr, typePattern);
878  }
879  return pTSS;
880 }
881 
882 StatusCode TrigInDetTrackFitter::getUnbiasedResiduals(const Trk::Track& pT,
883  std::vector<TrigL2HitResidual>& vResid, const EventContext& ctx) const {
884 
885  const Trk::TrackParameters* trackPars = pT.perigeeParameters();
886  if(trackPars==nullptr) {
887  ATH_MSG_WARNING("Fit Failed -- TrkTrack has no parameters");
888  return StatusCode::FAILURE;
889  }
890 
891  MagField::AtlasFieldCache fieldCache;
892 
893  SG::ReadCondHandle<AtlasFieldCacheCondObj> fieldCondObj{m_fieldCondObjInputKey, ctx};
894  if (!fieldCondObj.isValid()) {
895  ATH_MSG_ERROR("Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
896  return StatusCode::FAILURE;
897  }
898 
899  fieldCondObj->getInitializedCache (fieldCache);
900  std::vector<Trk::TrkBaseNode*> vpTrkNodes;
901  std::vector<Trk::TrkTrackState*> vpTrackStates;
902  vResid.clear();
903  double trk_theta = trackPars->parameters()[Trk::theta];
904  double trk_qOverP = trackPars->parameters()[Trk::qOverP];
905  double Pt = sin(trk_theta)/trk_qOverP;
906  if(fabs(Pt)<100.0)
907  {
908  ATH_MSG_DEBUG("TrigL2ResidualCalculator failed -- Estimated Pt is too low "<<Pt);
909  return StatusCode::FAILURE;
910  }
911 
912  // 1. Create filtering nodes
913 
914  bool trackResult = m_trackMaker->createDkfTrack(pT,vpTrkNodes, m_DChi2);
915  if(!trackResult)
916  {
917  ATH_MSG_DEBUG("TrigDkfTrackMaker failed");
918  return StatusCode::FAILURE;
919  }
920 
921  bool OK=true;
922  std::vector<Trk::TrkBaseNode*>::iterator pnIt,pnEnd(vpTrkNodes.end());
923 
924  for(std::vector<Trk::TrkBaseNode*>::iterator pNodeIt=vpTrkNodes.begin();pNodeIt!=vpTrkNodes.end();
925  ++pNodeIt)
926  {
927  Trk::TrkBaseNode* pMaskedNode=(*pNodeIt);
928  Trk::TrkTrackState* pMaskedState=nullptr;
929 
930  // 2. Create initial track state:
931 
932  double Rk[5];
933  Rk[0] = trackPars->parameters()[Trk::d0];
934  Rk[1] = trackPars->parameters()[Trk::z0];
935  Rk[2] = trackPars->parameters()[Trk::phi0];
936  if(Rk[2]>M_PI) Rk[2]-=2*M_PI;
937  if(Rk[2]<-M_PI) Rk[2]+=2*M_PI;
938  trk_theta = trackPars->parameters()[Trk::theta];
939  Rk[3] = trk_theta;
940  Rk[4] = 1000.0*trackPars->parameters()[Trk::qOverP];//MeV->GeV
941  //No need to correct scale back - not returning track
942 
943  // 3. Main algorithm: filter and smoother (Rauch-Tung-Striebel)
944 
945  Trk::TrkTrackState* pTS = new Trk::TrkTrackState(Rk);
946  double Gk[5][5] = {{100.0, 0, 0, 0, 0},
947  {0, 100.0, 0, 0, 0},
948  {0, 0, 0.01, 0, 0},
949  {0, 0, 0, 0.01, 0},
950  {0, 0, 0, 0, 0.1}};
951  pTS->setTrackCovariance(Gk);
952  if(m_doMultScatt)
953  pTS->setScatteringMode(1);
954  if(m_doBremm)
955  pTS->setScatteringMode(2);
956  vpTrackStates.push_back(pTS);
957 
958  ATH_MSG_DEBUG("Initial params: locT="<<Rk[0]<<" locL="<<Rk[1]<<" phi="<<Rk[2]
959  <<" theta="<<Rk[3]<<" Q="<<Rk[4]<<" pT="<<sin(Rk[3])/Rk[4]);
960 
961  OK=true;
962  Trk::TrkPlanarSurface *pSB=nullptr,*pSE;
963 
964  for(pnIt=vpTrkNodes.begin();pnIt!=pnEnd;++pnIt)
965  {
966  pSE=(*pnIt)->getSurface();
967  Trk::TrkTrackState* pNS=extrapolate(pTS,pSB,pSE,fieldCache);
968 
969  pSB=pSE;
970  if(pNS!=nullptr)
971  {
972  vpTrackStates.push_back(pNS);
973 
974  (*pnIt)->validateMeasurement(pNS);
975  ATH_MSG_DEBUG("dChi2="<<(*pnIt)->getChi2());
976  if((*pnIt)!=pMaskedNode)
977  {
978  (*pnIt)->updateTrackState(pNS);
979  }
980  else
981  {
982  pMaskedState=pNS;
983  }
984  pTS=pNS;
985  Pt=sin(pTS->getTrackState(3))/pTS->getTrackState(4);
986  if(fabs(Pt)<0.2)
987  {
988  ATH_MSG_DEBUG("Estimated Pt is too low "<<Pt<<" - skipping fit");
989  OK=false;break;
990  }
991  }
992  else
993  {
994  OK=false;break;
995  }
996  }
997  if(OK)
998  {
999  std::vector<Trk::TrkTrackState*>::reverse_iterator ptsrIt(vpTrackStates.rbegin()),
1000  ptsrEnd(vpTrackStates.rend());
1001 
1002  for(;ptsrIt!=ptsrEnd;++ptsrIt)
1003  {
1004  (*ptsrIt)->runSmoother();
1005  }
1006 
1007  pMaskedNode->validateMeasurement(pMaskedState);
1008 
1009  double r[2],V[2][2];
1010 
1011  int nSize=pMaskedNode->getResiduals(r);
1012  nSize=pMaskedNode->getInverseResidualVariance(V);
1013  const Trk::PrepRawData* pPRD = pMaskedNode->getPrepRawData();
1014 
1015  Identifier id = pPRD->identify();
1016 
1017  Region region = Region::Undefined;
1018  if(m_idHelper->is_pixel(id))
1019  {
1020  region=(m_pixelId->is_barrel(id)) ? Region::PixBarrel: Region::PixEndcap;
1021  if (m_pixelId->is_blayer(id)) {
1022  region = Region::IBL;
1023  }
1024  }
1025  if(m_idHelper->is_sct(id))
1026  {
1027  region=(m_sctId->is_barrel(id)) ? Region::SctBarrel : Region::SctEndcap;
1028  }
1029  if(nSize==1) {
1030  if(V[0][0]>0.0) {
1031  vResid.push_back(TrigL2HitResidual(id,region,r[0],r[0]*sqrt(V[0][0])));
1032  }
1033  else {
1034  OK=false;
1035  break;
1036  }
1037  }
1038  else {
1039  if((V[0][0]>0.0) && (V[1][1]>0.0)) {
1040  vResid.push_back(TrigL2HitResidual(id,region,r[0],r[0]*sqrt(V[0][0]),
1041  r[1],r[1]*sqrt(V[1][1])));
1042  }
1043  else {
1044  OK=false;
1045  break;
1046  }
1047  }
1048  }
1049  else
1050  {
1051  ATH_MSG_DEBUG("Forward Kalman filter: extrapolation failure ");
1052  vResid.clear();
1053  }
1054  for(std::vector<Trk::TrkTrackState*>::iterator ptsIt=vpTrackStates.begin();
1055  ptsIt!=vpTrackStates.end();++ptsIt) delete (*ptsIt);
1056  vpTrackStates.clear();
1057  if(!OK) break;
1058  }
1059  pnIt=vpTrkNodes.begin();pnEnd=vpTrkNodes.end();
1060  for(;pnIt!=pnEnd;++pnIt)
1061  {
1062  delete((*pnIt)->getSurface());
1063  delete (*pnIt);
1064  }
1065  vpTrkNodes.clear();
1066 
1067  if(OK) return StatusCode::SUCCESS;
1068  else return StatusCode::FAILURE;
1069 
1070 }
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
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
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:257
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
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
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:79
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.
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:71
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:560
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:70
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
TrigInDetTrackFitter::createTrackStateOnSurface
Trk::TrackStateOnSurface * createTrackStateOnSurface(Trk::TrkBaseNode *pN, const bool) const
Definition: TrigInDetTrackFitter.cxx:811
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:52
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
dqt_zlumi_alleff_HIST.A
A
Definition: dqt_zlumi_alleff_HIST.py:110
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
TrkTrackState.h
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
python.TrigInDetConfig.inputTracks
inputTracks
Definition: TrigInDetConfig.py:168
lumiFormat.i
int i
Definition: lumiFormat.py:92
Identifier
Definition: DetectorDescription/Identifier/Identifier/Identifier.h:32
beamspotman.n
n
Definition: beamspotman.py:731
TrkBaseNode.h
Trk::theta
@ theta
Definition: ParamDefs.h:72
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
TrigInDetTrackFitter::fitTrack
std::pair< Trk::Track *, Trk::Track * > fitTrack(const Trk::Track &, MagField::AtlasFieldCache &, const Trk::ParticleHypothesis &matEffects=Trk::pion, const bool addTPtoTSoS=false) const
Definition: TrigInDetTrackFitter.cxx:592
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:69
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
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
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
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:73
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
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:882
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
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:71
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...
Generated on Sun Jun 30 2024 21:30:39 for ATLAS Offline Software by   doxygen 1.8.18