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 55 of file TrigInDetTrackFitter.cxx.

57 : 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" );
71}
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 108 of file TrigInDetTrackFitter.cxx.

108 {
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}
double getTrackState(int i)
void setTrackState(const double A[5])
double getTrackCovariance(int i, int j)
void setTrackCovariance(double A[5][5])
float j(const xAOD::IParticle &, const xAOD::TrackMeasurementValidation &hit, const Eigen::Matrix3d &jab_inv)

◆ createTrackStateOnSurface()

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

Definition at line 817 of file TrigInDetTrackFitter.cxx.

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;
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}
#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()
virtual char getNodeType()
int getNdof() const
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 130 of file TrigInDetTrackFitter.cxx.

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;
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)
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}
#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 *)
double getRotMatrix(int, int)
double getInvRotMatrix(int, int)
void transformPointToGlobal(const double *, double *)
void setPreviousState(TrkTrackState *)
void setSmootherGain(double A[5][5])
void attachToSurface(TrkPlanarSurface *)
int cost(std::vector< std::string > &files, node &n, const std::string &directory="", bool deleteref=false, bool relocate=false)
Definition hcg.cxx:926
struct color C
list descr
print "%s.properties()" % self.__name__

◆ finalize()

StatusCode TrigInDetTrackFitter::finalize ( )
virtual

Definition at line 88 of file TrigInDetTrackFitter.cxx.

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}
#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 566 of file TrigInDetTrackFitter.cxx.

567{
569 fit(inputTracks,fittedTracks,tmp,ctx,matEffects,false);
570}
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 572 of file TrigInDetTrackFitter.cxx.

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}
#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 598 of file TrigInDetTrackFitter.cxx.

598 {
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)
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");
684 OK=false;break;
685 }
686 }
687 else
688 {
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>();
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<
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}
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 101 of file TrigInDetTrackFitter.cxx.

101 {
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}
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 888 of file TrigInDetTrackFitter.cxx.

889 {
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}
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 73 of file TrigInDetTrackFitter.cxx.

74{
75 ATH_CHECK(m_trackMaker.retrieve());
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}
#define ATH_CHECK
Evaluate an expression and check for errors.
const ServiceHandle< StoreGateSvc > & detStore() const

◆ 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 >, and AthCheckedComponent<::AthAlgTool >.

◆ 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: