ATLAS Offline Software
Loading...
Searching...
No Matches
TrigInDetTrackFitter Class Reference

#include <TrigInDetTrackFitter.h>

Inheritance diagram for TrigInDetTrackFitter:
Collaboration diagram for TrigInDetTrackFitter:

Public Member Functions

 TrigInDetTrackFitter (const std::string &, const std::string &, const IInterface *)
virtual StatusCode initialize ()
virtual StatusCode finalize ()
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
void fit (const TrackCollection &, TrackCollection &, const EventContext &, const Trk::ParticleHypothesis &matEffects=Trk::pion) const
void fit (const TrackCollection &, TrackCollection &, TrackCollection &, const EventContext &, const Trk::ParticleHypothesis &matEffects=Trk::pion, const bool addTPtoTSoS=false) const
StatusCode getUnbiasedResiduals (const Trk::Track &, std::vector< TrigL2HitResidual > &, const EventContext &) const
ServiceHandle< StoreGateSvc > & evtStore ()
 The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
const ServiceHandle< StoreGateSvc > & detStore () const
 The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.
virtual StatusCode sysInitialize () override
 Perform system initialization for an algorithm.
virtual StatusCode sysStart () override
 Handle START transition.
virtual std::vector< Gaudi::DataHandle * > inputHandles () const override
 Return this algorithm's input handles.
virtual std::vector< Gaudi::DataHandle * > outputHandles () const override
 Return this algorithm's output handles.
Gaudi::Details::PropertyBase & declareProperty (Gaudi::Property< T, V, H > &t)
void updateVHKA (Gaudi::Details::PropertyBase &)
MsgStream & msg () const
bool msgLvl (const MSG::Level lvl) const

Static Public Member Functions

static const InterfaceID & interfaceID ()

Protected Member Functions

void renounceArray (SG::VarHandleKeyArray &handlesArray)
 remove all handles from I/O resolution
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > renounce (T &h)
void extraDeps_update_handler (Gaudi::Details::PropertyBase &ExtraDeps)
 Add StoreName to extra input/output deps as needed.

Private Types

typedef ServiceHandle< StoreGateSvcStoreGateSvc_t

Private Member Functions

Trk::TrkTrackStateextrapolate (Trk::TrkTrackState *, Trk::TrkPlanarSurface *, Trk::TrkPlanarSurface *, MagField::AtlasFieldCache &) const
void getMagneticField (double[3], double *, MagField::AtlasFieldCache &) const
void correctScale (Trk::TrkTrackState *) const
Trk::TrackStateOnSurfacecreateTrackStateOnSurface (Trk::TrkBaseNode *pN, const bool, const EventContext &ctx) const
Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
 specialization for handling Gaudi::Property<SG::VarHandleKey>

Private Attributes

std::atomic< size_t > m_nTracksTotal
std::atomic< size_t > m_fitErrorsUnresolved
std::atomic< size_t > m_fitErrorsDivergence
std::atomic< size_t > m_fitErrorsLowPt
double m_DChi2
bool m_doMultScatt
bool m_doBremm
bool m_correctClusterPos
ToolHandle< ITrigDkfTrackMakerToolm_trackMaker
ToolHandle< Trk::IRIO_OnTrackCreatorm_ROTcreator
SG::ReadCondHandleKey< AtlasFieldCacheCondObjm_fieldCondObjInputKey {this, "AtlasFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"}
const PixelIDm_pixelId = nullptr
const SCT_IDm_sctId = nullptr
const AtlasDetectorIDm_idHelper = nullptr
StoreGateSvc_t m_evtStore
 Pointer to StoreGate (event store by default)
StoreGateSvc_t m_detStore
 Pointer to StoreGate (detector store by default)
std::vector< SG::VarHandleKeyArray * > m_vhka
bool m_varHandleArraysDeclared

Detailed Description

Definition at line 33 of file TrigInDetTrackFitter.h.

Member Typedef Documentation

◆ StoreGateSvc_t

typedef ServiceHandle<StoreGateSvc> AthCommonDataStore< AthCommonMsg< AlgTool > >::StoreGateSvc_t
privateinherited

Definition at line 388 of file AthCommonDataStore.h.

Constructor & Destructor Documentation

◆ TrigInDetTrackFitter()

TrigInDetTrackFitter::TrigInDetTrackFitter ( const std::string & t,
const std::string & n,
const IInterface * p )

Definition at line 58 of file TrigInDetTrackFitter.cxx.

60 : 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" );
74}
AthAlgTool()
Default constructor:
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
std::atomic< size_t > m_nTracksTotal
ToolHandle< Trk::IRIO_OnTrackCreator > m_ROTcreator
std::atomic< size_t > m_fitErrorsUnresolved
ToolHandle< ITrigDkfTrackMakerTool > m_trackMaker
std::atomic< size_t > m_fitErrorsLowPt
std::atomic< size_t > m_fitErrorsDivergence

Member Function Documentation

◆ correctScale()

void TrigInDetTrackFitter::correctScale ( Trk::TrkTrackState * pTS) const
private

Definition at line 111 of file TrigInDetTrackFitter.cxx.

111 {
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}
void setTrackCovariance(double A[5][5])
void setTrackState(const double A[5])

◆ createTrackStateOnSurface()

Trk::TrackStateOnSurface * TrigInDetTrackFitter::createTrackStateOnSurface ( Trk::TrkBaseNode * pN,
const bool addTPtoTSoS,
const EventContext & ctx ) const
private

Definition at line 820 of file TrigInDetTrackFitter.cxx.

821{
822 Trk::TrackStateOnSurface* pTSS=nullptr;
823 char type=pN->getNodeType();
824 std::unique_ptr<Trk::TrackParameters> pTP{};
825 if(type==0) return pTSS;
826
827
828 Trk::TrkTrackState* pTS=pN->getTrackState();
829 auto pM = AmgSymMatrix(5){};
830 for(int i=0;i<5;i++) {
831 for(int j=0;j<5;j++) {
832 (pM)(i,j)=pTS->getTrackCovariance(i,j);
833 }
834 }
835 const Trk::PrepRawData* pPRD=pN->getPrepRawData();
836
837 if((type==1)||(type==2))
838 {
839 const Trk::Surface& rS = pPRD->detectorElement()->surface();
840 const Trk::PlaneSurface* pPS = dynamic_cast<const Trk::PlaneSurface*>(&rS);
841 if(pPS==nullptr) return pTSS;
842
843 pTP=std::make_unique<Trk::AtaPlane>(pTS->getTrackState(0),
844 pTS->getTrackState(1),
845 pTS->getTrackState(2),
846 pTS->getTrackState(3),
847 pTS->getTrackState(4),*pPS,
848 std::move(pM));
849 }
850 else if(type==3)
851 {
852 const Trk::Surface& rS = pPRD->detectorElement()->surface(pPRD->identify());
853 const Trk::StraightLineSurface* pLS=dynamic_cast<const Trk::StraightLineSurface*>(&rS);
854 if(pLS==nullptr) return pTSS;
855
856 if((pTS->getTrackState(2)<-M_PI) ||(pTS->getTrackState(2)>M_PI)) {
857 ATH_MSG_WARNING("Phi out of range when correcting Trk::TrackStateOnSurface");
858 }
859
860
861 pTP=std::make_unique<Trk::AtaStraightLine>(pTS->getTrackState(0),
862 pTS->getTrackState(1),
863 pTS->getTrackState(2),
864 pTS->getTrackState(3),
865 pTS->getTrackState(4),
866 *pLS,
867 std::move(pM));
868 }
869 if(pTP==nullptr) return nullptr;
870 std::unique_ptr<Trk::RIO_OnTrack> pRIO{m_ROTcreator->correct(*pPRD,*pTP,ctx)};
871 if(pRIO==nullptr) {
872 return nullptr;
873 }
874 std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern;
877 auto pFQ=Trk::FitQualityOnSurface(pN->getChi2(),pN->getNdof());
878 if( addTPtoTSoS ) {
879 std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePatternwTP;
880 typePatternwTP.set(Trk::TrackStateOnSurface::Measurement);
881 typePatternwTP.set(Trk::TrackStateOnSurface::Scatterer);
882 auto pFQwTP=Trk::FitQualityOnSurface(pN->getChi2(),pN->getNdof());
883 pTSS = new Trk::TrackStateOnSurface(pFQwTP, pRIO->uniqueClone(), std::move(pTP), nullptr, typePatternwTP);
884 }
885 else {
886 pTSS = new Trk::TrackStateOnSurface(pFQ, std::move(pRIO), nullptr, nullptr, typePattern);
887 }
888 return pTSS;
889}
#define M_PI
#define ATH_MSG_WARNING(x)
#define AmgSymMatrix(dim)
virtual const TrkDetElementBase * detectorElement() const =0
return the detector element corresponding to this PRD The pointer will be zero if the det el is not d...
Identifier identify() const
return the identifier
@ Measurement
This is a measurement, and will at least contain a Trk::MeasurementBase.
@ Scatterer
This represents a scattering point on the track, and so will contain TrackParameters and MaterialEffe...
double getChi2() const
TrkTrackState * getTrackState()
virtual const PrepRawData * getPrepRawData()
int getNdof() const
virtual char getNodeType()
virtual const Surface & surface() const =0
Return surface associated with this detector element.

◆ declareGaudiProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty ( Gaudi::Property< T, V, H > & hndl,
const SG::VarHandleKeyType &  )
inlineprivateinherited

specialization for handling Gaudi::Property<SG::VarHandleKey>

Definition at line 156 of file AthCommonDataStore.h.

158 {
160 hndl.value(),
161 hndl.documentation());
162
163 }

◆ declareProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( Gaudi::Property< T, V, H > & t)
inlineinherited

Definition at line 145 of file AthCommonDataStore.h.

145 {
146 typedef typename SG::HandleClassifier<T>::type htype;
148 }
Gaudi::Details::PropertyBase & declareGaudiProperty(Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
specialization for handling Gaudi::Property<SG::VarHandleKey>

◆ detStore()

const ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< AlgTool > >::detStore ( ) const
inlineinherited

The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 95 of file AthCommonDataStore.h.

◆ evtStore()

ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< AlgTool > >::evtStore ( )
inlineinherited

The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 85 of file AthCommonDataStore.h.

◆ extraDeps_update_handler()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::extraDeps_update_handler ( Gaudi::Details::PropertyBase & ExtraDeps)
protectedinherited

Add StoreName to extra input/output deps as needed.

use the logic of the VarHandleKey to parse the DataObjID keys supplied via the ExtraInputs and ExtraOuputs Properties to add the StoreName if it's not explicitly given

◆ extrapolate()

Trk::TrkTrackState * TrigInDetTrackFitter::extrapolate ( Trk::TrkTrackState * pTS,
Trk::TrkPlanarSurface * pSB,
Trk::TrkPlanarSurface * pSE,
MagField::AtlasFieldCache & fieldCache ) const
private

Definition at line 133 of file TrigInDetTrackFitter.cxx.

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;
153 }
154 }
155
156 if(!samePlane) {
157
158 double gP[3],gPi[3],lP[3],gV[3],a,b,c,s,J0[7][5],descr,CQ,Ac,Av,Cc;
159 double V[3],P[3],M[3][3],D[4],Jm[7][7],
160 J1[5][7],gB[3],gBi[3],gBf[3],dBds[3],Buf[5][7],DVx,DVy,DVz;
161 int nStep,nStepMax;
162 double sl,ds;
163
164 double sint,cost,sinf,cosf;
165 sint=sin(pTS->getTrackState(3));cosf=cos(pTS->getTrackState(2));
166 sinf=sin(pTS->getTrackState(2));cost=cos(pTS->getTrackState(3));
167 gV[0]=sint*cosf;gV[1]=sint*sinf;gV[2]=cost;CQ=C*pTS->getTrackState(4);
168
169 memset(&J0[0][0],0,sizeof(J0));
170
171 if(pSB!=nullptr)
172 {
173 double L[3][3];
174 lP[0]=pTS->getTrackState(0);lP[1]=pTS->getTrackState(1);lP[2]=0.0;
175 pSB->transformPointToGlobal(lP,gP);
176 for(i=0;i<3;i++) for(j=0;j<3;j++) L[i][j]=pSB->getInvRotMatrix(i,j);
177
178 J0[0][0]=L[0][0];J0[0][1]=L[0][1];
179 J0[1][0]=L[1][0];J0[1][1]=L[1][1];
180 J0[2][0]=L[2][0];J0[2][1]=L[2][1];
181 J0[3][2]=-sinf*sint;J0[3][3]=cosf*cost;
182 J0[4][2]= cosf*sint;J0[4][3]=sinf*cost;
183 J0[5][3]=-sint;
184 J0[6][4]=1.0;
185 }
186 else
187 {
188 gP[0]=-pTS->getTrackState(0)*sinf;
189 gP[1]= pTS->getTrackState(0)*cosf;
190 gP[2]= pTS->getTrackState(1);
191 J0[0][0]=-sinf;J0[0][2]=-pTS->getTrackState(0)*cosf;
192 J0[1][0]= cosf;J0[1][2]=-pTS->getTrackState(0)*sinf;
193 J0[2][1]=1.0;
194 J0[3][2]=-sinf*sint;J0[3][3]=cosf*cost;
195 J0[4][2]= cosf*sint;J0[4][3]=sinf*cost;
196 J0[5][3]=-sint;
197 J0[6][4]=1.0;
198 }
199 for(i=0;i<4;i++) D[i]=pSE->getPar(i);
200 for(i=0;i<3;i++) gPi[i]=gP[i];
201
202 getMagneticField(gP,gB, fieldCache);
203
204 for(i=0;i<3;i++) gBi[i]=gB[i];
205
206 c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
207 b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
208 a=0.5*CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
209 gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
210 gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
211
212 descr=b*b-4.0*a*c;
213
214 if(descr<0.0)
215 {
216 // printf("D<0 - extrapolation failed\n");
217 return nullptr;
218 }
219
220 bool useExpansion=true;
221 double ratio = 4*a*c/(b*b);
222
223 if(fabs(ratio)>0.1)
224 useExpansion = false;
225
226 if(useExpansion) {
227 sl=-c/b;
228 sl=sl*(1-a*sl/b);
229 }
230 else {
231 int signb = (b<0.0)?-1:1;
232 sl = (-b+signb*sqrt(descr))/(2*a);
233 }
234
235 if(fabs(sl)<minStep) nStepMax=1;
236 else
237 {
238 nStepMax=(int)(fabs(sl)/minStep)+1;
239 }
240 if((nStepMax<0)||(nStepMax>1000))
241 {
242 return nullptr;
243 }
244 Av=sl*CQ;
245 Ac=0.5*sl*Av;
246 DVx=gV[1]*gB[2]-gV[2]*gB[1];
247 DVy=gV[2]*gB[0]-gV[0]*gB[2];
248 DVz=gV[0]*gB[1]-gV[1]*gB[0];
249
250 P[0]=gP[0]+gV[0]*sl+Ac*DVx;
251 P[1]=gP[1]+gV[1]*sl+Ac*DVy;
252 P[2]=gP[2]+gV[2]*sl+Ac*DVz;
253 V[0]=gV[0]+Av*DVx;
254 V[1]=gV[1]+Av*DVy;
255 V[2]=gV[2]+Av*DVz;
256
257 getMagneticField(P,gB,fieldCache);
258
259 for(i=0;i<3;i++) gBf[i]=gB[i];
260 for(i=0;i<3;i++)
261 {
262 dBds[i]=(gBf[i]-gBi[i])/sl;
263 gB[i]=gBi[i];
264 }
265 nStep=nStepMax;
266 while(nStep>0)
267 {
268 c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
269 b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
270 a=0.5*CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
271 gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
272 gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
273
274 ratio = 4*a*c/(b*b);
275 if(fabs(ratio)>0.1)
276 useExpansion = false;
277 else useExpansion = true;
278
279 if(useExpansion) {
280 sl=-c/b;
281 sl=sl*(1-a*sl/b);
282 }
283 else {
284 descr=b*b-4.0*a*c;
285 if(descr<0.0)
286 {
287 // printf("D<0 - extrapolation failed\n");
288 return nullptr;
289 }
290 int signb = (b<0.0)?-1:1;
291 sl = (-b+signb*sqrt(descr))/(2*a);
292 }
293
294 ds=sl/nStep;
295 Av=ds*CQ;
296 Ac=0.5*ds*Av;
297 DVx=gV[1]*gB[2]-gV[2]*gB[1];
298 DVy=gV[2]*gB[0]-gV[0]*gB[2];
299 DVz=gV[0]*gB[1]-gV[1]*gB[0];
300
301 P[0]=gP[0]+gV[0]*ds+Ac*DVx;
302 P[1]=gP[1]+gV[1]*ds+Ac*DVy;
303 P[2]=gP[2]+gV[2]*ds+Ac*DVz;
304 V[0]=gV[0]+Av*DVx;
305 V[1]=gV[1]+Av*DVy;
306 V[2]=gV[2]+Av*DVz;
307 for(i=0;i<3;i++)
308 {
309 gV[i]=V[i];gP[i]=P[i];
310 }
311 for(i=0;i<3;i++) gB[i]+=dBds[i]*ds;
312 nStep--;
313 }
314 pSE->transformPointToLocal(gP,lP);
315 Rf[0]=lP[0];Rf[1]=lP[1];
316 Rf[2]=atan2(V[1],V[0]);
317
318 if(fabs(V[2])>1.0)
319 {
320 return nullptr;
321 }
322
323 Rf[3]=acos(V[2]);
324 Rf[4]=pTS->getTrackState(4);
325
326 gV[0]=sint*cosf;gV[1]=sint*sinf;gV[2]=cost;
327
328 for(i=0;i<4;i++) D[i]=pSE->getPar(i);
329 for(i=0;i<3;i++) gP[i]=gPi[i];
330
331 for(i=0;i<3;i++)
332 {
333 gB[i]=0.5*(gBi[i]+gBf[i]);
334 }
335
336 c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
337 b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
338 a=CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
339 gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
340 gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
341
342 ratio = 4*a*c/(b*b);
343 if(fabs(ratio)>0.1)
344 useExpansion = false;
345 else useExpansion = true;
346
347 if(useExpansion) {
348 s=-c/b;
349 s=s*(1-a*s/b);
350 }
351 else {
352 descr=b*b-4.0*a*c;
353 if(descr<0.0)
354 {
355 // printf("D<0 - extrapolation failed\n");
356 return nullptr;
357 }
358 int signb = (b<0.0)?-1:1;
359 s = (-b+signb*sqrt(descr))/(2*a);
360 }
361
362 Av=s*CQ;
363 Ac=0.5*s*Av;
364 Cc=0.5*s*s*C;
365
366 DVx=gV[1]*gB[2]-gV[2]*gB[1];
367 DVy=gV[2]*gB[0]-gV[0]*gB[2];
368 DVz=gV[0]*gB[1]-gV[1]*gB[0];
369
370 P[0]=gP[0]+gV[0]*s+Ac*DVx;
371 P[1]=gP[1]+gV[1]*s+Ac*DVy;
372 P[2]=gP[2]+gV[2]*s+Ac*DVz;
373
374 V[0]=gV[0]+Av*DVx;V[1]=gV[1]+Av*DVy;V[2]=gV[2]+Av*DVz;
375 if (std::abs(V[2]) > 1.0) {
376 return nullptr;
377 }
378
379 pSE->transformPointToLocal(P,lP);
380
381 memset(&Jm[0][0],0,sizeof(Jm));
382
383 for(i=0;i<3;i++) for(j=0;j<3;j++) M[i][j]=pSE->getRotMatrix(i,j);
384
385 double coeff[3], dadVx,dadVy,dadVz,dadQ,dsdx,dsdy,dsdz,dsdVx,dsdVy,dsdVz,dsdQ;
386 coeff[0]=-c*c/(b*b*b);
387 coeff[1]=c*(1.0+3.0*c*a/(b*b))/(b*b);
388 coeff[2]=-(1.0+2.0*c*a/(b*b))/b;
389
390 dadVx=0.5*CQ*(-D[1]*gB[2]+D[2]*gB[1]);
391 dadVy=0.5*CQ*( D[0]*gB[2]-D[2]*gB[0]);
392 dadVz=0.5*CQ*(-D[0]*gB[1]+D[1]*gB[0]);
393 dadQ=0.5*C*(D[0]*DVx+D[1]*DVy+D[2]*DVz);
394
395 dsdx=coeff[2]*D[0];
396 dsdy=coeff[2]*D[1];
397 dsdz=coeff[2]*D[2];
398 dsdVx=coeff[0]*dadVx+coeff[1]*D[0];
399 dsdVy=coeff[0]*dadVy+coeff[1]*D[1];
400 dsdVz=coeff[0]*dadVz+coeff[1]*D[2];
401 dsdQ=coeff[0]*dadQ;
402
403 Jm[0][0]=1.0+V[0]*dsdx;
404 Jm[0][1]= V[0]*dsdy;
405 Jm[0][2]= V[0]*dsdz;
406
407 Jm[0][3]= s+V[0]*dsdVx;
408 Jm[0][4]= V[0]*dsdVy+Ac*gB[2];
409 Jm[0][5]= V[0]*dsdVz-Ac*gB[1];
410 Jm[0][6]= V[0]*dsdQ+Cc*DVx;
411
412 Jm[1][0]= V[1]*dsdx;
413 Jm[1][1]=1.0+V[1]*dsdy;
414 Jm[1][2]= V[1]*dsdz;
415
416 Jm[1][3]= V[1]*dsdVx-Ac*gB[2];
417 Jm[1][4]= s+V[1]*dsdVy;
418 Jm[1][5]= V[1]*dsdVz+Ac*gB[0];
419 Jm[1][6]= V[1]*dsdQ+Cc*DVy;
420
421 Jm[2][0]= V[2]*dsdx;
422 Jm[2][1]= V[2]*dsdy;
423 Jm[2][2]=1.0+V[2]*dsdz;
424 Jm[2][3]= V[2]*dsdVx+Ac*gB[1];
425 Jm[2][4]= V[2]*dsdVy-Ac*gB[0];
426 Jm[2][5]= s+V[2]*dsdVz;
427 Jm[2][6]= V[2]*dsdQ+Cc*DVz;
428
429 Jm[3][0]=dsdx*CQ*DVx;
430 Jm[3][1]=dsdy*CQ*DVx;
431 Jm[3][2]=dsdz*CQ*DVx;
432
433 Jm[3][3]=1.0+dsdVx*CQ*DVx;
434 Jm[3][4]=CQ*(dsdVy*DVx+s*gB[2]);
435 Jm[3][5]=CQ*(dsdVz*DVx-s*gB[1]);
436
437 Jm[3][6]=(CQ*dsdQ+C*s)*DVx;
438
439 Jm[4][0]=dsdx*CQ*DVy;
440 Jm[4][1]=dsdy*CQ*DVy;
441 Jm[4][2]=dsdz*CQ*DVy;
442
443 Jm[4][3]=CQ*(dsdVx*DVy-s*gB[2]);
444 Jm[4][4]=1.0+dsdVy*CQ*DVy;
445 Jm[4][5]=CQ*(dsdVz*DVy+s*gB[0]);
446
447 Jm[4][6]=(CQ*dsdQ+C*s)*DVy;
448
449 Jm[5][0]=dsdx*CQ*DVz;
450 Jm[5][1]=dsdy*CQ*DVz;
451 Jm[5][2]=dsdz*CQ*DVz;
452 Jm[5][3]=CQ*(dsdVx*DVz+s*gB[1]);
453 Jm[5][4]=CQ*(dsdVy*DVz-s*gB[0]);
454 Jm[5][5]=1.0+dsdVz*CQ*DVz;
455 Jm[5][6]=(CQ*dsdQ+C*s)*DVz;
456
457 Jm[6][6]=1.0;
458
459 memset(&J1[0][0],0,sizeof(J1));
460
461 J1[0][0]=M[0][0];J1[0][1]=M[0][1];J1[0][2]=M[0][2];
462 J1[1][0]=M[1][0];J1[1][1]=M[1][1];J1[1][2]=M[1][2];
463 J1[2][3]=-V[1]/(V[0]*V[0]+V[1]*V[1]);
464 J1[2][4]= V[0]/(V[0]*V[0]+V[1]*V[1]);
465 J1[3][5]=-1.0/sqrt(1-V[2]*V[2]);
466 J1[4][6]=1.0;
467
468 for(i=0;i<7;i++)
469 {
470 for(j=0;j<2;j++)
471 Buf[j][i]=J1[j][0]*Jm[0][i]+J1[j][1]*Jm[1][i]+J1[j][2]*Jm[2][i];
472 Buf[2][i]=J1[2][3]*Jm[3][i]+J1[2][4]*Jm[4][i];
473 Buf[3][i]=J1[3][5]*Jm[5][i];
474 Buf[4][i]=Jm[6][i];
475 }
476
477 if(pSB!=nullptr)
478 {
479 for(i=0;i<5;i++)
480 {
481 J[i][0]=Buf[i][0]*J0[0][0]+Buf[i][1]*J0[1][0]+Buf[i][2]*J0[2][0];
482 J[i][1]=Buf[i][0]*J0[0][1]+Buf[i][1]*J0[1][1]+Buf[i][2]*J0[2][1];
483 J[i][2]=Buf[i][3]*J0[3][2]+Buf[i][4]*J0[4][2];
484 J[i][3]=Buf[i][3]*J0[3][3]+Buf[i][4]*J0[4][3]+Buf[i][5]*J0[5][3];
485 J[i][4]=Buf[i][6];
486 }
487 }
488 else
489 {
490 for(i=0;i<5;i++)
491 {
492 J[i][0]=Buf[i][0]*J0[0][0]+Buf[i][1]*J0[1][0];
493 J[i][1]=Buf[i][2];
494 J[i][2]=Buf[i][0]*J0[0][2]+Buf[i][1]*J0[1][2]+Buf[i][3]*J0[3][2]+Buf[i][4]*J0[4][2];
495 J[i][3]=Buf[i][3]*J0[3][3]+Buf[i][4]*J0[4][3]+Buf[i][5]*J0[5][3];
496 J[i][4]=Buf[i][6];
497 }
498 }
499 }
500 else {
501 Rf[0]=pTS->getTrackState(0);
502 Rf[1]=pTS->getTrackState(1);
503 Rf[2]=pTS->getTrackState(2);
504 Rf[3]=pTS->getTrackState(3);
505 Rf[4]=pTS->getTrackState(4);
506 memset(&J[0][0],0,sizeof(J));
507 for(i=0;i<5;i++) J[i][i]=1.0;
508 }
509
510 for(i=0;i<5;i++) for(j=0;j<5;j++)
511 {
512 AG[i][j]=0.0;for(m=0;m<5;m++) AG[i][j]+=J[i][m]*pTS->getTrackCovariance(m,j);
513 }
514 for(i=0;i<5;i++) for(j=i;j<5;j++)
515 {
516 Gf[i][j]=0.0;
517 for(m=0;m<5;m++) Gf[i][j]+=AG[i][m]*J[j][m];
518 Gf[j][i]=Gf[i][j];
519 }
520
521 Trk::TrkTrackState* pTE=new Trk::TrkTrackState(pTS);
522
523 //workaround to keep using existing TrkTrackState code
524 double Rtmp[5];
525
526 for(i=0;i<4;i++) Rtmp[i] = Rf[i];
527 Rtmp[4] = 0.001*Rf[4];//GeV->MeV
528
529 pTE->setTrackState(Rtmp);
530 pTE->setTrackCovariance(Gf);
531 pTE->attachToSurface(pSE);
532
533 if(m_doMultScatt)
535
536 pTE->setTrackState(Rf);//restore
537
538 if(m_doBremm)
539 pTE->applyEnergyLoss(1);
540
541 //quick check for covariance sanity
542
543 for(int idx=0;idx<5;idx++) {
544 if(pTE->getTrackCovariance(idx,idx) < 0) {
545 ATH_MSG_DEBUG("REGTEST: cov(" << idx << "," << idx << ") =" << pTE->getTrackCovariance(idx,idx) << " < 0, reject track");
546 delete pTE;
547 return nullptr;
548 }
549 }
550
551 AmgSymMatrix(5) Gi;
552 for(i=0;i<5;i++) for(j=i;j<5;j++)
553 {
554 Gi.fillSymmetric(i, j, pTE->getTrackCovariance(i,j));
555 }
556 Gi = Gi.inverse();
557
558 for(i=0;i<5;i++) for(j=0;j<5;j++)
559 {
560 A[i][j]=0.0;
561 for(m=0;m<5;m++) A[i][j]+=AG[m][i]*Gi(m,j);
562 }
563 pTE->setPreviousState(pTS);
564 pTE->setSmootherGain(A);
565
566 return pTE;
567}
#define ATH_MSG_DEBUG(x)
static Double_t a
static Double_t P(Double_t *tt, Double_t *par)
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
void getMagneticField(double[3], double *, MagField::AtlasFieldCache &) const
void transformPointToLocal(const double *, double *)
void transformPointToGlobal(const double *, double *)
double getInvRotMatrix(int, int)
double getRotMatrix(int, int)
void setPreviousState(TrkTrackState *)
void attachToSurface(TrkPlanarSurface *)
void setSmootherGain(double A[5][5])
int cost(std::vector< std::string > &files, node &n, const std::string &directory="", bool deleteref=false, bool relocate=false)
Definition hcg.cxx:922
struct color C
list descr
print "%s.properties()" % self.__name__

◆ finalize()

StatusCode TrigInDetTrackFitter::finalize ( )
virtual

Definition at line 91 of file TrigInDetTrackFitter.cxx.

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}
#define ATH_MSG_INFO(x)

◆ fit() [1/2]

void TrigInDetTrackFitter::fit ( const TrackCollection & inputTracks,
TrackCollection & fittedTracks,
const EventContext & ctx,
const Trk::ParticleHypothesis & matEffects = Trk::pion ) const
virtual

Implements ITrigInDetTrackFitter.

Definition at line 569 of file TrigInDetTrackFitter.cxx.

570{
572 fit(inputTracks,fittedTracks,tmp,ctx,matEffects,false);
573}
DataVector< Trk::Track > TrackCollection
This typedef represents a collection of Trk::Track objects.
void fit(const TrackCollection &, TrackCollection &, const EventContext &, const Trk::ParticleHypothesis &matEffects=Trk::pion) const

◆ fit() [2/2]

void TrigInDetTrackFitter::fit ( const TrackCollection & inputTracks,
TrackCollection & fittedTracks,
TrackCollection & fittedTrackswTP,
const EventContext & ctx,
const Trk::ParticleHypothesis & matEffects = Trk::pion,
const bool addTPtoTSoS = false ) const
virtual

Implements ITrigInDetTrackFitter.

Definition at line 575 of file TrigInDetTrackFitter.cxx.

576{
577 MagField::AtlasFieldCache fieldCache;
578
579 SG::ReadCondHandle<AtlasFieldCacheCondObj> fieldCondObj{m_fieldCondObjInputKey, ctx};
580 if (!fieldCondObj.isValid()) {
581 ATH_MSG_ERROR("Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
582 return;
583 }
584
585 fieldCondObj->getInitializedCache (fieldCache);
586 fittedTracks.reserve(inputTracks.size());
587 if( addTPtoTSoS ) fittedTrackswTP.reserve(inputTracks.size());
588 for(auto trIt = inputTracks.begin(); trIt != inputTracks.end(); ++trIt) {
589 Trk::Track* fittedTrack = nullptr;
590 Trk::Track* fittedTrackwTP = nullptr;
591 std::tie(fittedTrack,fittedTrackwTP) = fitTrack(**trIt, fieldCache, ctx, matEffects, addTPtoTSoS);
592 if (fittedTrack!=nullptr) {
593 fittedTracks.push_back(fittedTrack);
594 }
595 if (addTPtoTSoS && fittedTrackwTP!=nullptr) {
596 fittedTrackswTP.push_back(fittedTrackwTP);
597 }
598 }
599}
#define ATH_MSG_ERROR(x)
void reserve(size_type n)
Attempt to preallocate enough memory for a specified number of elements.
value_type push_back(value_type pElem)
Add an element to the end of the collection.
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
size_type size() const noexcept
Returns the number of elements in the collection.
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
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_fieldCondObjInputKey

◆ fitTrack()

std::pair< Trk::Track *, Trk::Track * > TrigInDetTrackFitter::fitTrack ( const Trk::Track & recoTrack,
MagField::AtlasFieldCache & fieldCache,
const EventContext & ctx,
const Trk::ParticleHypothesis & matEffects = Trk::pion,
const bool addTPtoTSoS = false ) const

Definition at line 601 of file TrigInDetTrackFitter.cxx.

601 {
602
603 const Trk::TrackParameters* trackPars = recoTrack.perigeeParameters();
604 if(trackPars==nullptr) {
605 ATH_MSG_WARNING("Fit Failed -- TrkTrack has no parameters");
606 return std::make_pair(nullptr,nullptr);
607 }
608
609 // 1. Create initial track state:
610 double Rk[5];
611 Rk[0] = trackPars->parameters()[Trk::d0];
612 Rk[1] = trackPars->parameters()[Trk::z0];
613 Rk[2] = trackPars->parameters()[Trk::phi0];
614 if(Rk[2]>M_PI) Rk[2]-=2*M_PI;
615 if(Rk[2]<-M_PI) Rk[2]+=2*M_PI;
616 double trk_theta = trackPars->parameters()[Trk::theta];
617 Rk[3] = trk_theta;
618 double trk_qOverP = trackPars->parameters()[Trk::qOverP];
619 Rk[4] = 1000.0*trk_qOverP;//MeV->GeV
620 double trk_Pt = sin(trk_theta)/trk_qOverP;
621
622 if(fabs(trk_Pt)<100.0)
623 {
624 ATH_MSG_DEBUG("Estimated Pt is too low "<<trk_Pt<<" - skipping fit");
625 return std::make_pair(nullptr,nullptr);
626 }
627
628 // 2. Create filtering nodes
629
630 std::vector<Trk::TrkBaseNode*> vpTrkNodes;
631 std::vector<Trk::TrkTrackState*> vpTrackStates;
632 vpTrackStates.reserve(vpTrkNodes.size() + 1);
633 bool trackResult = m_trackMaker->createDkfTrack(recoTrack,vpTrkNodes, m_DChi2);
634 int nHits=vpTrkNodes.size();
635 ATH_MSG_VERBOSE(nHits<<" filtering nodes created");
636
637 if(!trackResult) return std::make_pair(nullptr,nullptr);
638
639 // 3. Main algorithm: filter and smoother (Rauch-Tung-Striebel)
641 Trk::TrkTrackState* pTS = new Trk::TrkTrackState(Rk);
642 double Gk[5][5] = {{100.0, 0, 0, 0, 0},
643 {0, 100.0, 0, 0, 0},
644 {0, 0, 0.01, 0, 0},
645 {0, 0, 0, 0.01, 0},
646 {0, 0, 0, 0, 0.1}};
647 pTS->setTrackCovariance(Gk);
648 if(m_doMultScatt)
649 pTS->setScatteringMode(1);
650 if(m_doBremm)
651 pTS->setScatteringMode(2);
652 vpTrackStates.push_back(pTS);
653
654 //ATH_MSG_DEBUG("Initial chi2: "<<recoTrack.chi2()<<" track authorId: "<<recoTrack.algorithmId());
655 ATH_MSG_VERBOSE("Initial params: locT="<<Rk[0]<<" locL="<<Rk[1]<<" phi="<<Rk[2]
656 <<" theta="<<Rk[3]<<" Q="<<Rk[4]<<" pT="<<sin(Rk[3])/Rk[4]<<" GeV");
657
658 bool OK=true;
659
660 double chi2tot=0.0;
661 int ndoftot=-5;
662
663 Trk::TrkPlanarSurface* pSB=nullptr;
664 Trk::TrkPlanarSurface* pSE=nullptr;
665 for(auto pnIt = vpTrkNodes.begin(); pnIt!=vpTrkNodes.end(); ++pnIt) {
666 pSE=(*pnIt)->getSurface();
667 Trk::TrkTrackState* pNS=extrapolate(pTS,pSB,pSE,fieldCache);
668
669 pSB=pSE;
670 if(pNS!=nullptr) {
671 vpTrackStates.push_back(pNS);
672
673 (*pnIt)->validateMeasurement(pNS);
674 ATH_MSG_VERBOSE("dChi2="<<(*pnIt)->getChi2());
675 if((*pnIt)->isValidated())
676 {
677 chi2tot+=(*pnIt)->getChi2();
678 ndoftot+=(*pnIt)->getNdof();
679 }
680 (*pnIt)->updateTrackState(pNS);
681 pTS=pNS;
682 double est_Pt = 1000.0*sin(pTS->getTrackState(3))/pTS->getTrackState(4);
683 if(fabs(est_Pt)<200.0)
684 {
685 ATH_MSG_VERBOSE("Estimated Pt is too low "<<est_Pt<<" - skipping fit");
687 OK=false;break;
688 }
689 }
690 else
691 {
693 OK=false;break;
694 }
695 }
696 Trk::Track* fittedTrack = nullptr;
697 Trk::Track* fittedTrackwTP = nullptr;
698 if(OK)
699 {
700 for(auto ptsIt = vpTrackStates.rbegin();ptsIt!=vpTrackStates.rend();++ptsIt)
701 {
702 (*ptsIt)->runSmoother();
703 }
704 pTS=(*vpTrackStates.begin());
705 //correct GeV->MeV
706
707 correctScale(pTS);
708
709 double qOverP = pTS->getTrackState(4);
710 double pt=sin(pTS->getTrackState(3))/pTS->getTrackState(4);
711 double phi0 = pTS->getTrackState(2);
712 if(phi0>M_PI) phi0-=2*M_PI;
713 if(phi0<-M_PI) phi0+=2*M_PI;
714 double theta = pTS->getTrackState(3);
715 double z0 = pTS->getTrackState(1);
716 double d0 = pTS->getTrackState(0);
717 bool bad_cov = false;
718 auto cov = AmgSymMatrix(5){};
719 for(int i=0;i<5;i++) {
720 double cov_diag = pTS->getTrackCovariance(i,i);
721 if (cov_diag < 0) {
722 bad_cov = true;//Diagonal elements must be positive
723 ATH_MSG_DEBUG("REGTEST: cov(" << i << "," << i << ") =" << cov_diag << " < 0, reject track");
724 break;
725 }
726 (cov)(i, i) = pTS->getTrackCovariance(i,i);
727 for(int j=i+1;j<5;j++) {
728 cov.fillSymmetric(i, j, pTS->getTrackCovariance(i,j));
729 }
730 }
731
732 if((ndoftot<0) || (fabs(pt)<100.0) || (std::isnan(pt)) || bad_cov)
733 {
734 ATH_MSG_DEBUG("Fit failed - possibly floating point problem");
735 }
736 else
737 {
738 Trk::PerigeeSurface perigeeSurface;
739 auto perigee = std::make_unique<Trk::Perigee>(d0, z0, phi0, theta, qOverP, perigeeSurface, cov);
740 ATH_MSG_VERBOSE("perigee: " << *perigee);
741 std::unique_ptr<Trk::Perigee> perigeewTP = (addTPtoTSoS) ? std::make_unique<Trk::Perigee>(d0, z0, phi0, theta, qOverP, perigeeSurface, cov) : nullptr;
742
743 std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern;
744 typePattern.set(Trk::TrackStateOnSurface::Perigee);
745 auto pParVec = std::make_unique<Trk::TrackStates>();
746 auto pParVecwTP = std::make_unique<Trk::TrackStates>();
748 pParVec->reserve(vpTrkNodes.size()+1);
749 pParVec->push_back(new Trk::TrackStateOnSurface(nullptr, std::move(perigee),nullptr, typePattern));
750 if (addTPtoTSoS) {
751 std::bitset<
753 typePatternwTP;
754 typePatternwTP.set(Trk::TrackStateOnSurface::Perigee);
755 pParVecwTP->reserve(vpTrkNodes.size() + 1);
756 pParVecwTP->push_back(
757 new Trk::TrackStateOnSurface(nullptr, std::move(perigeewTP), nullptr, typePatternwTP));
758 }
759 for (auto pnIt = vpTrkNodes.begin(); pnIt != vpTrkNodes.end(); ++pnIt) {
760 if((*pnIt)->isValidated()) {
761 Trk::TrackStateOnSurface* pTSS = createTrackStateOnSurface(*pnIt,false, ctx);
762 Trk::TrackStateOnSurface* pTSSwTP = nullptr;
763 if( addTPtoTSoS ) pTSSwTP = createTrackStateOnSurface(*pnIt,true,ctx);
764 if(pTSS!=nullptr) {
765 pParVec->push_back(pTSS);
766 }
767 if(pTSSwTP!=nullptr) {
768 pParVecwTP->push_back(pTSSwTP);
769 }
770 }
771 }
772 }
773 else {
774 pParVec->reserve(recoTrack.trackStateOnSurfaces()->size());
775 pParVec->push_back(new Trk::TrackStateOnSurface(nullptr, std::move(perigee),nullptr, typePattern));
776
777 for (auto tSOS = recoTrack.trackStateOnSurfaces()->begin(); tSOS != recoTrack.trackStateOnSurfaces()->end(); ++tSOS) {
778 //Don't store perigee - new perigee created above
779 if ((*tSOS)->type(Trk::TrackStateOnSurface::Perigee) == false) {
780 pParVec->push_back((*tSOS)->clone());
781 }
782 }
783 }
784 ATH_MSG_VERBOSE("Total chi2 ="<<chi2tot<<" NDOF="<<ndoftot);
785 if(msgLvl(MSG::VERBOSE)) {
786 double eta = -log(tan(0.5*theta));
787 ATH_MSG_VERBOSE("Fitted parameters: d0="<<d0<<" phi0="<<phi0<<" z0="<<z0
788 <<" eta0="<<eta<<" pt="<<pt);
789 }
790 auto pFQ= std::make_unique<Trk::FitQuality>(chi2tot,ndoftot);
791 Trk::TrackInfo info(recoTrack.info());
792 info.setParticleHypothesis(matEffects);
793 fittedTrack = new Trk::Track(info, std::move(pParVec), std::move(pFQ));//fittedTrack now owns pParVec and pFQ
794 if( addTPtoTSoS ) {
795 auto pFQwTP=std::make_unique<Trk::FitQuality>(chi2tot,ndoftot);
796 Trk::TrackInfo infowTP(recoTrack.info());
797 infowTP.setParticleHypothesis(matEffects);
798 fittedTrackwTP = new Trk::Track(infowTP, std::move(pParVecwTP), std::move(pFQwTP));//fittedTrack now owns pParVecwTP and pFQwTP
799 }
800 }
801 }
802 else
803 {
804 ATH_MSG_DEBUG("Forward Kalman filter: extrapolation failure ");
805 }
806
807 for(auto pnIt = vpTrkNodes.begin(); pnIt!=vpTrkNodes.end(); ++pnIt) {
808 delete((*pnIt)->getSurface());
809 delete (*pnIt);
810 }
811 vpTrkNodes.clear();
812 for(auto ptsIt = vpTrackStates.begin();ptsIt!=vpTrackStates.end();++ptsIt) {
813 delete (*ptsIt);
814 }
815 vpTrackStates.clear();
816
817 return std::make_pair(fittedTrack,fittedTrackwTP);
818}
Scalar eta() const
pseudorapidity method
Scalar theta() const
theta method
#define ATH_MSG_VERBOSE(x)
static const uint32_t nHits
bool msgLvl(const MSG::Level lvl) const
Trk::TrkTrackState * extrapolate(Trk::TrkTrackState *, Trk::TrkPlanarSurface *, Trk::TrkPlanarSurface *, MagField::AtlasFieldCache &) const
Trk::TrackStateOnSurface * createTrackStateOnSurface(Trk::TrkBaseNode *pN, const bool, const EventContext &ctx) const
void correctScale(Trk::TrkTrackState *) const
@ Perigee
This represents a perigee, and so will contain a Perigee object only.
const Trk::TrackStates * trackStateOnSurfaces() const
return a pointer to a const DataVector of const TrackStateOnSurfaces.
const TrackInfo & info() const
Returns a const ref to info of a const tracks.
const Perigee * perigeeParameters() const
return Perigee.
@ qOverP
perigee
@ phi0
Definition ParamDefs.h:65
@ theta
Definition ParamDefs.h:66
@ qOverP
perigee
Definition ParamDefs.h:67
@ d0
Definition ParamDefs.h:63
@ z0
Definition ParamDefs.h:64
ParametersBase< TrackParametersDim, Charged > TrackParameters

◆ getMagneticField()

void TrigInDetTrackFitter::getMagneticField ( double r[3],
double * B,
MagField::AtlasFieldCache & fieldCache ) const
private

Definition at line 104 of file TrigInDetTrackFitter.cxx.

104 {
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}
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,...
int r
Definition globals.cxx:22

◆ getUnbiasedResiduals()

StatusCode TrigInDetTrackFitter::getUnbiasedResiduals ( const Trk::Track & pT,
std::vector< TrigL2HitResidual > & vResid,
const EventContext & ctx ) const
virtual

Implements ITrigInDetTrackFitter.

Definition at line 891 of file TrigInDetTrackFitter.cxx.

892 {
893
894 const Trk::TrackParameters* trackPars = pT.perigeeParameters();
895 if(trackPars==nullptr) {
896 ATH_MSG_WARNING("Fit Failed -- TrkTrack has no parameters");
897 return StatusCode::FAILURE;
898 }
899
900 MagField::AtlasFieldCache fieldCache;
901
902 SG::ReadCondHandle<AtlasFieldCacheCondObj> fieldCondObj{m_fieldCondObjInputKey, ctx};
903 if (!fieldCondObj.isValid()) {
904 ATH_MSG_ERROR("Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
905 return StatusCode::FAILURE;
906 }
907
908 fieldCondObj->getInitializedCache (fieldCache);
909 std::vector<Trk::TrkBaseNode*> vpTrkNodes;
910 std::vector<Trk::TrkTrackState*> vpTrackStates;
911 vResid.clear();
912 double trk_theta = trackPars->parameters()[Trk::theta];
913 double trk_qOverP = trackPars->parameters()[Trk::qOverP];
914 double Pt = sin(trk_theta)/trk_qOverP;
915 if(fabs(Pt)<100.0)
916 {
917 ATH_MSG_DEBUG("TrigL2ResidualCalculator failed -- Estimated Pt is too low "<<Pt);
918 return StatusCode::FAILURE;
919 }
920
921 // 1. Create filtering nodes
922
923 bool trackResult = m_trackMaker->createDkfTrack(pT,vpTrkNodes, m_DChi2);
924 if(!trackResult)
925 {
926 ATH_MSG_DEBUG("TrigDkfTrackMaker failed");
927 return StatusCode::FAILURE;
928 }
929
930 bool OK=true;
931 std::vector<Trk::TrkBaseNode*>::iterator pnIt,pnEnd(vpTrkNodes.end());
932
933 for(std::vector<Trk::TrkBaseNode*>::iterator pNodeIt=vpTrkNodes.begin();pNodeIt!=vpTrkNodes.end();
934 ++pNodeIt)
935 {
936 Trk::TrkBaseNode* pMaskedNode=(*pNodeIt);
937 Trk::TrkTrackState* pMaskedState=nullptr;
938
939 // 2. Create initial track state:
940
941 double Rk[5];
942 Rk[0] = trackPars->parameters()[Trk::d0];
943 Rk[1] = trackPars->parameters()[Trk::z0];
944 Rk[2] = trackPars->parameters()[Trk::phi0];
945 if(Rk[2]>M_PI) Rk[2]-=2*M_PI;
946 if(Rk[2]<-M_PI) Rk[2]+=2*M_PI;
947 trk_theta = trackPars->parameters()[Trk::theta];
948 Rk[3] = trk_theta;
949 Rk[4] = 1000.0*trackPars->parameters()[Trk::qOverP];//MeV->GeV
950 //No need to correct scale back - not returning track
951
952 // 3. Main algorithm: filter and smoother (Rauch-Tung-Striebel)
953
954 Trk::TrkTrackState* pTS = new Trk::TrkTrackState(Rk);
955 double Gk[5][5] = {{100.0, 0, 0, 0, 0},
956 {0, 100.0, 0, 0, 0},
957 {0, 0, 0.01, 0, 0},
958 {0, 0, 0, 0.01, 0},
959 {0, 0, 0, 0, 0.1}};
960 pTS->setTrackCovariance(Gk);
961 if(m_doMultScatt)
962 pTS->setScatteringMode(1);
963 if(m_doBremm)
964 pTS->setScatteringMode(2);
965 vpTrackStates.push_back(pTS);
966
967 ATH_MSG_DEBUG("Initial params: locT="<<Rk[0]<<" locL="<<Rk[1]<<" phi="<<Rk[2]
968 <<" theta="<<Rk[3]<<" Q="<<Rk[4]<<" pT="<<sin(Rk[3])/Rk[4]);
969
970 OK=true;
971 Trk::TrkPlanarSurface *pSB=nullptr,*pSE;
972
973 for(pnIt=vpTrkNodes.begin();pnIt!=pnEnd;++pnIt)
974 {
975 pSE=(*pnIt)->getSurface();
976 Trk::TrkTrackState* pNS=extrapolate(pTS,pSB,pSE,fieldCache);
977
978 pSB=pSE;
979 if(pNS!=nullptr)
980 {
981 vpTrackStates.push_back(pNS);
982
983 (*pnIt)->validateMeasurement(pNS);
984 ATH_MSG_DEBUG("dChi2="<<(*pnIt)->getChi2());
985 if((*pnIt)!=pMaskedNode)
986 {
987 (*pnIt)->updateTrackState(pNS);
988 }
989 else
990 {
991 pMaskedState=pNS;
992 }
993 pTS=pNS;
994 Pt=sin(pTS->getTrackState(3))/pTS->getTrackState(4);
995 if(fabs(Pt)<0.2)
996 {
997 ATH_MSG_DEBUG("Estimated Pt is too low "<<Pt<<" - skipping fit");
998 OK=false;break;
999 }
1000 }
1001 else
1002 {
1003 OK=false;break;
1004 }
1005 }
1006 if(OK)
1007 {
1008 std::vector<Trk::TrkTrackState*>::reverse_iterator ptsrIt(vpTrackStates.rbegin()),
1009 ptsrEnd(vpTrackStates.rend());
1010
1011 for(;ptsrIt!=ptsrEnd;++ptsrIt)
1012 {
1013 (*ptsrIt)->runSmoother();
1014 }
1015
1016 pMaskedNode->validateMeasurement(pMaskedState);
1017
1018 double r[2],V[2][2];
1019
1020 int nSize=pMaskedNode->getResiduals(r);
1021 nSize=pMaskedNode->getInverseResidualVariance(V);
1022 const Trk::PrepRawData* pPRD = pMaskedNode->getPrepRawData();
1023
1024 Identifier id = pPRD->identify();
1025
1026 Region region = Region::Undefined;
1027 if(m_idHelper->is_pixel(id))
1028 {
1029 region=(m_pixelId->is_barrel(id)) ? Region::PixBarrel: Region::PixEndcap;
1030 if (m_pixelId->is_blayer(id)) {
1031 region = Region::IBL;
1032 }
1033 }
1034 if(m_idHelper->is_sct(id))
1035 {
1036 region=(m_sctId->is_barrel(id)) ? Region::SctBarrel : Region::SctEndcap;
1037 }
1038 if(nSize==1) {
1039 if(V[0][0]>0.0) {
1040 vResid.push_back(TrigL2HitResidual(id,region,r[0],r[0]*sqrt(V[0][0])));
1041 }
1042 else {
1043 OK=false;
1044 break;
1045 }
1046 }
1047 else {
1048 if((V[0][0]>0.0) && (V[1][1]>0.0)) {
1049 vResid.push_back(TrigL2HitResidual(id,region,r[0],r[0]*sqrt(V[0][0]),
1050 r[1],r[1]*sqrt(V[1][1])));
1051 }
1052 else {
1053 OK=false;
1054 break;
1055 }
1056 }
1057 }
1058 else
1059 {
1060 ATH_MSG_DEBUG("Forward Kalman filter: extrapolation failure ");
1061 vResid.clear();
1062 }
1063 for(std::vector<Trk::TrkTrackState*>::iterator ptsIt=vpTrackStates.begin();
1064 ptsIt!=vpTrackStates.end();++ptsIt) delete (*ptsIt);
1065 vpTrackStates.clear();
1066 if(!OK) break;
1067 }
1068 pnIt=vpTrkNodes.begin();pnEnd=vpTrkNodes.end();
1069 for(;pnIt!=pnEnd;++pnIt)
1070 {
1071 delete((*pnIt)->getSurface());
1072 delete (*pnIt);
1073 }
1074 vpTrkNodes.clear();
1075
1076 if(OK) return StatusCode::SUCCESS;
1077 else return StatusCode::FAILURE;
1078
1079}
const AtlasDetectorID * m_idHelper
virtual int getInverseResidualVariance(double[2][2])=0
virtual void validateMeasurement(TrkTrackState *)=0
virtual int getResiduals(double[2])=0

◆ initialize()

StatusCode TrigInDetTrackFitter::initialize ( )
virtual

Definition at line 76 of file TrigInDetTrackFitter.cxx.

77{
78 ATH_CHECK(m_trackMaker.retrieve());
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}
#define ATH_CHECK
Evaluate an expression and check for errors.
const ServiceHandle< StoreGateSvc > & detStore() const
retrieve(aClass, aKey=None)
Definition PyKernel.py:110

◆ inputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< AlgTool > >::inputHandles ( ) const
overridevirtualinherited

Return this algorithm's input handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

◆ interfaceID()

const InterfaceID & ITrigInDetTrackFitter::interfaceID ( )
inlinestaticinherited

Definition at line 24 of file ITrigInDetTrackFitter.h.

24 {
26 }
static const InterfaceID IID_ITrigInDetTrackFitter("ITrigInDetTrackFitter", 1, 0)

◆ msg()

MsgStream & AthCommonMsg< AlgTool >::msg ( ) const
inlineinherited

Definition at line 24 of file AthCommonMsg.h.

24 {
25 return this->msgStream();
26 }

◆ msgLvl()

bool AthCommonMsg< AlgTool >::msgLvl ( const MSG::Level lvl) const
inlineinherited

Definition at line 30 of file AthCommonMsg.h.

30 {
31 return this->msgLevel(lvl);
32 }

◆ outputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< AlgTool > >::outputHandles ( ) const
overridevirtualinherited

Return this algorithm's output handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

◆ renounce()

std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > AthCommonDataStore< AthCommonMsg< AlgTool > >::renounce ( T & h)
inlineprotectedinherited

Definition at line 380 of file AthCommonDataStore.h.

381 {
382 h.renounce();
384 }
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > renounce(T &h)

◆ renounceArray()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::renounceArray ( SG::VarHandleKeyArray & handlesArray)
inlineprotectedinherited

remove all handles from I/O resolution

Definition at line 364 of file AthCommonDataStore.h.

364 {
366 }

◆ sysInitialize()

virtual StatusCode AthCommonDataStore< AthCommonMsg< AlgTool > >::sysInitialize ( )
overridevirtualinherited

Perform system initialization for an algorithm.

We override this to declare all the elements of handle key arrays at the end of initialization. See comments on updateVHKA.

Reimplemented in asg::AsgMetadataTool, AthCheckedComponent< AthAlgTool >, AthCheckedComponent<::AthAlgTool >, and DerivationFramework::CfAthAlgTool.

◆ sysStart()

virtual StatusCode AthCommonDataStore< AthCommonMsg< AlgTool > >::sysStart ( )
overridevirtualinherited

Handle START transition.

We override this in order to make sure that conditions handle keys can cache a pointer to the conditions container.

◆ updateVHKA()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::updateVHKA ( Gaudi::Details::PropertyBase & )
inlineinherited

Definition at line 308 of file AthCommonDataStore.h.

308 {
309 // debug() << "updateVHKA for property " << p.name() << " " << p.toString()
310 // << " size: " << m_vhka.size() << endmsg;
311 for (auto &a : m_vhka) {
313 for (auto k : keys) {
314 k->setOwner(this);
315 }
316 }
317 }
std::vector< SG::VarHandleKeyArray * > m_vhka

Member Data Documentation

◆ m_correctClusterPos

bool TrigInDetTrackFitter::m_correctClusterPos
private

Definition at line 65 of file TrigInDetTrackFitter.h.

◆ m_DChi2

double TrigInDetTrackFitter::m_DChi2
private

Definition at line 62 of file TrigInDetTrackFitter.h.

◆ m_detStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< AlgTool > >::m_detStore
privateinherited

Pointer to StoreGate (detector store by default)

Definition at line 393 of file AthCommonDataStore.h.

◆ m_doBremm

bool TrigInDetTrackFitter::m_doBremm
private

Definition at line 64 of file TrigInDetTrackFitter.h.

◆ m_doMultScatt

bool TrigInDetTrackFitter::m_doMultScatt
private

Definition at line 63 of file TrigInDetTrackFitter.h.

◆ m_evtStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< AlgTool > >::m_evtStore
privateinherited

Pointer to StoreGate (event store by default)

Definition at line 390 of file AthCommonDataStore.h.

◆ m_fieldCondObjInputKey

SG::ReadCondHandleKey<AtlasFieldCacheCondObj> TrigInDetTrackFitter::m_fieldCondObjInputKey {this, "AtlasFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"}
private

Definition at line 68 of file TrigInDetTrackFitter.h.

68{this, "AtlasFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};

◆ m_fitErrorsDivergence

std::atomic<size_t> TrigInDetTrackFitter::m_fitErrorsDivergence
mutableprivate

Definition at line 59 of file TrigInDetTrackFitter.h.

◆ m_fitErrorsLowPt

std::atomic<size_t> TrigInDetTrackFitter::m_fitErrorsLowPt
mutableprivate

Definition at line 60 of file TrigInDetTrackFitter.h.

◆ m_fitErrorsUnresolved

std::atomic<size_t> TrigInDetTrackFitter::m_fitErrorsUnresolved
mutableprivate

Definition at line 58 of file TrigInDetTrackFitter.h.

◆ m_idHelper

const AtlasDetectorID* TrigInDetTrackFitter::m_idHelper = nullptr
private

Definition at line 71 of file TrigInDetTrackFitter.h.

◆ m_nTracksTotal

std::atomic<size_t> TrigInDetTrackFitter::m_nTracksTotal
mutableprivate

Definition at line 57 of file TrigInDetTrackFitter.h.

◆ m_pixelId

const PixelID* TrigInDetTrackFitter::m_pixelId = nullptr
private

Definition at line 69 of file TrigInDetTrackFitter.h.

◆ m_ROTcreator

ToolHandle<Trk::IRIO_OnTrackCreator> TrigInDetTrackFitter::m_ROTcreator
private

Definition at line 67 of file TrigInDetTrackFitter.h.

◆ m_sctId

const SCT_ID* TrigInDetTrackFitter::m_sctId = nullptr
private

Definition at line 70 of file TrigInDetTrackFitter.h.

◆ m_trackMaker

ToolHandle<ITrigDkfTrackMakerTool> TrigInDetTrackFitter::m_trackMaker
private

Definition at line 66 of file TrigInDetTrackFitter.h.

◆ m_varHandleArraysDeclared

bool AthCommonDataStore< AthCommonMsg< AlgTool > >::m_varHandleArraysDeclared
privateinherited

Definition at line 399 of file AthCommonDataStore.h.

◆ m_vhka

std::vector<SG::VarHandleKeyArray*> AthCommonDataStore< AthCommonMsg< AlgTool > >::m_vhka
privateinherited

Definition at line 398 of file AthCommonDataStore.h.


The documentation for this class was generated from the following files: