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