FairRoot/PandaRoot
Public Member Functions | Protected Attributes | Private Member Functions | Private Attributes | List of all members
RKTrackRep Class Reference

Track Representation module based on a Runge-Kutta algorithm including a full material model. More...

#include <RKTrackRep.h>

Inheritance diagram for RKTrackRep:
GFAbsTrackRep

Public Member Functions

 RKTrackRep ()
 
 RKTrackRep (const TVector3 &pos, const TVector3 &mom, const TVector3 &poserr, const TVector3 &momerr, const int &PDGCode)
 
 RKTrackRep (const GFTrackCand *aGFTrackCandPtr)
 
 RKTrackRep (const TVector3 &pos, const TVector3 &mom, const int &PDGCode)
 
 RKTrackRep (const GFDetPlane &pl, const TVector3 &mom, const int &PDGCode)
 
virtual ~RKTrackRep ()
 
virtual GFAbsTrackRepclone () const
 
virtual GFAbsTrackRepprototype () const
 
double extrapolate (const GFDetPlane &, TMatrixT< double > &statePred, TMatrixT< double > &covPred)
 returns the tracklength spanned in this extrapolation More...
 
double extrapolate (const GFDetPlane &, TMatrixT< double > &statePred)
 returns the tracklength spanned in this extrapolation More...
 
void extrapolateToPoint (const TVector3 &pos, TVector3 &poca, TVector3 &dirInPoca)
 This method is to extrapolate the track to point of closest approach to a point in space. More...
 
void extrapolateToLine (const TVector3 &point1, const TVector3 &point2, TVector3 &poca, TVector3 &dirInPoca, TVector3 &poca_onwire)
 This method extrapolates to the point of closest approach to a line. More...
 
double stepalong (double h, TVector3 &point, TVector3 &dir)
 make step of h cm along the track, returns the tracklength spanned in this extrapolation More...
 
TVector3 getPos (const GFDetPlane &)
 Returns position of the track in the plane. More...
 
TVector3 getMom (const GFDetPlane &)
 Returns momentum of the track in the plane. More...
 
void getPosMom (const GFDetPlane &, TVector3 &pos, TVector3 &mom)
 Gets position and momentum in the plane. More...
 
double getCharge () const
 Returns charge. More...
 
int getPDG ()
 
void switchDirection ()
 deprecated More...
 
void setPDG (int)
 Set PDG particle code. More...
 
void setData (const TMatrixT< double > &st, const GFDetPlane &pl, const TMatrixT< double > *cov=NULL, const TMatrixT< double > *aux=NULL)
 Sets state, plane and (optionally) covariance. More...
 
const TMatrixT< double > * getAuxInfo (const GFDetPlane &pl)
 Get auxillary information from the track representation. More...
 
bool hasAuxInfo ()
 See if the track representation has auxillary information stored. More...
 
void getPosMomCov (const GFDetPlane &pl, TVector3 &pos, TVector3 &mom, TMatrixT< double > &cov)
 method which gets position, momentum and 6x6 covariance matrix More...
 
double extrapolate (const GFDetPlane &plane)
 This changes the state and cov and plane of the rep. More...
 
unsigned int getDim () const
 returns dimension of state vector More...
 
virtual void Print () const
 
TMatrixT< double > getState () const
 
TMatrixT< double > getCov () const
 
double getStateElem (int i) const
 
double getCovElem (int i, int j) const
 
TVector3 getPos ()
 
TVector3 getMom ()
 
void getPosMomCov (TVector3 &pos, TVector3 &mom, TMatrixT< double > &c)
 
TMatrixT< double > getFirstState () const
 
TMatrixT< double > getFirstCov () const
 
GFDetPlane getFirstPlane () const
 
TMatrixT< double > getLastState () const
 
TMatrixT< double > getLastCov () const
 
GFDetPlane getLastPlane () const
 
double getChiSqu () const
 
double getRedChiSqu () const
 returns chi2/ndf More...
 
unsigned int getNDF () const
 
void setCov (const TMatrixT< double > &aCov)
 
void setFirstState (const TMatrixT< double > &aState)
 
void setFirstCov (const TMatrixT< double > &aCov)
 
void setFirstPlane (const GFDetPlane &aPlane)
 
void setLastState (const TMatrixT< double > &aState)
 
void setLastCov (const TMatrixT< double > &aCov)
 
void setLastPlane (const GFDetPlane &aPlane)
 
const GFDetPlanegetReferencePlane () const
 
void setChiSqu (double aChiSqu)
 
void setNDF (unsigned int n)
 
void addChiSqu (double aChiSqu)
 
void addNDF (unsigned int n)
 
void setStatusFlag (int _val)
 
bool setInverted (bool f=true)
 Deprecated. Should be removed soon. More...
 
bool getStatusFlag ()
 
virtual void reset ()
 

Protected Attributes

unsigned int fDimension
 Dimensionality of track representation. More...
 
TMatrixT< double > fState
 The vector of track parameters. More...
 
TMatrixT< double > fCov
 The covariance matrix. More...
 
double fChiSqu
 chiSqu of the track fit More...
 
unsigned int fNdf
 
int fStatusFlag
 status of track representation: 0 means everything's OK More...
 
bool fInverted
 specifies the direction of flight of the particle More...
 
TMatrixT< double > fFirstState
 state, cov and plane for first and last point in fit More...
 
TMatrixT< double > fFirstCov
 
TMatrixT< double > fLastState
 
TMatrixT< double > fLastCov
 
GFDetPlane fFirstPlane
 
GFDetPlane fLastPlane
 
GFDetPlane fRefPlane
 

Private Member Functions

RKTrackRepoperator= (const RKTrackRep *)
 
bool RKutta (const GFDetPlane &plane, double *P, double &coveredDistance, std::vector< TVector3 > &points, std::vector< double > &pointLengths, const double &maxLen=-1, bool calcCov=true) const
 Propagates the particle through the magnetic field. More...
 
TVector3 poca2Line (const TVector3 &extr1, const TVector3 &extr2, const TVector3 &point) const
 
double Extrap (const GFDetPlane &plane, TMatrixT< double > *state, TMatrixT< double > *cov=NULL) const
 Handles propagation and material effects. More...
 

Private Attributes

bool fDirection
 
int fPdg
 PDG particle code. More...
 
double fMass
 Mass (in GeV) More...
 
double fCharge
 Charge. More...
 
GFDetPlane fCachePlane
 
double fCacheSpu
 
double fSpu
 
TMatrixT< double > fAuxInfo
 

Detailed Description

Track Representation module based on a Runge-Kutta algorithm including a full material model.

Author
Christian Höppner (Technische Universität München, original author)
Johannes Rauch (Technische Universität München, author)

The Runge Kutta implementation stems from GEANT3 originally (R. Brun et al.). Porting to C goes back to Igor Gavrilenko @ CERN. The code was taken from the Phast analysis package of the COMPASS experiment (Sergei Gerrassimov @ CERN).

Definition at line 49 of file tracking/GenfitTools/trackrep/RKTrackRep/RKTrackRep.h.

Constructor & Destructor Documentation

RKTrackRep::RKTrackRep ( )
RKTrackRep::RKTrackRep ( const TVector3 &  pos,
const TVector3 &  mom,
const TVector3 &  poserr,
const TVector3 &  momerr,
const int &  PDGCode 
)

Definition at line 78 of file RKTrackRep.cxx.

References fCharge, GFAbsTrackRep::fCov, GFAbsTrackRep::fRefPlane, fSpu, GFAbsTrackRep::fState, GFDetPlane::getO(), GFDetPlane::getU(), GFDetPlane::getV(), GFDetPlane::setNormal(), GFDetPlane::setO(), setPDG(), and v.

82  :
83  GFAbsTrackRep(5), fDirection(true), fCachePlane(), fCacheSpu(1), fAuxInfo(1,1) {
84  setPDG(PDGCode); // also sets charge and mass
85 
86 
89 
90  fState[0][0]=fCharge/mom.Mag();
91  //u' and v'
92  fState[1][0]=0.0;
93  fState[2][0]=0.0;
94  //u and v
95  fState[3][0]=0.0;
96  fState[4][0]=0.0;
97  //spu
98  fSpu=1.;
99 
100  TVector3 o=fRefPlane.getO();
101  TVector3 u=fRefPlane.getU();
102  TVector3 v=fRefPlane.getV();
103  TVector3 w=u.Cross(v);
104  double pu=0.;
105  double pv=0.;
106  double pw=mom.Mag();
107 
108  fCov[3][3] = poserr.X()*poserr.X() * u.X()*u.X() +
109  poserr.Y()*poserr.Y() * u.Y()*u.Y() +
110  poserr.Z()*poserr.Z() * u.Z()*u.Z();
111  fCov[4][4] = poserr.X()*poserr.X() * v.X()*v.X() +
112  poserr.Y()*poserr.Y() * v.Y()*v.Y() +
113  poserr.Z()*poserr.Z() * v.Z()*v.Z();
114  fCov[0][0] = fCharge*fCharge/pow(mom.Mag(),6.) *
115  (mom.X()*mom.X() * momerr.X()*momerr.X()+
116  mom.Y()*mom.Y() * momerr.Y()*momerr.Y()+
117  mom.Z()*mom.Z() * momerr.Z()*momerr.Z());
118  fCov[1][1] = pow((u.X()/pw - w.X()*pu/(pw*pw)),2.)*momerr.X()*momerr.X() +
119  pow((u.Y()/pw - w.Y()*pu/(pw*pw)),2.)*momerr.Y()*momerr.Y() +
120  pow((u.Z()/pw - w.Z()*pu/(pw*pw)),2.)*momerr.Z()*momerr.Z();
121  fCov[2][2] = pow((v.X()/pw - w.X()*pv/(pw*pw)),2.)*momerr.X()*momerr.X() +
122  pow((v.Y()/pw - w.Y()*pv/(pw*pw)),2.)*momerr.Y()*momerr.Y() +
123  pow((v.Z()/pw - w.Z()*pv/(pw*pw)),2.)*momerr.Z()*momerr.Z();
124 
125 }
TVector3 pos
TVector3 getV() const
Definition: GFDetPlane.h:77
void setO(const TVector3 &o)
Definition: GFDetPlane.cxx:92
Double_t mom
Definition: plot_dirc.C:14
__m128 v
Definition: P4_F32vec4.h:4
void setNormal(TVector3 n)
Definition: GFDetPlane.cxx:156
GFDetPlane fRefPlane
TVector3 getU() const
Definition: GFDetPlane.h:76
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
void setPDG(int)
Set PDG particle code.
Definition: RKTrackRep.cxx:274
TVector3 getO() const
Definition: GFDetPlane.h:75
TMatrixT< double > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:91
RKTrackRep::RKTrackRep ( const GFTrackCand aGFTrackCandPtr)

Definition at line 128 of file RKTrackRep.cxx.

References fCharge, GFAbsTrackRep::fCov, GFAbsTrackRep::fRefPlane, fSpu, GFAbsTrackRep::fState, GFTrackCand::getDirError(), GFTrackCand::getDirSeed(), GFDetPlane::getO(), GFTrackCand::getPdgCode(), GFTrackCand::getPosError(), GFTrackCand::getPosSeed(), GFDetPlane::getU(), GFDetPlane::getV(), mom, GFDetPlane::setNormal(), GFDetPlane::setO(), setPDG(), and v.

128  :
129  GFAbsTrackRep(5), fDirection(true), fCachePlane(), fCacheSpu(1), fAuxInfo(1,1) {
130  setPDG(aGFTrackCandPtr->getPdgCode()); // also sets charge and mass
131 
132  TVector3 mom = aGFTrackCandPtr->getDirSeed();
133  fRefPlane.setO(aGFTrackCandPtr->getPosSeed());
134  fRefPlane.setNormal(mom);
135  double pw=mom.Mag();
136  fState[0][0]=fCharge/pw;
137  //u' and v'
138  fState[1][0]=0.0;
139  fState[2][0]=0.0;
140  //u and v
141  fState[3][0]=0.0;
142  fState[4][0]=0.0;
143  //spu
144  fSpu=1.;
145 
146  TVector3 o=fRefPlane.getO();
147  TVector3 u=fRefPlane.getU();
148  TVector3 v=fRefPlane.getV();
149  TVector3 w=u.Cross(v);
150  double pu=0.;
151  double pv=0.;
152 
153  TVector3 poserr = aGFTrackCandPtr->getPosError();
154  TVector3 momerr = aGFTrackCandPtr->getDirError();
155  fCov[3][3] = poserr.X()*poserr.X() * u.X()*u.X() +
156  poserr.Y()*poserr.Y() * u.Y()*u.Y() +
157  poserr.Z()*poserr.Z() * u.Z()*u.Z();
158  fCov[4][4] = poserr.X()*poserr.X() * v.X()*v.X() +
159  poserr.Y()*poserr.Y() * v.Y()*v.Y() +
160  poserr.Z()*poserr.Z() * v.Z()*v.Z();
161  fCov[0][0] = fCharge*fCharge/pow(mom.Mag(),6.) *
162  (mom.X()*mom.X() * momerr.X()*momerr.X()+
163  mom.Y()*mom.Y() * momerr.Y()*momerr.Y()+
164  mom.Z()*mom.Z() * momerr.Z()*momerr.Z());
165  fCov[1][1] = pow((u.X()/pw - w.X()*pu/(pw*pw)),2.)*momerr.X()*momerr.X() +
166  pow((u.Y()/pw - w.Y()*pu/(pw*pw)),2.)*momerr.Y()*momerr.Y() +
167  pow((u.Z()/pw - w.Z()*pu/(pw*pw)),2.)*momerr.Z()*momerr.Z();
168  fCov[2][2] = pow((v.X()/pw - w.X()*pv/(pw*pw)),2.)*momerr.X()*momerr.X() +
169  pow((v.Y()/pw - w.Y()*pv/(pw*pw)),2.)*momerr.Y()*momerr.Y() +
170  pow((v.Z()/pw - w.Z()*pv/(pw*pw)),2.)*momerr.Z()*momerr.Z();
171 
172 }
TVector3 getV() const
Definition: GFDetPlane.h:77
void setO(const TVector3 &o)
Definition: GFDetPlane.cxx:92
Double_t mom
Definition: plot_dirc.C:14
__m128 v
Definition: P4_F32vec4.h:4
void setNormal(TVector3 n)
Definition: GFDetPlane.cxx:156
TVector3 getPosSeed() const
get the seed value for track: pos
Definition: GFTrackCand.h:131
TVector3 getDirSeed() const
get the seed value for track: direction
Definition: GFTrackCand.h:133
TVector3 getDirError() const
get the seed value for track: error on direction (standard deviation)
Definition: GFTrackCand.h:139
GFDetPlane fRefPlane
TVector3 getPosError() const
get the seed value for track: error on pos (standard deviation)
Definition: GFTrackCand.h:137
int getPdgCode() const
get the PDG code
Definition: GFTrackCand.h:141
TVector3 getU() const
Definition: GFDetPlane.h:76
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
void setPDG(int)
Set PDG particle code.
Definition: RKTrackRep.cxx:274
TVector3 getO() const
Definition: GFDetPlane.h:75
TMatrixT< double > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:91
RKTrackRep::RKTrackRep ( const TVector3 &  pos,
const TVector3 &  mom,
const int &  PDGCode 
)

Definition at line 174 of file RKTrackRep.cxx.

References fCharge, GFAbsTrackRep::fCov, GFAbsTrackRep::fRefPlane, fSpu, GFAbsTrackRep::fState, GFDetPlane::getO(), GFDetPlane::getU(), GFDetPlane::getV(), GFDetPlane::setNormal(), GFDetPlane::setO(), setPDG(), and v.

176  :
177  GFAbsTrackRep(5),fDirection(true), fCachePlane(), fCacheSpu(1), fAuxInfo(1,1){
178  setPDG(PDGCode); // also sets charge and mass
179 
180 
181  fRefPlane.setO(pos);
183 
184  fState[0][0]=fCharge/mom.Mag();
185  //u' and v'
186  fState[1][0]=0.0;
187  fState[2][0]=0.0;
188  //u and v
189  fState[3][0]=0.0;
190  fState[4][0]=0.0;
191  //spu
192  fSpu=1.;
193 
194  TVector3 o=fRefPlane.getO();
195  TVector3 u=fRefPlane.getU();
196  TVector3 v=fRefPlane.getV();
197  TVector3 w=u.Cross(v);
198  double pu=0.;
199  double pv=0.;
200  double pw=mom.Mag();
201 
202  static const TVector3 stdPosErr(1.,1.,1.);
203  static const TVector3 stdMomErr(1.,1.,1.);
204 
205  fCov[3][3] = stdPosErr.X()*stdPosErr.X() * u.X()*u.X() +
206  stdPosErr.Y()*stdPosErr.Y() * u.Y()*u.Y() +
207  stdPosErr.Z()*stdPosErr.Z() * u.Z()*u.Z();
208  fCov[4][4] = stdPosErr.X()*stdPosErr.X() * v.X()*v.X() +
209  stdPosErr.Y()*stdPosErr.Y() * v.Y()*v.Y() +
210  stdPosErr.Z()*stdPosErr.Z() * v.Z()*v.Z();
211  fCov[0][0] = fCharge*fCharge/pow(mom.Mag(),6.) *
212  (mom.X()*mom.X() * stdMomErr.X()*stdMomErr.X()+
213  mom.Y()*mom.Y() * stdMomErr.Y()*stdMomErr.Y()+
214  mom.Z()*mom.Z() * stdMomErr.Z()*stdMomErr.Z());
215  fCov[1][1] = pow((u.X()/pw - w.X()*pu/(pw*pw)),2.)*stdMomErr.X()*stdMomErr.X() +
216  pow((u.Y()/pw - w.Y()*pu/(pw*pw)),2.)*stdMomErr.Y()*stdMomErr.Y() +
217  pow((u.Z()/pw - w.Z()*pu/(pw*pw)),2.)*stdMomErr.Z()*stdMomErr.Z();
218  fCov[2][2] = pow((v.X()/pw - w.X()*pv/(pw*pw)),2.)*stdMomErr.X()*stdMomErr.X() +
219  pow((v.Y()/pw - w.Y()*pv/(pw*pw)),2.)*stdMomErr.Y()*stdMomErr.Y() +
220  pow((v.Z()/pw - w.Z()*pv/(pw*pw)),2.)*stdMomErr.Z()*stdMomErr.Z();
221 
222 }
TVector3 pos
TVector3 getV() const
Definition: GFDetPlane.h:77
void setO(const TVector3 &o)
Definition: GFDetPlane.cxx:92
Double_t mom
Definition: plot_dirc.C:14
__m128 v
Definition: P4_F32vec4.h:4
void setNormal(TVector3 n)
Definition: GFDetPlane.cxx:156
GFDetPlane fRefPlane
TVector3 getU() const
Definition: GFDetPlane.h:76
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
void setPDG(int)
Set PDG particle code.
Definition: RKTrackRep.cxx:274
TVector3 getO() const
Definition: GFDetPlane.h:75
TMatrixT< double > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:91
RKTrackRep::RKTrackRep ( const GFDetPlane pl,
const TVector3 &  mom,
const int &  PDGCode 
)

Definition at line 224 of file RKTrackRep.cxx.

References fCharge, GFAbsTrackRep::fCov, GFAbsTrackRep::fRefPlane, fSpu, GFAbsTrackRep::fState, GFDetPlane::getO(), GFDetPlane::getU(), GFDetPlane::getV(), setPDG(), and v.

226  :
227  GFAbsTrackRep(5),fDirection(true), fCachePlane(), fCacheSpu(1), fAuxInfo(1,1){
228  setPDG(PDGCode); // also sets charge and mass
229 
230 
231  fRefPlane = pl;
232  TVector3 o=fRefPlane.getO();
233  TVector3 u=fRefPlane.getU();
234  TVector3 v=fRefPlane.getV();
235  TVector3 w=u.Cross(v);
236 
237  double pu=mom*u;
238  double pv=mom*v;
239  double pw=mom*w;
240 
241  fState[0][0]=fCharge/mom.Mag();
242  //u' and v'
243  fState[1][0]=pu/pw;
244  fState[2][0]=pv/pw;
245  //u and v
246  fState[3][0]=0.0;
247  fState[4][0]=0.0;
248  //spu
249  fSpu=1.;
250 
251  static const TVector3 stdPosErr2(1.,1.,1.);
252  static const TVector3 stdMomErr2(1.,1.,1.);
253 
254  fCov[3][3] = stdPosErr2.X()*stdPosErr2.X() * u.X()*u.X() +
255  stdPosErr2.Y()*stdPosErr2.Y() * u.Y()*u.Y() +
256  stdPosErr2.Z()*stdPosErr2.Z() * u.Z()*u.Z();
257  fCov[4][4] = stdPosErr2.X()*stdPosErr2.X() * v.X()*v.X() +
258  stdPosErr2.Y()*stdPosErr2.Y() * v.Y()*v.Y() +
259  stdPosErr2.Z()*stdPosErr2.Z() * v.Z()*v.Z();
260  fCov[0][0] = fCharge*fCharge/pow(mom.Mag(),6.) *
261  (mom.X()*mom.X() * stdMomErr2.X()*stdMomErr2.X()+
262  mom.Y()*mom.Y() * stdMomErr2.Y()*stdMomErr2.Y()+
263  mom.Z()*mom.Z() * stdMomErr2.Z()*stdMomErr2.Z());
264  fCov[1][1] = pow((u.X()/pw - w.X()*pu/(pw*pw)),2.)*stdMomErr2.X()*stdMomErr2.X() +
265  pow((u.Y()/pw - w.Y()*pu/(pw*pw)),2.)*stdMomErr2.Y()*stdMomErr2.Y() +
266  pow((u.Z()/pw - w.Z()*pu/(pw*pw)),2.)*stdMomErr2.Z()*stdMomErr2.Z();
267  fCov[2][2] = pow((v.X()/pw - w.X()*pv/(pw*pw)),2.)*stdMomErr2.X()*stdMomErr2.X() +
268  pow((v.Y()/pw - w.Y()*pv/(pw*pw)),2.)*stdMomErr2.Y()*stdMomErr2.Y() +
269  pow((v.Z()/pw - w.Z()*pv/(pw*pw)),2.)*stdMomErr2.Z()*stdMomErr2.Z();
270 
271 }
TVector3 getV() const
Definition: GFDetPlane.h:77
Double_t mom
Definition: plot_dirc.C:14
__m128 v
Definition: P4_F32vec4.h:4
GFDetPlane fRefPlane
TVector3 getU() const
Definition: GFDetPlane.h:76
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
void setPDG(int)
Set PDG particle code.
Definition: RKTrackRep.cxx:274
TVector3 getO() const
Definition: GFDetPlane.h:75
TMatrixT< double > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:91
RKTrackRep::~RKTrackRep ( )
virtual

Definition at line 69 of file RKTrackRep.cxx.

69  {
70  //GFMaterialEffects is now a Singleton
71 }

Member Function Documentation

void GFAbsTrackRep::addChiSqu ( double  aChiSqu)
inlineinherited

Definition at line 303 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fChiSqu.

Referenced by GFKalman::processHit().

303  {
304  fChiSqu += aChiSqu;
305  }
double fChiSqu
chiSqu of the track fit
Definition: GFAbsTrackRep.h:94
void GFAbsTrackRep::addNDF ( unsigned int  n)
inlineinherited

Definition at line 306 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fNdf, and n.

Referenced by GFKalman::processHit().

306  {
307  fNdf += n;
308  }
int n
unsigned int fNdf
Definition: GFAbsTrackRep.h:95
virtual GFAbsTrackRep* RKTrackRep::clone ( ) const
inlinevirtual

Implements GFAbsTrackRep.

Definition at line 74 of file tracking/GenfitTools/trackrep/RKTrackRep/RKTrackRep.h.

References RKTrackRep().

74 {return new RKTrackRep(*this);}
double RKTrackRep::Extrap ( const GFDetPlane plane,
TMatrixT< double > *  state,
TMatrixT< double > *  cov = NULL 
) const
private

Handles propagation and material effects.

extrapolate(), extrapolateToPoint() and extrapolateToLine() call this function. Extrap() needs a plane as an argument, hence extrapolateToPoint() and extrapolateToLine() create virtual detector planes. In this function, RKutta() is called and the resulting points and point paths are filtered so that the direction doesn't change and tiny steps are filtered out. After the propagation the material effects in #fEffect are called. Extrap() will loop until the plane is reached, unless the propagation fails or the maximum number of iterations is exceeded.

Definition at line 991 of file RKTrackRep.cxx.

References GFDetPlane::dist(), GFDetPlane::distance(), GFMaterialEffects::effects(), fabs(), fCharge, fPdg, GFMaterialEffects::getInstance(), i, GFDetPlane::inActive(), MINSTEP, noise, P, RKutta(), and GFException::setFatal().

Referenced by extrapolate(), extrapolateToLine(), extrapolateToPoint(), and stepalong().

991  {
992 
993  static const int maxNumIt(2000);
994  int numIt(0);
995  bool calcCov(true);
996  if(cov==NULL) calcCov=false;
997 
998  double *P;
999  if(calcCov) {P = new double[56]; memset(P,0x00,56*sizeof(double));}
1000  else {P = new double[7];} // not needed memset(P,0x00,7*sizeof(double));};
1001 
1002  for(int i=0;i<7;++i){
1003  P[i] = (*state)[i][0];
1004  }
1005 
1006  TMatrixT<double> jac(7,7);
1007  TMatrixT<double> jacT(7,7);
1008  TMatrixT<double> oldCov(7,7);
1009  if(calcCov) oldCov=(*cov);
1010  double coveredDistance(0.);
1011  double sumDistance(0.);
1012 
1013  while(true){
1014  if(numIt++ > maxNumIt){
1015  GFException exc("RKTrackRep::Extrap ==> maximum number of iterations exceeded",
1016  __LINE__,__FILE__);
1017  exc.setFatal();
1018  delete[] P;
1019  throw exc;
1020  }
1021 
1022  if(calcCov){
1023  memset(&P[7],0x00,49*sizeof(double));
1024  for(int i=0; i<6; ++i){
1025  P[(i+1)*7+i] = 1.;
1026  }
1027  P[55] = (*state)[6][0];
1028  }
1029 
1030  //double dir(1.); //[R.K. 01/2017] unused variable?
1031  {
1032  TVector3 Pvect(P[0],P[1],P[2]); //position
1033  TVector3 Avect(P[3],P[4],P[5]); //direction
1034  TVector3 dist = plane.dist(Pvect); //from point to plane
1035  //if(dist*Avect<0.) dir=-1.; //[R.K. 01/2017] unused variable?
1036  }
1037 
1038  TVector3 directionBefore(P[3],P[4],P[5]); // direction before propagation
1039  directionBefore.SetMag(1.);
1040 
1041  // propagation
1042  std::vector<TVector3> points;
1043  std::vector<double> pointPaths;
1044  if( ! this->RKutta(plane, P, coveredDistance, points, pointPaths, -1., calcCov) ) { // maxLen currently not used
1045  GFException exc("RKTrackRep::Extrap ==> Runge Kutta propagation failed",__LINE__,__FILE__);
1046  exc.setFatal(); // stops propagation; faster, but some hits will be lost
1047  delete[] P;
1048  throw exc;
1049  }
1050 
1051  TVector3 directionAfter(P[3],P[4],P[5]); // direction after propagation
1052  directionAfter.SetMag(1.);
1053 
1054  sumDistance+=coveredDistance;
1055 
1056  // filter Points
1057  std::vector<TVector3> pointsFilt(1, points.at(0));
1058  std::vector<double> pointPathsFilt(1, 0.);
1059  // only if in right direction
1060  for(unsigned int i=1;i<points.size();++i){
1061  if (pointPaths.at(i) * coveredDistance > 0.) {
1062  pointsFilt.push_back(points.at(i));
1063  pointPathsFilt.push_back(pointPaths.at(i));
1064  }
1065  else {
1066  pointsFilt.back() = points.at(i);
1067  pointPathsFilt.back() += pointPaths.at(i);
1068  }
1069  // clean up tiny steps
1070  int position = pointsFilt.size()-1; // position starts with 0
1071  if (fabs(pointPathsFilt.back()) < MINSTEP && position > 1) {
1072  pointsFilt.at(position-1) = pointsFilt.at(position);
1073  pointsFilt.pop_back();
1074  pointPathsFilt.at(position-1) += pointPathsFilt.at(position);
1075  pointPathsFilt.pop_back();
1076  }
1077  }
1078 
1079  //consistency check
1080  double checkSum(0.);
1081  for(unsigned int i=0;i<pointPathsFilt.size();++i){
1082  checkSum+=pointPathsFilt.at(i);
1083  }
1084  if(fabs(checkSum-coveredDistance)>1.E-7){
1085  GFException exc("RKTrackRep::Extrap ==> fabs(checkSum-coveredDistance)>1.E-7",__LINE__,__FILE__);
1086  exc.setFatal();
1087  delete[] P;
1088  throw exc;
1089  }
1090 
1091  if(calcCov){ //calculate Jacobian jac
1092  for(int i=0;i<7;++i){
1093  for(int j=0;j<7;++j){
1094  if(i<6) jac[i][j] = P[ (i+1)*7+j ];
1095  else jac[i][j] = P[ (i+1)*7+j ]/P[6];
1096  }
1097  }
1098  jacT = jac;
1099  jacT.T();
1100  }
1101 
1102  TMatrixT<double> noise(7,7); // zero everywhere by default
1103 
1104  // call MatEffects
1105  double momLoss; // momLoss has a sign - negative loss means momentum gain
1106 
1107  momLoss = GFMaterialEffects::getInstance()->effects(pointsFilt,
1108  pointPathsFilt,
1109  fabs(fCharge/P[6]), // momentum
1110  fPdg,
1111  calcCov,
1112  &noise,
1113  &jac,
1114  &directionBefore,
1115  &directionAfter);
1116 
1117  if(fabs(P[6])>1.E-10){ // do momLoss only for defined 1/momentum .ne.0
1118  P[6] = fCharge/(fabs(fCharge/P[6])-momLoss);
1119  }
1120 
1121  if(calcCov){ //propagate cov and add noise
1122  if(!(oldCov < 1.E200)){
1123  GFException exc("RKTrackRep::Extrap ==> covariance matrix exceeds numerical limits",__LINE__,__FILE__);
1124  exc.setFatal();
1125  delete[] P;
1126  throw exc;
1127  }
1128  oldCov = *cov;
1129  *cov = jacT*((oldCov)*jac)+noise;
1130  }
1131 
1132 
1133  //we arrived at the destination plane, if we point to the active area
1134  //of the plane (if it is finite), and the distance is below threshold
1135  if( plane.inActive(TVector3(P[0],P[1],P[2]),TVector3(P[3],P[4],P[5]))) {
1136  if(plane.distance(P[0],P[1],P[2])<MINSTEP) break;
1137  }
1138  }
1139  (*state)[0][0] = P[0]; (*state)[1][0] = P[1];
1140  (*state)[2][0] = P[2]; (*state)[3][0] = P[3];
1141  (*state)[4][0] = P[4]; (*state)[5][0] = P[5];
1142  (*state)[6][0] = P[6];
1143 
1144  delete[] P;
1145  return sumDistance;
1146 }
Int_t i
Definition: run_full.C:25
double distance(TVector3 &) const
Definition: GFDetPlane.cxx:343
bool RKutta(const GFDetPlane &plane, double *P, double &coveredDistance, std::vector< TVector3 > &points, std::vector< double > &pointLengths, const double &maxLen=-1, bool calcCov=true) const
Propagates the particle through the magnetic field.
Definition: RKTrackRep.cxx:640
#define MINSTEP
Definition: RKTrackRep.cxx:40
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
GeV c P
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:45
bool inActive(const TVector3 &point, const TVector3 &dir) const
intersect in the active area? C.f. GFAbsFinitePlane
Definition: GFDetPlane.h:132
static GFMaterialEffects * getInstance()
double noise
double effects(const std::vector< TVector3 > &points, const std::vector< double > &pointPaths, const double &mom, const int &pdg, const bool &doNoise=false, TMatrixT< double > *noise=NULL, const TMatrixT< double > *jacobian=NULL, const TVector3 *directionBefore=NULL, const TVector3 *directionAfter=NULL)
Calculates energy loss in the travelled path, optional calculation of noise matrix.
TVector3 dist(const TVector3 &point) const
Definition: GFDetPlane.cxx:207
double RKTrackRep::extrapolate ( const GFDetPlane pl,
TMatrixT< double > &  statePred,
TMatrixT< double > &  covPred 
)
virtual

returns the tracklength spanned in this extrapolation

The covariance matrix is transformed from the plane coordinate system to the master reference system (for the propagation) and, after propagation, back to the plane coordinate system.
Also the parameter spu (which is +1 or -1 and indicates the direction of the particle) is calculated and stored in fCacheSpu. The plane is stored in fCachePlane.
Master reference system (MARS):

\begin{eqnarray*}x & = & O_{x}+uU_{x}+vV_{x}\\y & = & O_{y}+uU_{y}+vV_{y}\\z & = & O_{z}+uU_{z}+vV_{z}\\a_{x} & = & \frac{\mbox{spu}}{\widetilde{p}}\left(N_{x}+u\prime U_{x}+v\prime V_{x}\right)\\a_{y} & = & \frac{\mbox{spu}}{\widetilde{p}}\left(N_{y}+u\prime U_{y}+v\prime V_{y}\right)\\a_{z} & = & \frac{\mbox{spu}}{\widetilde{p}}\left(N_{z}+u\prime U_{z}+v\prime V_{z}\right)\\\frac{q}{p} & = & \frac{q}{p}\end{eqnarray*}

Plane coordinate system:

\begin{eqnarray*}u & = & \left(x-O_{x}\right)U_{x}+\left(y-O_{y}\right)U_{y}+\left(z-O_{z}\right)U_{z}\\v & = & \left(x-O_{x}\right)U_{x}+\left(y-O_{y}\right)U_{y}+\left(z-O_{z}\right)U_{z}\\u\prime & = & \frac{a_{x}U_{x}+a_{y}U_{y}+a_{z}U_{z}}{a_{x}N_{x}+a_{y}N_{y}+a_{z}N_{z}}\\v\prime & = & \frac{a_{x}V_{x}+a_{y}V_{y}+a_{z}V_{z}}{a_{x}N_{x}+a_{y}N_{y}+a_{z}N_{z}}\\\frac{q}{p} & = & \frac{q}{p}\\\mbox{spu} & = & \frac{a_{x}N_{x}+a_{y}N_{y}+a_{z}N_{z}}{|a_{x}^{2}N_{x}^{2}+a_{y}^{2}N_{y}^{2}+a_{z}^{2}N_{z}^{2}|}=\pm1\end{eqnarray*}

Jacobians:
\(J_{p,M}=\left(\begin{array}{ccccc}\frac{\partial x}{\partial u} & \frac{\partial x}{\partial v} & \frac{\partial x}{\partial u\prime} & \frac{\partial x}{\partial v\prime} & \frac{\partial x}{\partial\frac{q}{p}}\\\frac{\partial y}{\partial u} & \frac{\partial y}{\partial v} & \frac{\partial y}{\partial u\prime} & \frac{\partial y}{\partial v\prime} & \frac{\partial y}{\partial\frac{q}{p}}\\\frac{\partial z}{\partial u} & \frac{\partial z}{\partial v} & \frac{\partial z}{\partial u\prime} & \frac{\partial z}{\partial v\prime} & \frac{\partial z}{\partial\frac{q}{p}}\\\frac{\partial a_{x}}{\partial u} & \frac{\partial a_{x}}{\partial v} & \frac{\partial a_{x}}{\partial u\prime} & \frac{\partial a_{x}}{\partial v\prime} & \frac{\partial a_{x}}{\partial\frac{q}{p}}\\\frac{\partial a_{y}}{\partial u} & \frac{\partial a_{y}}{\partial v} & \frac{\partial a_{y}}{\partial u\prime} & \frac{\partial a_{y}}{\partial v\prime} & \frac{\partial a_{y}}{\partial\frac{q}{p}}\\\frac{\partial a_{z}}{\partial u} & \frac{\partial a_{z}}{\partial v} & \frac{\partial a_{z}}{\partial u\prime} & \frac{\partial a_{z}}{\partial v\prime} & \frac{\partial a_{z}}{\partial\frac{q}{p}}\\\frac{\partial\frac{q}{p}}{\partial u} & \frac{\partial\frac{q}{p}}{\partial v} & \frac{\partial\frac{q}{p}}{\partial u\prime} & \frac{\partial\frac{q}{p}}{\partial v\prime} & \frac{\partial\frac{q}{p}}{\partial\frac{q}{p}}\end{array}\right)\)

\(J_{p,M}=\left(\begin{array}{cccccc}U_{x} & V_{x} & 0 & 0 & 0\\U_{y} & V_{y} & 0 & 0 & 0\\U_{z} & V_{z} & 0 & 0 & 0\\0 & 0 & \left\{ \frac{\textrm{spu}}{\widetilde{p}}U_{x}-\left[\frac{\textrm{spu}}{\widetilde{p}^{3}}\widetilde{p}_{x}\left(U_{x}\widetilde{p}_{x}+U_{y}\widetilde{p}_{y}+U_{z}\widetilde{p}_{z}\right)\right]\right\} & \left\{ \frac{\textrm{spu}}{\widetilde{p}}V_{x}-\left[\frac{\textrm{spu}}{\widetilde{p}^{3}}\widetilde{p}_{x}\left(V_{x}\widetilde{p}_{x}+V_{y}\widetilde{p}_{y}+V_{z}\widetilde{p}_{z}\right)\right]\right\} & 0\\0 & 0 & \left\{ \frac{\textrm{spu}}{\widetilde{p}}U_{y}-\left[\frac{\textrm{spu}}{\widetilde{p}^{3}}\widetilde{p}_{y}\left(U_{x}\widetilde{p}_{x}+U_{y}\widetilde{p}_{y}+U_{z}\widetilde{p}_{z}\right)\right]\right\} & \left\{ \frac{\textrm{spu}}{\widetilde{p}}V_{y}-\left[\frac{\textrm{spu}}{\widetilde{p}^{3}}\widetilde{p}_{y}\left(V_{x}\widetilde{p}_{x}+V_{y}\widetilde{p}_{y}+V_{z}\widetilde{p}_{z}\right)\right]\right\} & 0\\0 & 0 & \left\{ \frac{\textrm{spu}}{\widetilde{p}}U_{z}-\left[\frac{\textrm{spu}}{\widetilde{p}^{3}}\widetilde{p}_{z}\left(U_{x}\widetilde{p}_{x}+U_{y}\widetilde{p}_{y}+U_{z}\widetilde{p}_{z}\right)\right]\right\} & \left\{ \frac{\textrm{spu}}{\widetilde{p}}V_{z}-\left[\frac{\textrm{spu}}{\widetilde{p}^{3}}\widetilde{p}_{z}\left(V_{x}\widetilde{p}_{x}+V_{y}\widetilde{p}_{y}+V_{z}\widetilde{p}_{z}\right)\right]\right\} & 0\\0 & 0 & 0 & 0 & 1\end{array}\right)\)

with

\begin{eqnarray*}\widetilde{p} & = & \sqrt{\widetilde{p}_{x}^{2}+\widetilde{p}_{y}^{2}+\widetilde{p}_{z}^{2}}\\\widetilde{p}_{x} & = & \textrm{spu}\left(N_{x}+u\prime U_{x}+v\prime V_{x}\right)\\\widetilde{p}_{y} & = & \textrm{spu}\left(N_{y}+u\prime U_{y}+v\prime V_{y}\right)\\\widetilde{p}_{z} & = & \textrm{spu}\left(N_{z}+u\prime U_{z}+v\prime V_{z}\right)\end{eqnarray*}

\(J_{M,p}=\left(\begin{array}{ccccccc}\frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} & \frac{\partial u}{\partial z} & \frac{\partial u}{\partial a_{x}} & \frac{\partial u}{\partial a_{y}} & \frac{\partial u}{\partial a_{z}} & \frac{\partial u}{\partial\frac{q}{p}}\\\frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} & \frac{\partial v}{\partial z} & \frac{\partial v}{\partial a_{x}} & \frac{\partial v}{\partial a_{y}} & \frac{\partial v}{\partial a_{z}} & \frac{\partial v}{\partial\frac{q}{p}}\\\frac{\partial u\prime}{\partial x} & \frac{\partial u\prime}{\partial y} & \frac{\partial u\prime}{\partial z} & \frac{\partial u\prime}{\partial a_{x}} & \frac{\partial u\prime}{\partial a_{y}} & \frac{\partial u\prime}{\partial a_{z}} & \frac{\partial u\prime}{\partial\frac{q}{p}}\\\frac{\partial v\prime}{\partial x} & \frac{\partial v\prime}{\partial y} & \frac{\partial v\prime}{\partial z} & \frac{\partial v\prime}{\partial a_{x}} & \frac{\partial v\prime}{\partial a_{y}} & \frac{\partial v\prime}{\partial a_{z}} & \frac{\partial v\prime}{\partial\frac{q}{p}}\\ \frac{\partial\frac{q}{p}}{\partial x} & \frac{\partial\frac{q}{p}}{\partial y} & \frac{\partial\frac{q}{p}}{\partial z} & \frac{\partial\frac{q}{p}}{\partial a_{x}} & \frac{\partial\frac{q}{p}}{\partial a_{y}} & \frac{\partial\frac{q}{p}}{\partial a_{z}} & \frac{\partial\frac{q}{p}}{\partial\frac{q}{p}}\\ \\\end{array}\right)\)

\(J_{M,p}=\left(\begin{array}{ccccccc} U_{x} & U_{y} & U_{z} & 0 & 0 & 0 & 0\\ V_{x} & V_{y} & V_{z} & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & \frac{U_{x}\left(a_{y}N_{y}+a_{z}N_{z}\right)-N_{x}\left(a_{y}U_{y}+a_{z}U_{z}\right)}{\left(a_{x}N_{x}+a_{y}N_{y}+a_{z}N_{z}\right)^{2}} & \frac{U_{y}\left(a_{x}N_{x}+a_{z}N_{z}\right)-N_{y}\left(a_{x}U_{x}+a_{z}U_{z}\right)}{\left(a_{x}N_{x}+a_{y}N_{y}+a_{z}N_{z}\right)^{2}} & \frac{U_{z}\left(a_{x}N_{x}+a_{y}N_{y}\right)-N_{z}\left(a_{x}U_{x}+a_{y}U_{y}\right)}{\left(a_{x}N_{x}+a_{y}N_{y}+a_{z}N_{z}\right)^{2}} & 0\\ 0 & 0 & 0 & \frac{V_{x}\left(a_{y}N_{y}+a_{z}N_{z}\right)-N_{x}\left(a_{y}V_{y}+a_{z}V_{z}\right)}{\left(a_{x}N_{x}+a_{y}N_{y}+a_{z}N_{z}\right)^{2}} & \frac{V_{y}\left(a_{x}N_{x}+a_{z}N_{z}\right)-N_{y}\left(a_{x}V_{x}+a_{z}V_{z}\right)}{\left(a_{x}N_{x}+a_{y}N_{y}+a_{z}N_{z}\right)^{2}} & \frac{V_{z}\left(a_{x}N_{x}+a_{y}N_{y}\right)-N_{z}\left(a_{x}V_{x}+a_{y}V_{y}\right)}{\left(a_{x}N_{x}+a_{y}N_{y}+a_{z}N_{z}\right)^{2}} & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 1\\ \\\end{array}\right)\)

Implements GFAbsTrackRep.

Definition at line 442 of file RKTrackRep.cxx.

References Extrap(), fabs(), fCachePlane, fCacheSpu, GFAbsTrackRep::fCov, GFAbsTrackRep::fRefPlane, fSpu, GFAbsTrackRep::fState, GFDetPlane::getNormal(), GFDetPlane::getO(), GFDetPlane::getU(), GFDetPlane::getV(), pos, v, W, X, Y, and Z.

Referenced by getMom(), getPos(), getPosMom(), and getPosMomCov().

444  {
445 
446  TMatrixT<double> cov7x7(7,7);
447  TMatrixT<double> J_pM(7,5);
448 
449  TVector3 o=fRefPlane.getO();
450  TVector3 u=fRefPlane.getU();
451  TVector3 v=fRefPlane.getV();
452  TVector3 w=u.Cross(v);
453 
454  J_pM[0][3] = u.X();J_pM[0][4]=v.X(); // dx/du
455  J_pM[1][3] = u.Y();J_pM[1][4]=v.Y();
456  J_pM[2][3] = u.Z();J_pM[2][4]=v.Z();
457 
458  TVector3 pTilde = fSpu * (w + fState[1][0] * u + fState[2][0] * v);
459  double pTildeMag = pTilde.Mag();
460 
461  //J_pM matrix is d(x,y,z,ax,ay,az,q/p) / d(q/p,u',v',u,v)
462 
463  // da_x/du'
464  J_pM[3][1] = fSpu/pTildeMag*(u.X()-pTilde.X()/(pTildeMag*pTildeMag)*u*pTilde);
465  J_pM[4][1] = fSpu/pTildeMag*(u.Y()-pTilde.Y()/(pTildeMag*pTildeMag)*u*pTilde);
466  J_pM[5][1] = fSpu/pTildeMag*(u.Z()-pTilde.Z()/(pTildeMag*pTildeMag)*u*pTilde);
467  // da_x/dv'
468  J_pM[3][2] = fSpu/pTildeMag*(v.X()-pTilde.X()/(pTildeMag*pTildeMag)*v*pTilde);
469  J_pM[4][2] = fSpu/pTildeMag*(v.Y()-pTilde.Y()/(pTildeMag*pTildeMag)*v*pTilde);
470  J_pM[5][2] = fSpu/pTildeMag*(v.Z()-pTilde.Z()/(pTildeMag*pTildeMag)*v*pTilde);
471  // dqOp/dqOp
472  J_pM[6][0] = 1.;
473 
474  TMatrixT<double> J_pM_transp(J_pM);
475  J_pM_transp.T();
476 
477  cov7x7 = J_pM*(fCov*J_pM_transp);
478 
479  TVector3 pos = o + fState[3][0]*u + fState[4][0]*v;
480  TMatrixT<double> state7(7,1);
481  state7[0][0] = pos.X();
482  state7[1][0] = pos.Y();
483  state7[2][0] = pos.Z();
484  state7[3][0] = pTilde.X()/pTildeMag;;
485  state7[4][0] = pTilde.Y()/pTildeMag;;
486  state7[5][0] = pTilde.Z()/pTildeMag;;
487  state7[6][0] = fState[0][0];
488 
489  double coveredDistance = this->Extrap(pl,&state7,&cov7x7);
490 
491  TVector3 O = pl.getO();
492  TVector3 U = pl.getU();
493  TVector3 V = pl.getV();
494  TVector3 W = pl.getNormal();
495 
496  double X = state7[0][0];
497  double Y = state7[1][0];
498  double Z = state7[2][0];
499  double AX = state7[3][0];
500  double AY = state7[4][0];
501  double AZ = state7[5][0];
502  double QOP = state7[6][0];
503  TVector3 A(AX,AY,AZ);
504  TVector3 Point(X,Y,Z);
505  TMatrixT<double> J_Mp(5,7);
506 
507  // J_Mp matrix is d(q/p,u',v',u,v) / d(x,y,z,ax,ay,az,q/p)
508  J_Mp[0][6] = 1.;
509  //du'/da_x
510  double AtW = A*W;
511  J_Mp[1][3] = (U.X()*(AtW)-W.X()*(A*U))/(AtW*AtW);
512  J_Mp[1][4] = (U.Y()*(AtW)-W.Y()*(A*U))/(AtW*AtW);
513  J_Mp[1][5] = (U.Z()*(AtW)-W.Z()*(A*U))/(AtW*AtW);
514  //dv'/da_x
515  J_Mp[2][3] = (V.X()*(AtW)-W.X()*(A*V))/(AtW*AtW);
516  J_Mp[2][4] = (V.Y()*(AtW)-W.Y()*(A*V))/(AtW*AtW);
517  J_Mp[2][5] = (V.Z()*(AtW)-W.Z()*(A*V))/(AtW*AtW);
518  //du/dx
519  J_Mp[3][0] = U.X();
520  J_Mp[3][1] = U.Y();
521  J_Mp[3][2] = U.Z();
522  //dv/dx
523  J_Mp[4][0] = V.X();
524  J_Mp[4][1] = V.Y();
525  J_Mp[4][2] = V.Z();
526 
527  TMatrixT<double> J_Mp_transp(J_Mp);
528  J_Mp_transp.T();
529 
530  covPred.ResizeTo(5,5);
531 
532  covPred = J_Mp*(cov7x7*J_Mp_transp);
533 
534  statePred.ResizeTo(5,1);
535  statePred[0][0] = QOP;
536  statePred[1][0] = (A*U)/(A*W);
537  statePred[2][0] = (A*V)/(A*W);
538  statePred[3][0] = (Point-O)*U;
539  statePred[4][0] = (Point-O)*V;
540 
541  fCachePlane = pl;
542  fCacheSpu = (A*W)/fabs(A*W);
543 
544  return coveredDistance;
545 }
TVector3 pos
TVector3 getV() const
Definition: GFDetPlane.h:77
double Extrap(const GFDetPlane &plane, TMatrixT< double > *state, TMatrixT< double > *cov=NULL) const
Handles propagation and material effects.
Definition: RKTrackRep.cxx:991
double Y
Definition: anaLmdDigi.C:68
__m128 v
Definition: P4_F32vec4.h:4
#define W
Definition: createSTT.C:76
GFDetPlane fRefPlane
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
double X
Definition: anaLmdDigi.C:68
TVector3 getU() const
Definition: GFDetPlane.h:76
double Z
Definition: anaLmdDigi.C:68
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
TVector3 getNormal() const
Definition: GFDetPlane.cxx:138
TVector3 getO() const
Definition: GFDetPlane.h:75
TMatrixT< double > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:91
double RKTrackRep::extrapolate ( const GFDetPlane pl,
TMatrixT< double > &  statePred 
)
virtual

returns the tracklength spanned in this extrapolation

Reimplemented from GFAbsTrackRep.

Definition at line 550 of file RKTrackRep.cxx.

References Extrap(), GFAbsTrackRep::fRefPlane, fSpu, GFAbsTrackRep::fState, GFDetPlane::getNormal(), GFDetPlane::getO(), GFDetPlane::getU(), GFDetPlane::getV(), pos, v, W, X, Y, and Z.

551  {
552 
553  TVector3 o=fRefPlane.getO();
554  TVector3 u=fRefPlane.getU();
555  TVector3 v=fRefPlane.getV();
556  TVector3 w=u.Cross(v);
557 
558 
559  TVector3 pTilde = fSpu* (w + fState[1][0] * u + fState[2][0] * v);
560  double pTildeMag = pTilde.Mag();
561 
562  TVector3 pos = o + fState[3][0]*u + fState[4][0]*v;
563 
564  TMatrixT<double> state7(7,1);
565  state7[0][0] = pos.X();
566  state7[1][0] = pos.Y();
567  state7[2][0] = pos.Z();
568  state7[3][0] = pTilde.X()/pTildeMag;
569  state7[4][0] = pTilde.Y()/pTildeMag;
570  state7[5][0] = pTilde.Z()/pTildeMag;
571  state7[6][0] = fState[0][0];
572 
573  TVector3 O = pl.getO();
574  TVector3 U = pl.getU();
575  TVector3 V = pl.getV();
576  TVector3 W = pl.getNormal();
577 
578  double coveredDistance = this->Extrap(pl,&state7);
579 
580  double X = state7[0][0];
581  double Y = state7[1][0];
582  double Z = state7[2][0];
583  double AX = state7[3][0];
584  double AY = state7[4][0];
585  double AZ = state7[5][0];
586  double QOP = state7[6][0];
587  TVector3 A(AX,AY,AZ);
588  TVector3 Point(X,Y,Z);
589 
590  statePred.ResizeTo(5,1);
591  statePred[0][0] = QOP;
592  statePred[1][0] = (A*U)/(A*W);
593  statePred[2][0] = (A*V)/(A*W);
594  statePred[3][0] = (Point-O)*U;
595  statePred[4][0] = (Point-O)*V;
596 
597  return coveredDistance;
598 }
TVector3 pos
TVector3 getV() const
Definition: GFDetPlane.h:77
double Extrap(const GFDetPlane &plane, TMatrixT< double > *state, TMatrixT< double > *cov=NULL) const
Handles propagation and material effects.
Definition: RKTrackRep.cxx:991
double Y
Definition: anaLmdDigi.C:68
__m128 v
Definition: P4_F32vec4.h:4
#define W
Definition: createSTT.C:76
GFDetPlane fRefPlane
double X
Definition: anaLmdDigi.C:68
TVector3 getU() const
Definition: GFDetPlane.h:76
double Z
Definition: anaLmdDigi.C:68
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
TVector3 getNormal() const
Definition: GFDetPlane.cxx:138
TVector3 getO() const
Definition: GFDetPlane.h:75
double GFAbsTrackRep::extrapolate ( const GFDetPlane plane)
inherited

This changes the state and cov and plane of the rep.

This method extrapolates to to the plane and sets the results of state, cov and also plane in itself.

Definition at line 33 of file GFAbsTrackRep.cxx.

References GFAbsTrackRep::extrapolate(), GFAbsTrackRep::fDimension, and GFAbsTrackRep::setData().

33  {
34  TMatrixT<double> statePred(fDimension,1);
35  TMatrixT<double> covPred(fDimension,fDimension);
36  double retVal = extrapolate(plane,statePred,covPred);
37  setData(statePred,plane,&covPred);
38  return retVal;
39 }
unsigned int fDimension
Dimensionality of track representation.
Definition: GFAbsTrackRep.h:85
virtual double extrapolate(const GFDetPlane &plane, TMatrixT< double > &statePred)
returns the tracklength spanned in this extrapolation
virtual void setData(const TMatrixT< double > &st, const GFDetPlane &pl, const TMatrixT< double > *cov=NULL, const TMatrixT< double > *aux=NULL)
Puts the track representation in a given state.
void RKTrackRep::extrapolateToLine ( const TVector3 &  point1,
const TVector3 &  point2,
TVector3 &  poca,
TVector3 &  dirInPoca,
TVector3 &  poca_onwire 
)
virtual

This method extrapolates to the point of closest approach to a line.

Reimplemented from GFAbsTrackRep.

Definition at line 390 of file RKTrackRep.cxx.

References Extrap(), fabs(), GFAbsTrackRep::fRefPlane, fSpu, GFAbsTrackRep::fState, GFDetPlane::getO(), GFDetPlane::getU(), GFDetPlane::getV(), MINSTEP, poca2Line(), point, GFDetPlane::setO(), GFDetPlane::setU(), GFDetPlane::setV(), and v.

394  {
395  static const int maxIt(30);
396 
397  TVector3 o=fRefPlane.getO();
398  TVector3 u=fRefPlane.getU();
399  TVector3 v=fRefPlane.getV();
400  TVector3 w=u.Cross(v);
401 
402  TVector3 pTilde = fSpu* (w + fState[1][0] * u + fState[2][0] * v);
403  pTilde.SetMag(1.);
404 
405  TVector3 point = o + fState[3][0]*u + fState[4][0]*v;
406 
407  TMatrixT<double> state7(7,1);
408  state7[0][0] = point.X();
409  state7[1][0] = point.Y();
410  state7[2][0] = point.Z();
411  state7[3][0] = pTilde.X();
412  state7[4][0] = pTilde.Y();
413  state7[5][0] = pTilde.Z();
414  state7[6][0] = fState[0][0];
415 
416  double coveredDistance(0.);
417 
418  GFDetPlane pl;
419  int iterations(0);
420 
421  while(true){
422  pl.setO(point1);
423  TVector3 currentDir(state7[3][0],state7[4][0],state7[5][0]);
424  pl.setU(currentDir.Cross(point2-point1));
425  pl.setV(point2-point1);
426  coveredDistance = this->Extrap(pl,&state7);
427 
428  if(fabs(coveredDistance)<MINSTEP) break;
429  if(++iterations == maxIt) {
430  GFException exc("RKTrackRep::extrapolateToLine ==> extrapolation to line failed, maximum number of iterations reached",__LINE__,__FILE__);
431  throw exc;
432  }
433  }
434  poca.SetXYZ(state7[0][0],state7[1][0],state7[2][0]);
435  dirInPoca.SetXYZ(state7[3][0],state7[4][0],state7[5][0]);
436  poca_onwire = poca2Line(point1,point2,poca);
437 }
TVector3 getV() const
Definition: GFDetPlane.h:77
double Extrap(const GFDetPlane &plane, TMatrixT< double > *state, TMatrixT< double > *cov=NULL) const
Handles propagation and material effects.
Definition: RKTrackRep.cxx:991
Detector plane genfit geometry class.
Definition: GFDetPlane.h:59
void setU(const TVector3 &u)
Definition: GFDetPlane.cxx:105
void setO(const TVector3 &o)
Definition: GFDetPlane.cxx:92
__m128 v
Definition: P4_F32vec4.h:4
GFDetPlane fRefPlane
#define MINSTEP
Definition: RKTrackRep.cxx:40
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
void setV(const TVector3 &v)
Definition: GFDetPlane.cxx:118
TVector3 poca2Line(const TVector3 &extr1, const TVector3 &extr2, const TVector3 &point) const
Definition: RKTrackRep.cxx:376
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:45
TVector3 getU() const
Definition: GFDetPlane.h:76
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
TVector3 getO() const
Definition: GFDetPlane.h:75
PndSdsMCPoint * point
Definition: anaLmdCluster.C:72
void RKTrackRep::extrapolateToPoint ( const TVector3 &  pos,
TVector3 &  poca,
TVector3 &  dirInPoca 
)
virtual

This method is to extrapolate the track to point of closest approach to a point in space.

Reimplemented from GFAbsTrackRep.

Definition at line 329 of file RKTrackRep.cxx.

References Extrap(), fabs(), GFAbsTrackRep::fRefPlane, fSpu, GFAbsTrackRep::fState, GFDetPlane::getO(), GFDetPlane::getU(), GFDetPlane::getV(), MINSTEP, point, GFDetPlane::setON(), and v.

Referenced by PndLmdBPRungeKuttaTask::Exec(), and PndLmdPerformanceTask::Exec().

331  {
332 
333  static const int maxIt(30);
334 
335  TVector3 o=fRefPlane.getO();
336  TVector3 u=fRefPlane.getU();
337  TVector3 v=fRefPlane.getV();
338  TVector3 w=u.Cross(v);
339 
340  TVector3 pTilde = fSpu* (w + fState[1][0] * u + fState[2][0] * v);
341  pTilde.SetMag(1.);
342 
343  TVector3 point = o + fState[3][0]*u + fState[4][0]*v;
344 
345  TMatrixT<double> state7(7,1);
346  state7[0][0] = point.X();
347  state7[1][0] = point.Y();
348  state7[2][0] = point.Z();
349  state7[3][0] = pTilde.X();
350  state7[4][0] = pTilde.Y();
351  state7[5][0] = pTilde.Z();
352  state7[6][0] = fState[0][0];
353 
354  double coveredDistance(0.);
355 
356  GFDetPlane pl;
357  int iterations(0);
358 
359  while(true){
360  pl.setON(pos,TVector3(state7[3][0],state7[4][0],state7[5][0]));
361  coveredDistance = this->Extrap(pl,&state7);
362 
363  if(fabs(coveredDistance)<MINSTEP) break;
364  if(++iterations == maxIt) {
365  GFException exc("RKTrackRep::extrapolateToPoint ==> extrapolation to point failed, maximum number of iterations reached",__LINE__,__FILE__);
366  throw exc;
367  }
368  }
369  poca.SetXYZ(state7[0][0],state7[1][0],state7[2][0]);
370  dirInPoca.SetXYZ(state7[3][0],state7[4][0],state7[5][0]);
371 }
TVector3 pos
TVector3 getV() const
Definition: GFDetPlane.h:77
double Extrap(const GFDetPlane &plane, TMatrixT< double > *state, TMatrixT< double > *cov=NULL) const
Handles propagation and material effects.
Definition: RKTrackRep.cxx:991
Detector plane genfit geometry class.
Definition: GFDetPlane.h:59
void setON(const TVector3 &o, const TVector3 &n)
Definition: GFDetPlane.cxx:145
__m128 v
Definition: P4_F32vec4.h:4
GFDetPlane fRefPlane
#define MINSTEP
Definition: RKTrackRep.cxx:40
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:45
TVector3 getU() const
Definition: GFDetPlane.h:76
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
TVector3 getO() const
Definition: GFDetPlane.h:75
PndSdsMCPoint * point
Definition: anaLmdCluster.C:72
const TMatrixT< double > * RKTrackRep::getAuxInfo ( const GFDetPlane )
virtual

Get auxillary information from the track representation.

AuxInfo is a mechanism which allows creators of track repersentations to hand out any information they like (as long as it is compatible with a TMatrixT<double>). It should be used if setData requires additional information to update the representation, but it can also be used for debugging information if needed. See also the documentation of GFAbsTrackRep::setData().

Reimplemented from GFAbsTrackRep.

Definition at line 57 of file RKTrackRep.cxx.

References fAuxInfo, fCachePlane, and fCacheSpu.

57  {
58 
59  if(pl!=fCachePlane) {
60  std::cerr << "RKTrackRep::getAuxInfo() - Fatal error: Trying to get auxillary information with planes mismatch (Information returned does not belong to requested plane)! -> abort in line " << __LINE__ << std::endl;
61  throw;
62  }
63  fAuxInfo.ResizeTo(1,1);
64  fAuxInfo(0,0) = fCacheSpu;
65  return &fAuxInfo;
66 
67 }
double RKTrackRep::getCharge ( ) const
inlinevirtual

Returns charge.

Implements GFAbsTrackRep.

Definition at line 144 of file tracking/GenfitTools/trackrep/RKTrackRep/RKTrackRep.h.

References fCharge.

double GFAbsTrackRep::getChiSqu ( ) const
inlineinherited

Definition at line 244 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fChiSqu.

Referenced by GFTrack::getChiSqu(), and GFAbsTrackRep::getRedChiSqu().

244  {
245  return fChiSqu;
246  }
double fChiSqu
chiSqu of the track fit
Definition: GFAbsTrackRep.h:94
TMatrixT<double> GFAbsTrackRep::getCov ( ) const
inlineinherited

Definition at line 198 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fCov.

Referenced by GFTrack::blowUpCovs(), GFKalman::processHit(), and GFKalman::processTrack().

198  {
199  return fCov;
200  }
TMatrixT< double > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:91
double GFAbsTrackRep::getCovElem ( int  i,
int  j 
) const
inlineinherited

Definition at line 203 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fCov.

203 {return fCov(i,j);}
Int_t i
Definition: run_full.C:25
TMatrixT< double > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:91
unsigned int GFAbsTrackRep::getDim ( ) const
inlineinherited

returns dimension of state vector

Definition at line 191 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fDimension.

Referenced by GeaneTrackRep::extrapolateToPoint(), GFKalman::getChi2Hit(), GFAbsTrackRep::getNDF(), GFTrack::getResiduals(), and GFKalman::processHit().

191 {return fDimension;}
unsigned int fDimension
Dimensionality of track representation.
Definition: GFAbsTrackRep.h:85
TMatrixT<double> GFAbsTrackRep::getFirstCov ( ) const
inlineinherited

Definition at line 229 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fFirstCov.

Referenced by GenfitTrack2PndTrack().

229  {
230  return fFirstCov;
231  }
TMatrixT< double > fFirstCov
GFDetPlane GFAbsTrackRep::getFirstPlane ( ) const
inlineinherited

Definition at line 232 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fFirstPlane.

Referenced by GenfitTrack2PndTrack().

232  {
233  return fFirstPlane;
234  }
GFDetPlane fFirstPlane
TMatrixT<double> GFAbsTrackRep::getFirstState ( ) const
inlineinherited

Definition at line 226 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fFirstState.

Referenced by GenfitTrack2PndTrack().

226  {
227  return fFirstState;
228  }
TMatrixT< double > fFirstState
state, cov and plane for first and last point in fit
TMatrixT<double> GFAbsTrackRep::getLastCov ( ) const
inlineinherited

Definition at line 238 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fLastCov.

Referenced by GenfitTrack2PndTrack().

238  {
239  return fLastCov;
240  }
TMatrixT< double > fLastCov
GFDetPlane GFAbsTrackRep::getLastPlane ( ) const
inlineinherited

Definition at line 241 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fLastPlane.

Referenced by GenfitTrack2PndTrack().

241  {
242  return fLastPlane;
243  }
GFDetPlane fLastPlane
TMatrixT<double> GFAbsTrackRep::getLastState ( ) const
inlineinherited

Definition at line 235 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fLastState.

Referenced by GenfitTrack2PndTrack().

235  {
236  return fLastState;
237  }
TMatrixT< double > fLastState
TVector3 RKTrackRep::getMom ( const GFDetPlane pl)
virtual

Returns momentum of the track in the plane.

If GFDetPlane equals the reference plane fRefPlane, returns current momentum; otherwise it extrapolates the track to the plane and returns the momentum.

Implements GFAbsTrackRep.

Definition at line 297 of file RKTrackRep.cxx.

References extrapolate(), fabs(), fCacheSpu, GFAbsTrackRep::fRefPlane, fSpu, GFAbsTrackRep::fState, GFDetPlane::getNormal(), GFDetPlane::getU(), and GFDetPlane::getV().

297  {
298  TMatrixT<double> statePred(fState);
299  TVector3 retmom;
300  if(pl!=fRefPlane) {
301  extrapolate(pl,statePred);
302  retmom = fCacheSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
303  }
304  else{
305  retmom = fSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
306  }
307  retmom.SetMag(1./fabs(statePred[0][0]));
308  return retmom;
309 }
TVector3 getV() const
Definition: GFDetPlane.h:77
GFDetPlane fRefPlane
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
TVector3 getU() const
Definition: GFDetPlane.h:76
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
TVector3 getNormal() const
Definition: GFDetPlane.cxx:138
double extrapolate(const GFDetPlane &, TMatrixT< double > &statePred, TMatrixT< double > &covPred)
returns the tracklength spanned in this extrapolation
Definition: RKTrackRep.cxx:442
TVector3 GFAbsTrackRep::getMom ( )
inlineinherited
unsigned int GFAbsTrackRep::getNDF ( ) const
inlineinherited

Definition at line 252 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fNdf, and GFAbsTrackRep::getDim().

Referenced by GFTrack::getNDF(), and GFAbsTrackRep::getRedChiSqu().

252  {
253  if(fNdf>getDim()) return fNdf-getDim();
254  return 0;
255  }
unsigned int fNdf
Definition: GFAbsTrackRep.h:95
unsigned int getDim() const
returns dimension of state vector
int RKTrackRep::getPDG ( )
inline

Definition at line 146 of file tracking/GenfitTools/trackrep/RKTrackRep/RKTrackRep.h.

References fPdg.

146 {return fPdg;};
TVector3 RKTrackRep::getPos ( const GFDetPlane pl)
virtual

Returns position of the track in the plane.

If GFDetPlane equals the reference plane fRefPlane, returns current position; otherwise it extrapolates the track to the plane and returns the position.

Implements GFAbsTrackRep.

Definition at line 287 of file RKTrackRep.cxx.

References extrapolate(), GFAbsTrackRep::fRefPlane, GFAbsTrackRep::fState, GFDetPlane::getO(), GFDetPlane::getU(), GFDetPlane::getV(), and s.

287  {
288  if(pl!=fRefPlane){
289  TMatrixT<double> s(5,1);
290  extrapolate(pl,s);
291  return pl.getO()+s[3][0]*pl.getU()+s[4][0]*pl.getV();
292  }
293  return fRefPlane.getO()+fState[3][0]*fRefPlane.getU()+fState[4][0]*fRefPlane.getV();
294 }
TVector3 getV() const
Definition: GFDetPlane.h:77
TLorentzVector s
Definition: Pnd2DStar.C:50
GFDetPlane fRefPlane
TVector3 getU() const
Definition: GFDetPlane.h:76
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
double extrapolate(const GFDetPlane &, TMatrixT< double > &statePred, TMatrixT< double > &covPred)
returns the tracklength spanned in this extrapolation
Definition: RKTrackRep.cxx:442
TVector3 getO() const
Definition: GFDetPlane.h:75
TVector3 GFAbsTrackRep::getPos ( )
inlineinherited

Definition at line 220 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fRefPlane, and GFAbsTrackRep::getPos().

Referenced by GeaneTrackRep::extrapolateToLine(), and GFAbsTrackRep::getPos().

220 {return getPos(fRefPlane);}
GFDetPlane fRefPlane
TVector3 getPos()
void RKTrackRep::getPosMom ( const GFDetPlane pl,
TVector3 &  pos,
TVector3 &  mom 
)
virtual

Gets position and momentum in the plane.

If GFDetPlane equals the reference plane fRefPlane, it gets current position and momentum; otherwise it extrapolates the track to the plane and gets the position and momentum.

Implements GFAbsTrackRep.

Definition at line 312 of file RKTrackRep.cxx.

References extrapolate(), fabs(), fCacheSpu, GFAbsTrackRep::fRefPlane, fSpu, GFAbsTrackRep::fState, GFDetPlane::getNormal(), GFDetPlane::getO(), GFDetPlane::getU(), and GFDetPlane::getV().

313  {
314  TMatrixT<double> statePred(fState);
315  if(pl!=fRefPlane) {
316  extrapolate(pl,statePred);
317  mom = fCacheSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
318  }
319  else{
320  mom = fSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
321  }
322  mom.SetMag(1./fabs(statePred[0][0]));
323  pos = pl.getO()+(statePred[3][0]*pl.getU())+(statePred[4][0]*pl.getV());
324 }
TVector3 pos
TVector3 getV() const
Definition: GFDetPlane.h:77
Double_t mom
Definition: plot_dirc.C:14
GFDetPlane fRefPlane
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
TVector3 getU() const
Definition: GFDetPlane.h:76
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
TVector3 getNormal() const
Definition: GFDetPlane.cxx:138
double extrapolate(const GFDetPlane &, TMatrixT< double > &statePred, TMatrixT< double > &covPred)
returns the tracklength spanned in this extrapolation
Definition: RKTrackRep.cxx:442
TVector3 getO() const
Definition: GFDetPlane.h:75
void RKTrackRep::getPosMomCov ( const GFDetPlane pl,
TVector3 &  pos,
TVector3 &  mom,
TMatrixT< double > &  cov 
)
virtual

method which gets position, momentum and 6x6 covariance matrix

default implementation in cxx file, if a ConcreteTrackRep can not implement this functionality

Reimplemented from GFAbsTrackRep.

Definition at line 1148 of file RKTrackRep.cxx.

References extrapolate(), fCacheSpu, fCharge, GFAbsTrackRep::fCov, GFAbsTrackRep::fRefPlane, fSpu, GFAbsTrackRep::fState, GFDetPlane::getO(), GFDetPlane::getU(), GFDetPlane::getV(), p, and W.

Referenced by PndLmdBPRungeKuttaTask::Exec().

1148  {
1149 
1150  const TVector3 U=pl.getU();
1151  const TVector3 V=pl.getV();
1152  const TVector3 W=U.Cross(V);
1153 
1154  TMatrixT<double> statePred=fState;
1155  TMatrixT<double> covPred=fCov;
1156  double spu=fSpu;
1157  if(pl!=fRefPlane)
1158  {
1159  extrapolate(pl, statePred, covPred);
1160  spu=fCacheSpu;
1161  }
1162 
1163  const double p=fCharge/statePred[0][0]; //Magnitude of the momentum.
1164  const TVector3 pTilde=spu*(W+statePred[1][0]*U+statePred[2][0]*V);
1165  const double pTildeMagnitude=pTilde.Mag();
1166 
1167  TMatrixT<double> jacobian(6, 5); //Jacobian matrix \frac{\partial\left(x, y, z, p_x, p_y, p_z\right)}{\partial\left(\frac{q}{p}, u^\prime, v^\prime, u, v\right)}
1168  jacobian[0][3]=U.X();
1169  jacobian[0][4]=V.X();
1170  jacobian[1][3]=U.Y();
1171  jacobian[1][4]=V.Y();
1172  jacobian[2][3]=U.Z();
1173  jacobian[2][4]=V.Z();
1174  jacobian[3][1]=p*spu/pTildeMagnitude*(U.X()-pTilde.X()*(U*pTilde)/(pTildeMagnitude*pTildeMagnitude));
1175  jacobian[3][2]=p*spu/pTildeMagnitude*(V.X()-pTilde.X()*(V*pTilde)/(pTildeMagnitude*pTildeMagnitude));
1176  jacobian[4][1]=p*spu/pTildeMagnitude*(U.Y()-pTilde.Y()*(U*pTilde)/(pTildeMagnitude*pTildeMagnitude));
1177  jacobian[4][2]=p*spu/pTildeMagnitude*(V.Y()-pTilde.Y()*(V*pTilde)/(pTildeMagnitude*pTildeMagnitude));
1178  jacobian[5][1]=p*spu/pTildeMagnitude*(U.Z()-pTilde.Z()*(U*pTilde)/(pTildeMagnitude*pTildeMagnitude));
1179  jacobian[5][2]=p*spu/pTildeMagnitude*(V.Z()-pTilde.Z()*(V*pTilde)/(pTildeMagnitude*pTildeMagnitude));
1180 
1181  TMatrixT<double> transposedJacobian=jacobian;
1182  transposedJacobian.T();
1183 
1184  cov.ResizeTo(6, 6);
1185 
1186  pos=pl.getO()+statePred[3][0]*U+statePred[4][0]*V;
1187  mom=p/pTildeMagnitude*pTilde;
1188  cov=jacobian*covPred*transposedJacobian;
1189 }
TVector3 pos
Double_t p
Definition: anasim.C:58
TVector3 getV() const
Definition: GFDetPlane.h:77
Double_t mom
Definition: plot_dirc.C:14
#define W
Definition: createSTT.C:76
GFDetPlane fRefPlane
TVector3 getU() const
Definition: GFDetPlane.h:76
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
double extrapolate(const GFDetPlane &, TMatrixT< double > &statePred, TMatrixT< double > &covPred)
returns the tracklength spanned in this extrapolation
Definition: RKTrackRep.cxx:442
TVector3 getO() const
Definition: GFDetPlane.h:75
TMatrixT< double > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:91
void GFAbsTrackRep::getPosMomCov ( TVector3 &  pos,
TVector3 &  mom,
TMatrixT< double > &  c 
)
inlineinherited

Definition at line 222 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fRefPlane, and GFAbsTrackRep::getPosMomCov().

222  {
224  }
TVector3 pos
Double_t mom
Definition: plot_dirc.C:14
virtual void getPosMomCov(const GFDetPlane &pl, TVector3 &pos, TVector3 &mom, TMatrixT< double > &cov)
method which gets position, momentum and 6x6 covariance matrix
GFDetPlane fRefPlane
double GFAbsTrackRep::getRedChiSqu ( ) const
inlineinherited

returns chi2/ndf

Definition at line 248 of file GFAbsTrackRep.h.

References GFAbsTrackRep::getChiSqu(), and GFAbsTrackRep::getNDF().

Referenced by GFTrack::getRedChiSqu().

248  {
249  if(getNDF()>0) return getChiSqu()/getNDF();
250  return 0;
251  }
double getChiSqu() const
unsigned int getNDF() const
const GFDetPlane& GFAbsTrackRep::getReferencePlane ( ) const
inlineinherited

Definition at line 295 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fRefPlane.

Referenced by GFKalman::processHit(), and GFKalman::processTrack().

295 {return fRefPlane;}
GFDetPlane fRefPlane
TMatrixT<double> GFAbsTrackRep::getState ( ) const
inlineinherited

Definition at line 195 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fState.

Referenced by GFTrack::fillGeoTrack(), GFTrack::getResiduals(), GFKalman::processHit(), and GFKalman::processTrack().

195  {
196  return fState;
197  }
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
double GFAbsTrackRep::getStateElem ( int  i) const
inlineinherited

Definition at line 202 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fState.

202 {return fState(i,0);}
Int_t i
Definition: run_full.C:25
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
bool GFAbsTrackRep::getStatusFlag ( )
inlineinherited

Definition at line 318 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fStatusFlag.

Referenced by GFTrack::blowUpCovs(), PndHypMicroWriter::Exec(), PndHypMicroIdealWriter::Exec(), PndHypDKalmanTask::Exec(), GFKalman::fittingPass(), for(), and GFDaf::processTrack().

318  {
319  return fStatusFlag;
320  }
int fStatusFlag
status of track representation: 0 means everything&#39;s OK
Definition: GFAbsTrackRep.h:98
bool RKTrackRep::hasAuxInfo ( )
inlinevirtual

See if the track representation has auxillary information stored.

See if auxillary information is stored in the track representation. See the documentation of GFAbsTrackRep::getAuxInfo() for details.

Reimplemented from GFAbsTrackRep.

Definition at line 162 of file tracking/GenfitTools/trackrep/RKTrackRep/RKTrackRep.h.

162 { return true; }
RKTrackRep& RKTrackRep::operator= ( const RKTrackRep )
inlineprivate

Definition at line 168 of file tracking/GenfitTools/trackrep/RKTrackRep/RKTrackRep.h.

168 {return *this;}; // rhs // [R.K.03/2017] unused variable(s)
TVector3 RKTrackRep::poca2Line ( const TVector3 &  extr1,
const TVector3 &  extr2,
const TVector3 &  point 
) const
private

Definition at line 376 of file RKTrackRep.cxx.

References t.

Referenced by extrapolateToLine().

376  {
377 
378  TVector3 theWire = extr2-extr1;
379  if(theWire.Mag()<1.E-8){
380  GFException exc("RKTrackRep::poca2Line ==> try to find poca between line and point, but the line is really just a point",__LINE__,__FILE__);
381  throw exc;
382  }
383  double t = 1./(theWire*theWire)*(point*theWire+extr1*extr1-extr1*extr2);
384  return (extr1+t*theWire);
385 }
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:45
TTree * t
Definition: bump_analys.C:13
PndSdsMCPoint * point
Definition: anaLmdCluster.C:72
void GFAbsTrackRep::Print ( ) const
virtualinherited

Definition at line 97 of file GFAbsTrackRep.cxx.

References GFAbsTrackRep::fChiSqu, GFAbsTrackRep::fCov, GFAbsTrackRep::fRefPlane, GFAbsTrackRep::fState, and GFDetPlane::Print().

Referenced by PndLmdBPRungeKuttaTask::Exec(), PndHypDKalmanTask::Exec(), and GFTrack::Print().

97  {
98  std::cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
99  std::cout<<"GFAbsTrackRep::Parameters at reference plane ";
100  fRefPlane.Print();
101  std::cout<<"GFAbsTrackRep::State"<<std::endl;
102  fState.Print();
103  std::cout<<"GFAbsTrackRep::Covariances"<<std::endl;
104  fCov.Print();
105  std::cout<<"GFAbsTrackRep::chi^2"<<std::endl;
106  std::cout<<fChiSqu<<std::endl;
107  std::cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
108 }
GFDetPlane fRefPlane
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
double fChiSqu
chiSqu of the track fit
Definition: GFAbsTrackRep.h:94
void Print() const
Definition: GFDetPlane.cxx:238
TMatrixT< double > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:91
virtual GFAbsTrackRep* RKTrackRep::prototype ( ) const
inlinevirtual

Implements GFAbsTrackRep.

Definition at line 75 of file tracking/GenfitTools/trackrep/RKTrackRep/RKTrackRep.h.

References RKTrackRep().

75 {return new RKTrackRep();}
void GFAbsTrackRep::reset ( )
virtualinherited

Definition at line 84 of file GFAbsTrackRep.cxx.

References GFAbsTrackRep::fCov, GFAbsTrackRep::fFirstCov, GFAbsTrackRep::fFirstState, GFAbsTrackRep::fLastCov, GFAbsTrackRep::fLastState, GFAbsTrackRep::fRefPlane, GFAbsTrackRep::fState, and GFDetPlane::set().

84  {
85  std::cout<<"GFAbsTrackRep::reset"<<std::endl;
86  TVector3 nullVec(0.,0.,0.);
87  fRefPlane.set(nullVec,nullVec,nullVec);
88  fState.Zero();
89  fCov.Zero();
90  fFirstState.Zero();
91  fFirstCov.Zero();
92  fLastState.Zero();
93  fLastCov.Zero();
94 }
void set(const TVector3 &o, const TVector3 &u, const TVector3 &v)
Definition: GFDetPlane.cxx:83
TMatrixT< double > fLastState
TMatrixT< double > fFirstCov
TMatrixT< double > fFirstState
state, cov and plane for first and last point in fit
GFDetPlane fRefPlane
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
TMatrixT< double > fLastCov
TMatrixT< double > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:91
bool RKTrackRep::RKutta ( const GFDetPlane plane,
double *  P,
double &  coveredDistance,
std::vector< TVector3 > &  points,
std::vector< double > &  pointLengths,
const double &  maxLen = -1,
bool  calcCov = true 
) const
private

Propagates the particle through the magnetic field.

If the propagation is successfull and the plane is reached, the function returns true. The argument P has to contain the state (P[0] - P[6]) and a unity matrix (P[7] - P[55]) with the last column multiplied wit q/p (hence P[55] is not 1 but q/p). Propagated state and the jacobian (with the last column multiplied wit q/p) of the extrapolation are written to P. In the main loop of the Runge Kutta algorithm, the steppers in #fEffect are called and may reduce the estimated stepsize so that a maximum momentum loss will not be exceeded. If this is the case, RKutta() will only propagate the reduced distance and then return. This is to ensure that material effects, which are calculated after the propagation, are taken into account properly.

Definition at line 640 of file RKTrackRep.cxx.

References GFDetPlane::distance(), error(), fabs(), fCharge, fPdg, GFFieldManager::getFieldVal(), GFMaterialEffects::getInstance(), GFDetPlane::getNormal(), GFDetPlane::getO(), i, GFDetPlane::inActive(), MINSTEP, pos, r, R, sqrt(), GFMaterialEffects::stepper(), and W.

Referenced by Extrap().

646  {
647 
648  static const double EC = .000149896229; // c/(2*10^12) resp. c/2Tera
649  static const double DLT = .0002; // max. deviation for approximation-quality test
650  static const double DLT32 = DLT/32.; //
651  static const double P3 = 1./3.; // 1/3
652  static const double Smax = 100.; // max. step allowed > 0
653  static const double Wmax = 3000.; // max. way allowed
654  static const double Pmin = 4.E-3; // minimum momentum for propagation [GeV]
655  static const int ND = 56; // number of variables for derivatives calculation
656  static const int ND1 = ND-7; // = 49
657  double* R = &P[0]; // Start coordinates in cm ( x, y, z)
658  double* A = &P[3]; // Start directions (ax, ay, az); ax^2+ay^2+az^2=1
659  double SA[3] = {0.,0.,0.}; // Start directions derivatives
660  double Pinv = P[6]*EC; // P[6] is charge/momentum in e/(Gev/c)
661  double Way = 0.; // Total way of the trajectory
662  double Way2 = 0.; // Total way of the trajectory with correct signs
663  bool error = false; // Error of propogation
664  bool stopBecauseOfMaterial = false; // does not go through main loop again when stepsize is reduced by stepper
665 
666  points.clear();
667  pointPaths.clear();
668  //std::cout<<"coords "<<R[0]<<" "<<R[1]<<" "<<R[2]<< std::endl;
669  //std::cout<<"R "<<sqrt(pow(R[0],2)+pow(R[1],2))<< std::endl;
670  //std::cout<<"momentum "<<fabs(fCharge/P[6])<< std::endl;
671  if(fabs(fCharge/P[6])<Pmin){
672  std::cerr << "RKTrackRep::RKutta ==> momentum too low: " << fabs(fCharge/P[6])*1000. << " MeV" << std::endl;
673  return (false);
674  }
675 
676  double SU[4];
677  TVector3 O = plane.getO();
678  TVector3 W = plane.getNormal();
679  if(W*O > 0){ // make SU vector point away from origin
680  SU[0] = W.X();
681  SU[1] = W.Y();
682  SU[2] = W.Z();
683  }
684  else{
685  SU[0] = -1.*W.X();
686  SU[1] = -1.*W.Y();
687  SU[2] = -1.*W.Z();
688  }
689  SU[3] = plane.distance(0.,0.,0.);
690 
691  //
692  // Step estimation until surface
693  //
694  double Step,An,Dist,Dis,S,Sl=0;
695 
696  points.push_back(TVector3(R[0],R[1],R[2]));
697  pointPaths.push_back(0.);
698 
699  An=A[0]*SU[0]+A[1]*SU[1]+A[2]*SU[2]; // An = A * N; component of A normal to surface
700  if(fabs(An) < 1.E-6) {
701  std::cerr << "RKTrackRep::RKutta ==> cannot propagate perpendicular to plane " << std::endl;
702  return false; // no propagation if A perpendicular to surface
703  }
704  if( plane.inActive(TVector3(R[0],R[1],R[2]),TVector3(A[0],A[1],A[2]))) { // if direction is not pointing to active part of surface
705  Dist=SU[3]-R[0]*SU[0]-R[1]*SU[1]-R[2]*SU[2]; // Distance between start coordinates and surface
706  Step=Dist/An;
707  }
708  else{ // make sure to extrapolate towards surface
709  if( (O.X()-R[0])*A[0] + (O.Y()-R[1])*A[1] + (O.Z()-R[2])*A[2] >0 ){ // if direction A pointing from start coordinates R towards surface
710  Dist = sqrt((R[0]-O.X())*(R[0]-O.X())+ // |R-O|; Distance between start coordinates and origin of surface
711  (R[1]-O.Y())*(R[1]-O.Y())+
712  (R[2]-O.Z())*(R[2]-O.Z()));
713  }
714  else{ // if direction pointing away from surface
715  Dist = -1.*sqrt((R[0]-O.X())*(R[0]-O.X())+
716  (R[1]-O.Y())*(R[1]-O.Y())+
717  (R[2]-O.Z())*(R[2]-O.Z()));
718  }
719  Step=Dist;
720  }
721 
722  if(fabs(Step)>Wmax) {
723  std::cerr<<"RKTrackRep::RKutta ==> Too long extrapolation requested : "<<Step<<" cm !"<<std::endl;
724  std::cerr<<"X = "<<R[0]<<" Y = "<<R[1]<<" Z = "<<R[2]
725  <<" COSx = "<<A[0]<<" COSy = "<<A[1]<<" COSz = "<<A[2]<<std::endl;
726  std::cout<<"Destination X = "<<SU[0]*SU[3]<<std::endl;
727  return(false);
728  }
729 
730  // reduce maximum stepsize S to Smax
731  Step>Smax ? S=Smax : Step<-Smax ? S=-Smax : S=Step;
732 
733  //
734  // Main cycle of Runge-Kutta method
735  //
736 
737  //for saving points only when the direction didnt change
738  int Ssign=1;
739  if(S<0) Ssign = -1;
740 
741  while(fabs(Step)>MINSTEP && !stopBecauseOfMaterial) {
742 
743  // call stepper and reduce stepsize
744  double stepperLen;
745  stepperLen = GFMaterialEffects::getInstance()->stepper(fabs(S),
746  R[0],R[1],R[2],
747  Ssign*A[0],Ssign*A[1],Ssign*A[2],
748  fabs(fCharge/P[6]),
749  fPdg);
750  if (stepperLen<MINSTEP) stepperLen=MINSTEP; // prevents tiny stepsizes that can slow down the program
751  if (S > stepperLen) {
752  S = stepperLen;
753  stopBecauseOfMaterial = true;
754  }
755  else if (S < -stepperLen) {
756  S = -stepperLen;
757  stopBecauseOfMaterial = true;
758  }
759 
760  double H0[12],H1[12],H2[12],r[3];
761  double S3=P3*S, S4=.25*S, PS2=Pinv*S;
762 
763  //
764  // First point
765  //
766  r[0]=R[0] ; r[1]=R[1] ; r[2]=R[2] ;
767  TVector3 pos(r[0],r[1],r[2]); // vector of start coordinates R0 (x, y, z)
768  TVector3 H0vect = GFFieldManager::getFieldVal(pos); // magnetic field in 10^-4 T = kGauss
769  H0[0]=PS2*H0vect.X(); H0[1]=PS2*H0vect.Y(); H0[2]=PS2*H0vect.Z(); // H0 is PS2*(Hx, Hy, Hz) @ R0
770  double A0=A[1]*H0[2]-A[2]*H0[1], B0=A[2]*H0[0]-A[0]*H0[2], C0=A[0]*H0[1]-A[1]*H0[0]; // (ax, ay, az) x H0
771  double A2=A[0]+A0 , B2=A[1]+B0 , C2=A[2]+C0 ; // (A0, B0, C0) + (ax, ay, az)
772  double A1=A2+A[0] , B1=B2+A[1] , C1=C2+A[2] ; // (A0, B0, C0) + 2*(ax, ay, az)
773 
774  //
775  // Second point
776  //
777  r[0]+=A1*S4 ; r[1]+=B1*S4 ; r[2]+=C1*S4 ; //setup.Field(r,H1);
778  pos.SetXYZ(r[0],r[1],r[2]);
779  TVector3 H1vect = GFFieldManager::getFieldVal(pos);
780  H1[0]=H1vect.X()*PS2; H1[1]=H1vect.Y()*PS2;H1[2]=H1vect.Z()*PS2; // H1 is PS2*(Hx, Hy, Hz) @ (x, y, z) + 0.25*S * [(A0, B0, C0) + 2*(ax, ay, az)]
781  double A3,B3,C3,A4,B4,C4,A5,B5,C5;
782  A3 = B2*H1[2]-C2*H1[1]+A[0]; B3=C2*H1[0]-A2*H1[2]+A[1]; C3=A2*H1[1]-B2*H1[0]+A[2]; // (A2, B2, C2) x H1 + (ax, ay, az)
783  A4 = B3*H1[2]-C3*H1[1]+A[0]; B4=C3*H1[0]-A3*H1[2]+A[1]; C4=A3*H1[1]-B3*H1[0]+A[2]; // (A3, B3, C3) x H1 + (ax, ay, az)
784  A5 = A4-A[0]+A4 ; B5=B4-A[1]+B4 ; C5=C4-A[2]+C4 ; // 2*(A4, B4, C4) - (ax, ay, az)
785 
786  //
787  // Last point
788  //
789  r[0]=R[0]+S*A4 ; r[1]=R[1]+S*B4 ; r[2]=R[2]+S*C4 ; //setup.Field(r,H2);
790  pos.SetXYZ(r[0],r[1],r[2]);
791  TVector3 H2vect = GFFieldManager::getFieldVal(pos);
792  H2[0]=H2vect.X()*PS2; H2[1]=H2vect.Y()*PS2;H2[2]=H2vect.Z()*PS2; // H2 is PS2*(Hx, Hy, Hz) @ (x, y, z) + 0.25*S * (A4, B4, C4)
793  double A6=B5*H2[2]-C5*H2[1], B6=C5*H2[0]-A5*H2[2], C6=A5*H2[1]-B5*H2[0]; // (A5, B5, C5) x H2
794 
795  //
796  // Test approximation quality on given step and possible step reduction
797  //
798  double EST = fabs((A1+A6)-(A3+A4))+fabs((B1+B6)-(B3+B4))+fabs((C1+C6)-(C3+C4)); // EST = ||(ABC1+ABC6)-(ABC3+ABC4)||_1 = ||(axzy x H0 + ABC5 x H2) - (ABC2 x H1 + ABC3 x H1)||_1
799  if(EST>DLT) {
800  S*=0.5;
801  stopBecauseOfMaterial = false;
802  continue;
803  }
804 
805  //
806  // Derivatives of track parameters in last point
807  //
808  if(calcCov){
809  for(int i=7; i!=ND; i+=7) { // i = 7, 14, 21, 28, 35, 42, 49; ND = 56; ND1 = 49; rows of Jacobian
810 
811  double* dR = &P[i]; // dR = (dX/dpN, dY/dpN, dZ/dpN)
812  double* dA = &P[i+3]; // dA = (dAx/dpN, dAy/dpN, dAz/dpN); N = X,Y,Z,Ax,Ay,Az,q/p
813 
814  //first point
815  double dA0 = H0[ 2]*dA[1]-H0[ 1]*dA[2]; // dA0/dp }
816  double dB0 = H0[ 0]*dA[2]-H0[ 2]*dA[0]; // dB0/dp } = dA x H0
817  double dC0 = H0[ 1]*dA[0]-H0[ 0]*dA[1]; // dC0/dp }
818 
819  if(i==ND1) {dA0+=A0; dB0+=B0; dC0+=C0;} // if last row: (dA0, dB0, dC0) := (dA0, dB0, dC0) + (A0, B0, C0)
820 
821  double dA2 = dA0+dA[0]; // }
822  double dB2 = dB0+dA[1]; // } = (dA0, dB0, dC0) + dA
823  double dC2 = dC0+dA[2]; // }
824 
825  //second point
826  double dA3 = dA[0]+dB2*H1[2]-dC2*H1[1]; // dA3/dp }
827  double dB3 = dA[1]+dC2*H1[0]-dA2*H1[2]; // dB3/dp } = dA + (dA2, dB2, dC2) x H1
828  double dC3 = dA[2]+dA2*H1[1]-dB2*H1[0]; // dC3/dp }
829 
830  if(i==ND1) {dA3+=A3-A[0]; dB3+=B3-A[1]; dC3+=C3-A[2];} // if last row: (dA3, dB3, dC3) := (dA3, dB3, dC3) + (A3, B3, C3) - (ax, ay, az)
831 
832  double dA4 = dA[0]+dB3*H1[2]-dC3*H1[1]; // dA4/dp }
833  double dB4 = dA[1]+dC3*H1[0]-dA3*H1[2]; // dB4/dp } = dA + (dA3, dB3, dC3) x H1
834  double dC4 = dA[2]+dA3*H1[1]-dB3*H1[0]; // dC4/dp }
835 
836  if(i==ND1) {dA4+=A4-A[0]; dB4+=B4-A[1]; dC4+=C4-A[2];} // if last row: (dA4, dB4, dC4) := (dA4, dB4, dC4) + (A4, B4, C4) - (ax, ay, az)
837 
838  //last point
839  double dA5 = dA4+dA4-dA[0]; // }
840  double dB5 = dB4+dB4-dA[1]; // } = 2*(dA4, dB4, dC4) - dA
841  double dC5 = dC4+dC4-dA[2]; // }
842 
843  double dA6 = dB5*H2[2]-dC5*H2[1]; // dA6/dp }
844  double dB6 = dC5*H2[0]-dA5*H2[2]; // dB6/dp } = (dA5, dB5, dC5) x H2
845  double dC6 = dA5*H2[1]-dB5*H2[0]; // dC6/dp }
846 
847  if(i==ND1) {dA6+=A6; dB6+=B6; dC6+=C6;} // if last row: (dA6, dB6, dC6) := (dA6, dB6, dC6) + (A6, B6, C6)
848 
849  dR[0]+=(dA2+dA3+dA4)*S3; dA[0] = (dA0+dA3+dA3+dA5+dA6)*P3; // dR := dR + S3*[(dA2, dB2, dC2) + (dA3, dB3, dC3) + (dA4, dB4, dC4)]
850  dR[1]+=(dB2+dB3+dB4)*S3; dA[1] = (dB0+dB3+dB3+dB5+dB6)*P3; // dA := 1/3*[(dA0, dB0, dC0) + 2*(dA3, dB3, dC3) + (dA5, dB5, dC5) + (dA6, dB6, dC6)]
851  dR[2]+=(dC2+dC3+dC4)*S3; dA[2] = (dC0+dC3+dC3+dC5+dC6)*P3;
852  }
853  }
854 
855  Way2 += S; // add stepsize to way (signed)
856  if((Way+=fabs(S))>Wmax){
857  std::cerr<<"RKTrackRep::RKutta ==> Trajectory is longer than length limit : "<<Way<<" cm !"
858  << " p/q = "<<1./P[6]<< " GeV"<<std::endl;
859  return(false);
860  }
861 
862  //
863  // Track parameters in last point
864  //
865  R[0]+=(A2+A3+A4)*S3; A[0]+=(SA[0]=(A0+A3+A3+A5+A6)*P3-A[0]); // R = R0 + S3*[(A2, B2, C2) + (A3, B3, C3) + (A4, B4, C4)]
866  R[1]+=(B2+B3+B4)*S3; A[1]+=(SA[1]=(B0+B3+B3+B5+B6)*P3-A[1]); // A = 1/3*[(A0, B0, C0) + 2*(A3, B3, C3) + (A5, B5, C5) + (A6, B6, C6)]
867  R[2]+=(C2+C3+C4)*S3; A[2]+=(SA[2]=(C0+C3+C3+C5+C6)*P3-A[2]); // SA = A_new - A_old
868  Sl=S; // last S used
869 
870  // if extrapolation has changed direction, delete the last point, because it is
871  // not a consecutive point to be used for material estimations
872  if(Ssign*S<0.) {
873  pointPaths.at(pointPaths.size()-1)+=S;
874  points.erase(points.end());
875  }
876  else{
877  pointPaths.push_back(S);
878  }
879 
880  points.push_back(TVector3(R[0],R[1],R[2]));
881 
882  double CBA = 1./sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2]); // 1/|A|
883  A[0]*=CBA; A[1]*=CBA; A[2]*=CBA; // normalize A
884 
885  // Step estimation until surface and test conditions for stop of propogation
886  if(fabs(Way2)>Wmax) {
887  Dis=0.;
888  Dist=0.;
889  S=0;
890  Step=0.;
891  break;
892  }
893 
894 
895  An=A[0]*SU[0]+A[1]*SU[1]+A[2]*SU[2];
896 
897  if(fabs(An) < 1.E-6) {
898  error=true;
899  Step=0;
900  break;
901  }
902 
903  if( plane.inActive(TVector3(R[0],R[1],R[2]),TVector3(A[0],A[1],A[2]))) {
904  Dis=SU[3]-R[0]*SU[0]-R[1]*SU[1]-R[2]*SU[2];
905  Step=Dis/An;
906  }
907  else{
908  if( (O.X()-R[0])*A[0] + (O.Y()-R[1])*A[1] + (O.Z()-R[2])*A[2] >0 ){
909  Dis = sqrt((R[0]-O.X())*(R[0]-O.X())+
910  (R[1]-O.Y())*(R[1]-O.Y())+
911  (R[2]-O.Z())*(R[2]-O.Z()));
912  }
913  else{
914  Dis = -1.*sqrt((R[0]-O.X())*(R[0]-O.X())+
915  (R[1]-O.Y())*(R[1]-O.Y())+
916  (R[2]-O.Z())*(R[2]-O.Z()));
917  }
918  Step = Dis; // signed distance to surface
919  }
920 
921  if (Dis*Dist>0 && fabs(Dis)>fabs(Dist)) { // did not get closer to surface
922  error=true;
923  Step=0;
924  break;
925  }
926  Dist=Dis;
927 
928  //
929  // reset & check step size
930  //
931  // reset S to Step if extrapolation too long or in wrong direction
932  if (S*Step<0. || fabs(S)>fabs(Step)) S=Step;
933  else if (EST<DLT32 && fabs(2.*S)<=Smax) S*=2.;
934 
935  } //end of main loop
936 
937  //
938  // Output information preparation for main track parameteres
939  //
940 
941  if (!stopBecauseOfMaterial) { // linear extrapolation to surface
942  if(fabs(Sl) > 1.E-12) Sl=1./Sl; // Sl = inverted last Stepsize Sl
943  A [0]+=(SA[0]*=Sl)*Step; // Step = distance to surface
944  A [1]+=(SA[1]*=Sl)*Step; // SA*Sl = delta A / delta way; local derivative of A with respect to the length of the way
945  A [2]+=(SA[2]*=Sl)*Step; // A = A + Step * SA*Sl
946 
947  P[0] = R[0]+Step*(A[0]-.5*Step*SA[0]); // P = R + Step*(A - 1/2*Step*SA); approximation for final point on surface
948  P[1] = R[1]+Step*(A[1]-.5*Step*SA[1]);
949  P[2] = R[2]+Step*(A[2]-.5*Step*SA[2]);
950 
951  points.push_back(TVector3(P[0],P[1],P[2]));
952  pointPaths.push_back(Step);
953  }
954 
955  double CBA = 1./sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2]);
956 
957  P[3] = A[0]*CBA; // normalize A
958  P[4] = A[1]*CBA;
959  P[5] = A[2]*CBA;
960 
961  //
962  // Output derivatives of track parameters preparation
963  //
964  An = A[0]*SU[0]+A[1]*SU[1]+A[2]*SU[2];
965  fabs(An) < 1.E-6 ? An=1./An : An = 0; // 1/A_normal
966 
967  if(calcCov && !stopBecauseOfMaterial){
968  for(int i=7; i!=ND; i+=7) {
969  double* dR = &P[i]; double* dA = &P[i+3];
970  S = (dR[0]*SU[0]+dR[1]*SU[1]+dR[2]*SU[2])*An; // dR_normal / A_normal
971  dR[0]-=S*A [0]; dR[1]-=S*A [1]; dR[2]-=S*A [2];
972  dA[0]-=S*SA[0]; dA[1]-=S*SA[1]; dA[2]-=S*SA[2];
973  }
974  }
975 
976  if(error){
977  std::cerr << "RKTrackRep::RKutta ==> Do not get closer. Path = " << Way << " cm" << " p/q = " << 1./P[6] << " GeV" << std::endl;
978  return(false);
979  }
980 
981  // calculate covered distance
982  if (!stopBecauseOfMaterial) coveredDistance=Way2+Step;
983  else coveredDistance=Way2;
984 
985  return(true);
986 }
TVector3 pos
double r
Definition: RiemannTest.C:14
Int_t i
Definition: run_full.C:25
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
double distance(TVector3 &) const
Definition: GFDetPlane.cxx:343
static TVector3 getFieldVal(const TVector3 &x)
static void error(int no)
Definition: ranlxd.cxx:419
#define W
Definition: createSTT.C:76
#define MINSTEP
Definition: RKTrackRep.cxx:40
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
double stepper(const double &maxDist, const double &posx, const double &posy, const double &posz, const double &dirx, const double &diry, const double &dirz, const double &mom, const int &pdg)
Returns maximum length so that a specified momentum loss will not be exceeded.
GeV c P
bool inActive(const TVector3 &point, const TVector3 &dir) const
intersect in the active area? C.f. GFAbsFinitePlane
Definition: GFDetPlane.h:132
static GFMaterialEffects * getInstance()
TVector3 getNormal() const
Definition: GFDetPlane.cxx:138
Double_t R
Definition: checkhelixhit.C:61
TVector3 getO() const
Definition: GFDetPlane.h:75
void GFAbsTrackRep::setChiSqu ( double  aChiSqu)
inlineinherited

Definition at line 297 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fChiSqu.

Referenced by GFKalman::fittingPass().

297  {
298  fChiSqu = aChiSqu;
299  }
double fChiSqu
chiSqu of the track fit
Definition: GFAbsTrackRep.h:94
void GFAbsTrackRep::setCov ( const TMatrixT< double > &  aCov)
inlineinherited

Definition at line 273 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fCov.

Referenced by GFTrack::blowUpCovs().

273  {
274  fCov=aCov;
275  }
TMatrixT< double > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:91
void RKTrackRep::setData ( const TMatrixT< double > &  st,
const GFDetPlane pl,
const TMatrixT< double > *  cov = NULL,
const TMatrixT< double > *  aux = NULL 
)
virtual

Sets state, plane and (optionally) covariance.

This function also sets the parameter fSpu to the value stored in fCacheSpu. Therefore it has to be ensured that the plane #pl is the same as the plane of the last extrapolation (i.e. fCachePlane), where fCacheSpu was calculated. Hence, if the argument #pl is not equal to fCachePlane, an error message is shown an an exception is thrown.

Reimplemented from GFAbsTrackRep.

Definition at line 43 of file RKTrackRep.cxx.

References fCachePlane, fCacheSpu, fCharge, fSpu, GFAbsTrackRep::fState, and GFAbsTrackRep::setData().

43  {
44  if(aux != NULL) {
45  fCacheSpu = (*aux)(0,0);
46  } else {
47  if(pl!=fCachePlane){
48  std::cerr << "RKTrackRep::setData() - a fatal error ocurred! It was called with a reference plane which is not the same as the one from the last extrapolate(plane,state,cov)-> abort in line " << __LINE__ << std::endl;
49  throw;
50  }
51  }
52  GFAbsTrackRep::setData(st,pl,cov);
53  if (fCharge*fState[0][0] < 0) fCharge *= -1; // set charge accordingly! (fState[0][0] = q/p)
55 }
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
virtual void setData(const TMatrixT< double > &st, const GFDetPlane &pl, const TMatrixT< double > *cov=NULL, const TMatrixT< double > *aux=NULL)
Puts the track representation in a given state.
void GFAbsTrackRep::setFirstCov ( const TMatrixT< double > &  aCov)
inlineinherited

Definition at line 279 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fFirstCov.

279  {
280  fFirstCov = aCov;
281  }
TMatrixT< double > fFirstCov
void GFAbsTrackRep::setFirstPlane ( const GFDetPlane aPlane)
inlineinherited

Definition at line 282 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fFirstPlane.

282  {
283  fFirstPlane = aPlane;;
284  }
GFDetPlane fFirstPlane
void GFAbsTrackRep::setFirstState ( const TMatrixT< double > &  aState)
inlineinherited

Definition at line 276 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fFirstState.

276  {
277  fFirstState = aState;
278  }
TMatrixT< double > fFirstState
state, cov and plane for first and last point in fit
bool GFAbsTrackRep::setInverted ( bool  f = true)
inlineinherited

Deprecated. Should be removed soon.

Definition at line 316 of file GFAbsTrackRep.h.

References f, and GFAbsTrackRep::fInverted.

316 {fInverted=f; return true;}
bool fInverted
specifies the direction of flight of the particle
TFile * f
Definition: bump_analys.C:12
void GFAbsTrackRep::setLastCov ( const TMatrixT< double > &  aCov)
inlineinherited

Definition at line 288 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fLastCov.

288  {
289  fLastCov = aCov;
290  }
TMatrixT< double > fLastCov
void GFAbsTrackRep::setLastPlane ( const GFDetPlane aPlane)
inlineinherited

Definition at line 291 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fLastPlane.

291  {
292  fLastPlane = aPlane;;
293  }
GFDetPlane fLastPlane
void GFAbsTrackRep::setLastState ( const TMatrixT< double > &  aState)
inlineinherited

Definition at line 285 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fLastState.

285  {
286  fLastState = aState;
287  }
TMatrixT< double > fLastState
void GFAbsTrackRep::setNDF ( unsigned int  n)
inlineinherited

Definition at line 300 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fNdf, and n.

Referenced by GFKalman::fittingPass().

300  {
301  fNdf = n;
302  }
int n
unsigned int fNdf
Definition: GFAbsTrackRep.h:95
void RKTrackRep::setPDG ( int  i)

Set PDG particle code.

Definition at line 274 of file RKTrackRep.cxx.

References exit(), fCharge, fMass, fPdg, and i.

Referenced by RKTrackRep().

274  {
275  fPdg = i;
276  TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(fPdg);
277  if(part == 0){
278  std::cerr << "RKTrackRep::setPDG particle " << i
279  << " not known to TDatabasePDG -> abort" << std::endl;
280  exit(1);
281  }
282  fMass = part->Mass();
283  fCharge = part->Charge()/(3.);
284 }
Int_t i
Definition: run_full.C:25
exit(0)
void GFAbsTrackRep::setStatusFlag ( int  _val)
inlineinherited

Definition at line 309 of file GFAbsTrackRep.h.

References GFAbsTrackRep::fStatusFlag.

Referenced by GFKalman::fittingPass(), and GFDaf::processTrack().

309  {
310  fStatusFlag = _val;
311  }
int fStatusFlag
status of track representation: 0 means everything&#39;s OK
Definition: GFAbsTrackRep.h:98
double RKTrackRep::stepalong ( double  h,
TVector3 &  point,
TVector3 &  dir 
)
virtual

make step of h cm along the track, returns the tracklength spanned in this extrapolation

Also returns the position and direction by reference. It does NOT alter the state of the trackrep and starts extrapolating from fRefPlane.

Reimplemented from GFAbsTrackRep.

Definition at line 1192 of file RKTrackRep.cxx.

References Extrap(), fabs(), GFAbsTrackRep::fRefPlane, fSpu, GFAbsTrackRep::fState, GFDetPlane::getO(), GFDetPlane::getU(), GFDetPlane::getV(), MINSTEP, point, GFDetPlane::setON(), and v.

1192  {
1193 
1194  TVector3 dest;
1195 
1196  static const int maxIt(30);
1197 
1198  TVector3 o=fRefPlane.getO();
1199  TVector3 u=fRefPlane.getU();
1200  TVector3 v=fRefPlane.getV();
1201  TVector3 w=u.Cross(v);
1202 
1203  TVector3 pTilde = fSpu* (w + fState[1][0] * u + fState[2][0] * v);
1204  pTilde.SetMag(1.);
1205 
1206  TVector3 point = o + fState[3][0]*u + fState[4][0]*v;
1207 
1208  TMatrixT<double> state7(7,1);
1209  state7[0][0] = point.X();
1210  state7[1][0] = point.Y();
1211  state7[2][0] = point.Z();
1212  state7[3][0] = pTilde.X();
1213  state7[4][0] = pTilde.Y();
1214  state7[5][0] = pTilde.Z();
1215  state7[6][0] = fState[0][0];
1216 
1217  double coveredDistance(0.);
1218 
1219  GFDetPlane pl;
1220  int iterations(0);
1221 
1222  while(true){
1223  pos.SetXYZ(state7[0][0], state7[1][0], state7[2][0]);
1224  dir.SetXYZ(state7[3][0], state7[4][0], state7[5][0]);
1225  dir.SetMag(1.);
1226 
1227  dest = pos + (h - coveredDistance) * dir;
1228 
1229  pl.setON(dest, dir);
1230  coveredDistance += this->Extrap(pl,&state7);
1231 
1232  if(fabs(h - coveredDistance)<MINSTEP) break;
1233  if(++iterations == maxIt) {
1234  GFException exc("RKTrackRep::stepalong ==> maximum number of iterations reached",__LINE__,__FILE__);
1235  throw exc;
1236  }
1237  }
1238 
1239  pos.SetXYZ(state7[0][0], state7[1][0], state7[2][0]);
1240  dir.SetXYZ(state7[3][0], state7[4][0], state7[5][0]);
1241  dir.SetMag(1.);
1242 
1243  return coveredDistance;
1244 
1245 }
TVector3 pos
TVector3 getV() const
Definition: GFDetPlane.h:77
double Extrap(const GFDetPlane &plane, TMatrixT< double > *state, TMatrixT< double > *cov=NULL) const
Handles propagation and material effects.
Definition: RKTrackRep.cxx:991
Detector plane genfit geometry class.
Definition: GFDetPlane.h:59
void setON(const TVector3 &o, const TVector3 &n)
Definition: GFDetPlane.cxx:145
__m128 v
Definition: P4_F32vec4.h:4
GFDetPlane fRefPlane
#define MINSTEP
Definition: RKTrackRep.cxx:40
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:45
TVector3 getU() const
Definition: GFDetPlane.h:76
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
TVector3 getO() const
Definition: GFDetPlane.h:75
PndSdsMCPoint * point
Definition: anaLmdCluster.C:72
void RKTrackRep::switchDirection ( )
inlinevirtual

Member Data Documentation

TMatrixT<double> RKTrackRep::fAuxInfo
private

Definition at line 211 of file tracking/GenfitTools/trackrep/RKTrackRep/RKTrackRep.h.

Referenced by getAuxInfo().

GFDetPlane RKTrackRep::fCachePlane
private
double RKTrackRep::fCacheSpu
private
double RKTrackRep::fCharge
private
double GFAbsTrackRep::fChiSqu
protectedinherited

chiSqu of the track fit

Definition at line 94 of file GFAbsTrackRep.h.

Referenced by GFAbsTrackRep::addChiSqu(), GFAbsTrackRep::getChiSqu(), GFAbsTrackRep::Print(), and GFAbsTrackRep::setChiSqu().

TMatrixT<double> GFAbsTrackRep::fCov
protectedinherited
unsigned int GFAbsTrackRep::fDimension
protectedinherited

Dimensionality of track representation.

Definition at line 85 of file GFAbsTrackRep.h.

Referenced by GeaneTrackRep::extrapolate(), GFAbsTrackRep::extrapolate(), and GFAbsTrackRep::getDim().

bool RKTrackRep::fDirection
private
TMatrixT<double> GFAbsTrackRep::fFirstCov
protectedinherited
GFDetPlane GFAbsTrackRep::fFirstPlane
protectedinherited

Definition at line 108 of file GFAbsTrackRep.h.

Referenced by GFAbsTrackRep::getFirstPlane(), and GFAbsTrackRep::setFirstPlane().

TMatrixT<double> GFAbsTrackRep::fFirstState
protectedinherited

state, cov and plane for first and last point in fit

Definition at line 103 of file GFAbsTrackRep.h.

Referenced by GFAbsTrackRep::getFirstState(), GFAbsTrackRep::reset(), and GFAbsTrackRep::setFirstState().

bool GFAbsTrackRep::fInverted
protectedinherited

specifies the direction of flight of the particle

Definition at line 100 of file GFAbsTrackRep.h.

Referenced by GFAbsTrackRep::setInverted().

TMatrixT<double> GFAbsTrackRep::fLastCov
protectedinherited
GFDetPlane GFAbsTrackRep::fLastPlane
protectedinherited

Definition at line 109 of file GFAbsTrackRep.h.

Referenced by GFAbsTrackRep::getLastPlane(), and GFAbsTrackRep::setLastPlane().

TMatrixT<double> GFAbsTrackRep::fLastState
protectedinherited
double RKTrackRep::fMass
private

Mass (in GeV)

Definition at line 204 of file tracking/GenfitTools/trackrep/RKTrackRep/RKTrackRep.h.

Referenced by setPDG().

unsigned int GFAbsTrackRep::fNdf
protectedinherited
int RKTrackRep::fPdg
private

PDG particle code.

Definition at line 202 of file tracking/GenfitTools/trackrep/RKTrackRep/RKTrackRep.h.

Referenced by Extrap(), getPDG(), RKutta(), and setPDG().

GFDetPlane GFAbsTrackRep::fRefPlane
protectedinherited
double RKTrackRep::fSpu
private
TMatrixT<double> GFAbsTrackRep::fState
protectedinherited
int GFAbsTrackRep::fStatusFlag
protectedinherited

status of track representation: 0 means everything's OK

Definition at line 98 of file GFAbsTrackRep.h.

Referenced by GFAbsTrackRep::getStatusFlag(), and GFAbsTrackRep::setStatusFlag().


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