LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
genf::RKTrackRep Class Reference

#include "RKTrackRep.h"

Inheritance diagram for genf::RKTrackRep:
genf::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_t > &statePred, TMatrixT< Double_t > &covPred)
 returns the tracklength spanned in this extrapolation More...
 
double extrapolate (const GFDetPlane &, TMatrixT< Double_t > &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...
 
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...
 
TVector3 getMomLast (const GFDetPlane &)
 
void getPosMom (const GFDetPlane &, TVector3 &pos, TVector3 &mom)
 Gets position and momentum in the plane by exrapolating or not. More...
 
double getCharge () const
 Returns charge. More...
 
void switchDirection ()
 deprecated More...
 
void setPDG (int)
 Set PDG particle code. More...
 
int getPDG ()
 
void rescaleCovOffDiags ()
 
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)
 
bool hasAuxInfo ()
 
double extrapolate (const GFDetPlane &plane)
 This changes the state and cov and plane of the rep. More...
 
virtual void stepalong (double h)
 make step of h cm along the track More...
 
unsigned int getDim () const
 returns dimension of state vector More...
 
virtual void Print (std::ostream &out=std::cout) const
 
const TMatrixT< Double_t > & getState () const
 
const TMatrixT< Double_t > & getCov () const
 
double getStateElem (int i) const
 
double getCovElem (int i, int j) const
 
TVector3 getPos ()
 
TVector3 getMom ()
 
virtual void getPosMomCov (const GFDetPlane &pl, TVector3 &pos, TVector3 &mom, TMatrixT< Double_t > &cov)
 method which gets position, momentum and 6x6 covariance matrix More...
 
void getPosMomCov (TVector3 &pos, TVector3 &mom, TMatrixT< Double_t > &c)
 
TMatrixT< Double_t > getFirstState () const
 
TMatrixT< Double_t > getFirstCov () const
 
GFDetPlane getFirstPlane () const
 
TMatrixT< Double_t > getLastState () const
 
TMatrixT< Double_t > getLastCov () const
 
GFDetPlane getLastPlane () const
 
double getChiSqu () const
 
double getRedChiSqu () const
 returns chi2/ndf More...
 
unsigned int getNDF () const
 
virtual void setData (const TMatrixT< Double_t > &st, const GFDetPlane &pl, const TMatrixT< Double_t > *cov=NULL)
 
void setCov (const TMatrixT< Double_t > &aCov)
 
void setFirstState (const TMatrixT< Double_t > &aState)
 
void setFirstCov (const TMatrixT< Double_t > &aCov)
 
void setFirstPlane (const GFDetPlane &aPlane)
 
void setLastState (const TMatrixT< Double_t > &aState)
 
void setLastCov (const TMatrixT< Double_t > &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_t > fState
 The vector of track parameters. More...
 
TMatrixT< Double_t > 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_t > fFirstState
 state, cov and plane for first and last point in fit More...
 
TMatrixT< Double_t > fFirstCov
 
TMatrixT< Double_t > fLastState
 
TMatrixT< Double_t > fLastCov
 
GFDetPlane fFirstPlane
 
GFDetPlane fLastPlane
 
GFDetPlane fRefPlane
 

Private Member Functions

RKTrackRepoperator= (const RKTrackRep *)
 
 RKTrackRep (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
 Contains all material effects. More...
 
TVector3 poca2Line (const TVector3 &extr1, const TVector3 &extr2, const TVector3 &point) const
 
double Extrap (const GFDetPlane &plane, TMatrixT< Double_t > *state, TMatrixT< Double_t > *cov=NULL) const
 Handles propagation and material effects. More...
 

Private Attributes

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

Detailed Description

Definition at line 54 of file RKTrackRep.h.

Constructor & Destructor Documentation

genf::RKTrackRep::RKTrackRep ( )

Definition at line 81 of file RKTrackRep.cxx.

Referenced by clone(), and prototype().

82  : GFAbsTrackRep(5)
83  , fCachePlane()
84  , fCacheSpu(1)
85  , fSpu(1)
86  , fAuxInfo(1, 1)
87  , fDirection(true)
88  , fPdg(0)
89  , fMass(0.)
90  , fCharge(-1)
91 {}
int fPdg
PDG particle code.
Definition: RKTrackRep.h:170
double fMass
Mass (in GeV)
Definition: RKTrackRep.h:172
GFDetPlane fCachePlane
Definition: RKTrackRep.h:160
TMatrixT< double > fAuxInfo
Definition: RKTrackRep.h:163
double fCharge
Charge.
Definition: RKTrackRep.h:174
genf::RKTrackRep::RKTrackRep ( const TVector3 &  pos,
const TVector3 &  mom,
const TVector3 &  poserr,
const TVector3 &  momerr,
const int &  PDGCode 
)

Definition at line 93 of file RKTrackRep.cxx.

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

98  : GFAbsTrackRep(5), fCachePlane(), fCacheSpu(1), fAuxInfo(1, 1), fDirection(true)
99 {
100  setPDG(PDGCode); // also sets charge and mass
101 
102  //fEffect = new GFMaterialEffects();
103 
104  fRefPlane.setO(pos);
105  fRefPlane.setNormal(mom);
106 
107  fState[0][0] = fCharge / mom.Mag();
108  //u' and v'
109  fState[1][0] = 0.0;
110  fState[2][0] = 0.0;
111  //u and v
112  fState[3][0] = 0.0;
113  fState[4][0] = 0.0;
114  //spu
115  fSpu = 1.;
116 
117  TVector3 o = fRefPlane.getO();
118  TVector3 u = fRefPlane.getU();
119  TVector3 v = fRefPlane.getV();
120  TVector3 w = u.Cross(v);
121  double pu = 0.;
122  double pv = 0.;
123  double pw = mom.Mag();
124 
125  fCov[3][3] = poserr.X() * poserr.X() * u.X() * u.X() + poserr.Y() * poserr.Y() * u.Y() * u.Y() +
126  poserr.Z() * poserr.Z() * u.Z() * u.Z();
127  fCov[4][4] = poserr.X() * poserr.X() * v.X() * v.X() + poserr.Y() * poserr.Y() * v.Y() * v.Y() +
128  poserr.Z() * poserr.Z() * v.Z() * v.Z();
129  fCov[0][0] =
130  fCharge * fCharge / pow(mom.Mag(), 6.) *
131  (mom.X() * mom.X() * momerr.X() * momerr.X() + mom.Y() * mom.Y() * momerr.Y() * momerr.Y() +
132  mom.Z() * mom.Z() * momerr.Z() * momerr.Z());
133  fCov[1][1] = pow((u.X() / pw - w.X() * pu / (pw * pw)), 2.) * momerr.X() * momerr.X() +
134  pow((u.Y() / pw - w.Y() * pu / (pw * pw)), 2.) * momerr.Y() * momerr.Y() +
135  pow((u.Z() / pw - w.Z() * pu / (pw * pw)), 2.) * momerr.Z() * momerr.Z();
136  fCov[2][2] = pow((v.X() / pw - w.X() * pv / (pw * pw)), 2.) * momerr.X() * momerr.X() +
137  pow((v.Y() / pw - w.Y() * pv / (pw * pw)), 2.) * momerr.Y() * momerr.Y() +
138  pow((v.Z() / pw - w.Z() * pv / (pw * pw)), 2.) * momerr.Z() * momerr.Z();
139 }
GFDetPlane fCachePlane
Definition: RKTrackRep.h:160
void setO(const TVector3 &o)
Definition: GFDetPlane.cxx:94
TMatrixT< double > fAuxInfo
Definition: RKTrackRep.h:163
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:92
TVector3 getO() const
Definition: GFDetPlane.h:77
void setPDG(int)
Set PDG particle code.
Definition: RKTrackRep.cxx:284
TVector3 getU() const
Definition: GFDetPlane.h:78
void setNormal(TVector3 n)
Definition: GFDetPlane.cxx:151
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
TVector3 getV() const
Definition: GFDetPlane.h:79
double fCharge
Charge.
Definition: RKTrackRep.h:174
Float_t w
Definition: plot.C:20
genf::RKTrackRep::RKTrackRep ( const GFTrackCand aGFTrackCandPtr)

Definition at line 143 of file RKTrackRep.cxx.

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

144  : GFAbsTrackRep(5), fCachePlane(), fCacheSpu(1), fAuxInfo(1, 1), fDirection(true)
145 {
146  setPDG(aGFTrackCandPtr->getPdgCode()); // also sets charge and mass
147 
148  TVector3 mom = aGFTrackCandPtr->getDirSeed();
149  fRefPlane.setO(aGFTrackCandPtr->getPosSeed());
150  fRefPlane.setNormal(mom);
151  double pw = mom.Mag();
152  fState[0][0] = fCharge / pw;
153  //u' and v'
154  fState[1][0] = 0.0;
155  fState[2][0] = 0.0;
156  //u and v
157  fState[3][0] = 0.0;
158  fState[4][0] = 0.0;
159  //spu
160  fSpu = 1.;
161 
162  TVector3 o = fRefPlane.getO();
163  TVector3 u = fRefPlane.getU();
164  TVector3 v = fRefPlane.getV();
165  TVector3 w = u.Cross(v);
166  double pu = 0.;
167  double pv = 0.;
168 
169  TVector3 poserr = aGFTrackCandPtr->getPosError();
170  TVector3 momerr = aGFTrackCandPtr->getDirError();
171  fCov[3][3] = poserr.X() * poserr.X() * u.X() * u.X() + poserr.Y() * poserr.Y() * u.Y() * u.Y() +
172  poserr.Z() * poserr.Z() * u.Z() * u.Z();
173  fCov[4][4] = poserr.X() * poserr.X() * v.X() * v.X() + poserr.Y() * poserr.Y() * v.Y() * v.Y() +
174  poserr.Z() * poserr.Z() * v.Z() * v.Z();
175  fCov[0][0] =
176  fCharge * fCharge / pow(mom.Mag(), 6.) *
177  (mom.X() * mom.X() * momerr.X() * momerr.X() + mom.Y() * mom.Y() * momerr.Y() * momerr.Y() +
178  mom.Z() * mom.Z() * momerr.Z() * momerr.Z());
179  fCov[1][1] = pow((u.X() / pw - w.X() * pu / (pw * pw)), 2.) * momerr.X() * momerr.X() +
180  pow((u.Y() / pw - w.Y() * pu / (pw * pw)), 2.) * momerr.Y() * momerr.Y() +
181  pow((u.Z() / pw - w.Z() * pu / (pw * pw)), 2.) * momerr.Z() * momerr.Z();
182  fCov[2][2] = pow((v.X() / pw - w.X() * pv / (pw * pw)), 2.) * momerr.X() * momerr.X() +
183  pow((v.Y() / pw - w.Y() * pv / (pw * pw)), 2.) * momerr.Y() * momerr.Y() +
184  pow((v.Z() / pw - w.Z() * pv / (pw * pw)), 2.) * momerr.Z() * momerr.Z();
185 }
GFDetPlane fCachePlane
Definition: RKTrackRep.h:160
void setO(const TVector3 &o)
Definition: GFDetPlane.cxx:94
TMatrixT< double > fAuxInfo
Definition: RKTrackRep.h:163
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:92
TVector3 getO() const
Definition: GFDetPlane.h:77
void setPDG(int)
Set PDG particle code.
Definition: RKTrackRep.cxx:284
TVector3 getU() const
Definition: GFDetPlane.h:78
void setNormal(TVector3 n)
Definition: GFDetPlane.cxx:151
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
TVector3 getV() const
Definition: GFDetPlane.h:79
double fCharge
Charge.
Definition: RKTrackRep.h:174
Float_t w
Definition: plot.C:20
genf::RKTrackRep::RKTrackRep ( const TVector3 &  pos,
const TVector3 &  mom,
const int &  PDGCode 
)

Definition at line 187 of file RKTrackRep.cxx.

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

188  : GFAbsTrackRep(5), fCachePlane(), fCacheSpu(1), fAuxInfo(1, 1), fDirection(true)
189 {
190  setPDG(PDGCode); // also sets charge and mass
191 
192  //fEffect = new GFMaterialEffects();
193 
194  fRefPlane.setO(pos);
195  fRefPlane.setNormal(mom);
196 
197  fState[0][0] = fCharge / mom.Mag();
198  //u' and v'
199  fState[1][0] = 0.0;
200  fState[2][0] = 0.0;
201  //u and v
202  fState[3][0] = 0.0;
203  fState[4][0] = 0.0;
204  //spu
205  fSpu = 1.;
206 
207  TVector3 o = fRefPlane.getO();
208  TVector3 u = fRefPlane.getU();
209  TVector3 v = fRefPlane.getV();
210  TVector3 w = u.Cross(v);
211  double pu = 0.;
212  double pv = 0.;
213  double pw = mom.Mag();
214 
215  static const TVector3 stdPosErr(1., 1., 1.);
216  static const TVector3 stdMomErr(1., 1., 1.);
217 
218  fCov[3][3] = stdPosErr.X() * stdPosErr.X() * u.X() * u.X() +
219  stdPosErr.Y() * stdPosErr.Y() * u.Y() * u.Y() +
220  stdPosErr.Z() * stdPosErr.Z() * u.Z() * u.Z();
221  fCov[4][4] = stdPosErr.X() * stdPosErr.X() * v.X() * v.X() +
222  stdPosErr.Y() * stdPosErr.Y() * v.Y() * v.Y() +
223  stdPosErr.Z() * stdPosErr.Z() * v.Z() * v.Z();
224  fCov[0][0] = fCharge * fCharge / pow(mom.Mag(), 6.) *
225  (mom.X() * mom.X() * stdMomErr.X() * stdMomErr.X() +
226  mom.Y() * mom.Y() * stdMomErr.Y() * stdMomErr.Y() +
227  mom.Z() * mom.Z() * stdMomErr.Z() * stdMomErr.Z());
228  fCov[1][1] = pow((u.X() / pw - w.X() * pu / (pw * pw)), 2.) * stdMomErr.X() * stdMomErr.X() +
229  pow((u.Y() / pw - w.Y() * pu / (pw * pw)), 2.) * stdMomErr.Y() * stdMomErr.Y() +
230  pow((u.Z() / pw - w.Z() * pu / (pw * pw)), 2.) * stdMomErr.Z() * stdMomErr.Z();
231  fCov[2][2] = pow((v.X() / pw - w.X() * pv / (pw * pw)), 2.) * stdMomErr.X() * stdMomErr.X() +
232  pow((v.Y() / pw - w.Y() * pv / (pw * pw)), 2.) * stdMomErr.Y() * stdMomErr.Y() +
233  pow((v.Z() / pw - w.Z() * pv / (pw * pw)), 2.) * stdMomErr.Z() * stdMomErr.Z();
234 }
GFDetPlane fCachePlane
Definition: RKTrackRep.h:160
void setO(const TVector3 &o)
Definition: GFDetPlane.cxx:94
TMatrixT< double > fAuxInfo
Definition: RKTrackRep.h:163
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:92
TVector3 getO() const
Definition: GFDetPlane.h:77
void setPDG(int)
Set PDG particle code.
Definition: RKTrackRep.cxx:284
TVector3 getU() const
Definition: GFDetPlane.h:78
void setNormal(TVector3 n)
Definition: GFDetPlane.cxx:151
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
TVector3 getV() const
Definition: GFDetPlane.h:79
double fCharge
Charge.
Definition: RKTrackRep.h:174
Float_t w
Definition: plot.C:20
genf::RKTrackRep::RKTrackRep ( const GFDetPlane pl,
const TVector3 &  mom,
const int &  PDGCode 
)

Definition at line 236 of file RKTrackRep.cxx.

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

237  : GFAbsTrackRep(5), fCachePlane(), fCacheSpu(1), fAuxInfo(1, 1), fDirection(true)
238 {
239  setPDG(PDGCode); // also sets charge and mass
240 
241  //fEffect = new GFMaterialEffects();
242 
243  fRefPlane = pl;
244  TVector3 o = fRefPlane.getO();
245  TVector3 u = fRefPlane.getU();
246  TVector3 v = fRefPlane.getV();
247  TVector3 w = u.Cross(v);
248 
249  double pu = mom * u;
250  double pv = mom * v;
251  double pw = mom * w;
252 
253  fState[0][0] = fCharge / mom.Mag();
254  //u' and v'
255  fState[1][0] = pu / pw;
256  fState[2][0] = pv / pw;
257  //u and v
258  fState[3][0] = 0.0;
259  fState[4][0] = 0.0;
260  //spu
261  fSpu = 1.;
262 
263  static const TVector3 stdPosErr2(1., 1., 1.);
264  static const TVector3 stdMomErr2(1., 1., 1.);
265 
266  fCov[3][3] = stdPosErr2.X() * stdPosErr2.X() * u.X() * u.X() +
267  stdPosErr2.Y() * stdPosErr2.Y() * u.Y() * u.Y() +
268  stdPosErr2.Z() * stdPosErr2.Z() * u.Z() * u.Z();
269  fCov[4][4] = stdPosErr2.X() * stdPosErr2.X() * v.X() * v.X() +
270  stdPosErr2.Y() * stdPosErr2.Y() * v.Y() * v.Y() +
271  stdPosErr2.Z() * stdPosErr2.Z() * v.Z() * v.Z();
272  fCov[0][0] = fCharge * fCharge / pow(mom.Mag(), 6.) *
273  (mom.X() * mom.X() * stdMomErr2.X() * stdMomErr2.X() +
274  mom.Y() * mom.Y() * stdMomErr2.Y() * stdMomErr2.Y() +
275  mom.Z() * mom.Z() * stdMomErr2.Z() * stdMomErr2.Z());
276  fCov[1][1] = pow((u.X() / pw - w.X() * pu / (pw * pw)), 2.) * stdMomErr2.X() * stdMomErr2.X() +
277  pow((u.Y() / pw - w.Y() * pu / (pw * pw)), 2.) * stdMomErr2.Y() * stdMomErr2.Y() +
278  pow((u.Z() / pw - w.Z() * pu / (pw * pw)), 2.) * stdMomErr2.Z() * stdMomErr2.Z();
279  fCov[2][2] = pow((v.X() / pw - w.X() * pv / (pw * pw)), 2.) * stdMomErr2.X() * stdMomErr2.X() +
280  pow((v.Y() / pw - w.Y() * pv / (pw * pw)), 2.) * stdMomErr2.Y() * stdMomErr2.Y() +
281  pow((v.Z() / pw - w.Z() * pv / (pw * pw)), 2.) * stdMomErr2.Z() * stdMomErr2.Z();
282 }
GFDetPlane fCachePlane
Definition: RKTrackRep.h:160
TMatrixT< double > fAuxInfo
Definition: RKTrackRep.h:163
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:92
TVector3 getO() const
Definition: GFDetPlane.h:77
void setPDG(int)
Set PDG particle code.
Definition: RKTrackRep.cxx:284
TVector3 getU() const
Definition: GFDetPlane.h:78
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
TVector3 getV() const
Definition: GFDetPlane.h:79
double fCharge
Charge.
Definition: RKTrackRep.h:174
Float_t w
Definition: plot.C:20
genf::RKTrackRep::~RKTrackRep ( )
virtual

Definition at line 75 of file RKTrackRep.cxx.

76 {
77  // delete fEffect;
78  //GFMaterialEffects is now a Singleton
79 }
genf::RKTrackRep::RKTrackRep ( const RKTrackRep )
inlineprivate

Definition at line 166 of file RKTrackRep.h.

Member Function Documentation

void genf::GFAbsTrackRep::addChiSqu ( double  aChiSqu)
inlineinherited

Definition at line 264 of file GFAbsTrackRep.h.

Referenced by genf::GFKalman::processHit().

264 { fChiSqu += aChiSqu; }
double fChiSqu
chiSqu of the track fit
Definition: GFAbsTrackRep.h:95
void genf::GFAbsTrackRep::addNDF ( unsigned int  n)
inlineinherited

Definition at line 265 of file GFAbsTrackRep.h.

References n.

Referenced by genf::GFKalman::processHit().

265 { fNdf += n; }
unsigned int fNdf
Definition: GFAbsTrackRep.h:96
Char_t n[5]
virtual GFAbsTrackRep* genf::RKTrackRep::clone ( ) const
inlinevirtual

Implements genf::GFAbsTrackRep.

Definition at line 73 of file RKTrackRep.h.

References RKTrackRep().

73 { return new RKTrackRep(*this); }
double genf::RKTrackRep::Extrap ( const GFDetPlane plane,
TMatrixT< Double_t > *  state,
TMatrixT< Double_t > *  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 1134 of file RKTrackRep.cxx.

References genf::GFDetPlane::dist(), larg4::dist(), genf::GFDetPlane::distance(), E, genf::GFMaterialEffects::effects(), fCharge, trkf::fill(), fPdg, genf::GFMaterialEffects::getInstance(), genf::GFDetPlane::inActive(), MINSTEP, noise(), RKutta(), GFException::setFatal(), and value.

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

1137 {
1138 
1139  static const int maxNumIt(2000);
1140  int numIt(0);
1141  bool calcCov(true);
1142  if (cov == NULL) calcCov = false;
1143 
1144  double P[56]; // only 7 used if no covariance matrix
1145  if (calcCov) std::fill(P, P + std::extent<decltype(P)>::value, 0);
1146  // else {} // not needed std::fill(P, P + 7, 0);
1147 
1148  for (int i = 0; i < 7; ++i) {
1149  P[i] = (*state)[i][0];
1150  }
1151 
1152  TMatrixT<Double_t> jac(7, 7);
1153  TMatrixT<Double_t> jacT(7, 7);
1154  TMatrixT<Double_t> oldCov(7, 7);
1155  if (calcCov) oldCov = (*cov);
1156  double coveredDistance(0.);
1157  double sumDistance(0.);
1158 
1159  while (true) {
1160  if (numIt++ > maxNumIt) {
1161  throw GFException(
1162  "RKTrackRep::Extrap ==> maximum number of iterations exceeded", __LINE__, __FILE__)
1163  .setFatal();
1164  }
1165 
1166  if (calcCov) {
1167  memset(&P[7], 0x00, 49 * sizeof(double));
1168  for (int i = 0; i < 6; ++i) {
1169  P[(i + 1) * 7 + i] = 1.;
1170  }
1171  P[55] = (*state)[6][0];
1172  }
1173 
1174  // double dir(1.);
1175  {
1176  TVector3 Pvect(P[0], P[1], P[2]); //position
1177  TVector3 Avect(P[3], P[4], P[5]); //direction
1178  TVector3 dist = plane.dist(Pvect); //from point to plane
1179  // std::cout << "RKTrackRep::Extrap(): Position, Direction, Plane, dist " << std::endl;
1180  //Pvect.Print();
1181  //Avect.Print();
1182  //plane.Print();
1183  //dist.Print();
1184  // if(dist*Avect<0.) dir=-1.;
1185  }
1186 
1187  TVector3 directionBefore(P[3], P[4], P[5]); // direction before propagation
1188  directionBefore.SetMag(1.);
1189 
1190  // propagation
1191  std::vector<TVector3> points;
1192  std::vector<double> pointPaths;
1193  if (!this->RKutta(plane,
1194  P,
1195  coveredDistance,
1196  points,
1197  pointPaths,
1198  -1.,
1199  calcCov)) { // maxLen currently not used
1200  //GFException exc("RKTrackRep::Extrap ==> Runge Kutta propagation failed",__LINE__,__FILE__);
1201 
1202  //exc.setFatal(); // stops propagation; faster, but some hits will be lost
1203  mf::LogInfo("RKTrackRep::RKutta(): ") << "Throw cet exception here, ... ";
1204  throw GFException("RKTrackRep::RKutta(): Runge Kutta propagation failed", __LINE__, __FILE__)
1205  .setFatal();
1206  }
1207 
1208  TVector3 directionAfter(P[3], P[4], P[5]); // direction after propagation
1209  directionAfter.SetMag(1.);
1210 
1211  sumDistance += coveredDistance;
1212 
1213  // filter Points
1214  std::vector<TVector3> pointsFilt(1, points.at(0));
1215  std::vector<double> pointPathsFilt(1, 0.);
1216  // only if in right direction
1217  for (unsigned int i = 1; i < points.size(); ++i) {
1218  if (pointPaths.at(i) * coveredDistance > 0.) {
1219  pointsFilt.push_back(points.at(i));
1220  pointPathsFilt.push_back(pointPaths.at(i));
1221  }
1222  else {
1223  pointsFilt.back() = points.at(i);
1224  pointPathsFilt.back() += pointPaths.at(i);
1225  }
1226  // clean up tiny steps
1227  int position = pointsFilt.size() - 1; // position starts with 0
1228  if (fabs(pointPathsFilt.back()) < MINSTEP && position > 1) {
1229  pointsFilt.at(position - 1) = pointsFilt.at(position);
1230  pointsFilt.pop_back();
1231  pointPathsFilt.at(position - 1) += pointPathsFilt.at(position);
1232  pointPathsFilt.pop_back();
1233  }
1234  }
1235 
1236  //consistency check
1237  double checkSum(0.);
1238  for (unsigned int i = 0; i < pointPathsFilt.size(); ++i) {
1239  checkSum += pointPathsFilt.at(i);
1240  }
1241  //assert(fabs(checkSum-coveredDistance)<1.E-7);
1242  if (fabs(checkSum - coveredDistance) > 1.E-7) {
1243  throw GFException(
1244  "RKTrackRep::Extrap ==> fabs(checkSum-coveredDistance)>1.E-7", __LINE__, __FILE__)
1245  .setFatal();
1246  }
1247 
1248  if (calcCov) { //calculate Jacobian jac
1249  for (int i = 0; i < 7; ++i) {
1250  for (int j = 0; j < 7; ++j) {
1251  if (i < 6)
1252  jac[i][j] = P[(i + 1) * 7 + j];
1253  else
1254  jac[i][j] = P[(i + 1) * 7 + j] / P[6];
1255  }
1256  }
1257  jacT = jac;
1258  jacT.T();
1259  }
1260 
1261  TMatrixT<Double_t> noise(7, 7); // zero everywhere by default
1262 
1263  // call MatEffects
1264  double momLoss; // momLoss has a sign - negative loss means momentum gain
1265 
1266  /*
1267  momLoss = fEffect->effects(pointsFilt,
1268  pointPathsFilt,
1269  fabs(fCharge/P[6]), // momentum
1270  fPdg,
1271  calcCov,
1272  &noise,
1273  &jac,
1274  &directionBefore,
1275  &directionAfter);
1276  */
1277  momLoss = GFMaterialEffects::getInstance()->effects(pointsFilt,
1278  pointPathsFilt,
1279  fabs(fCharge / P[6]), // momentum
1280  fPdg,
1281  calcCov,
1282  &noise,
1283  &jac,
1284  &directionBefore,
1285  &directionAfter);
1286 
1287  if (fabs(P[6]) > 1.E-10) { // do momLoss only for defined 1/momentum .ne.0
1288  P[6] = fCharge / (fabs(fCharge / P[6]) - momLoss);
1289  }
1290 
1291  if (calcCov) { //propagate cov and add noise
1292  oldCov = *cov;
1293  *cov = jacT * ((oldCov)*jac) + noise;
1294  }
1295 
1296  //we arrived at the destination plane, if we point to the active area
1297  //of the plane (if it is finite), and the distance is below threshold
1298  if (plane.inActive(TVector3(P[0], P[1], P[2]), TVector3(P[3], P[4], P[5]))) {
1299  if (plane.distance(P[0], P[1], P[2]) < MINSTEP) break;
1300  }
1301  }
1302  (*state)[0][0] = P[0];
1303  (*state)[1][0] = P[1];
1304  (*state)[2][0] = P[2];
1305  (*state)[3][0] = P[3];
1306  (*state)[4][0] = P[4];
1307  (*state)[5][0] = P[5];
1308  (*state)[6][0] = P[6];
1309 
1310  return sumDistance;
1311 }
int fPdg
PDG particle code.
Definition: RKTrackRep.h:170
static GFMaterialEffects * getInstance()
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
unsigned int noise()
Definition: chem4.cc:261
Float_t E
Definition: plot.C:20
void fill(const art::PtrVector< recob::Hit > &hits, int only_plane)
#define MINSTEP
Definition: RKTrackRep.cxx:39
double value
Definition: spectrum.C:18
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:47
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
double effects(const std::vector< TVector3 > &points, const std::vector< double > &pointPaths, const double &mom, const int &pdg, const bool &doNoise=false, TMatrixT< Double_t > *noise=NULL, const TMatrixT< Double_t > *jacobian=NULL, const TVector3 *directionBefore=NULL, const TVector3 *directionAfter=NULL)
Calculates energy loss in the travelled path, optional calculation of noise matrix.
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:75
double fCharge
Charge.
Definition: RKTrackRep.h:174
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
Contains all material effects.
Definition: RKTrackRep.cxx:682
double genf::RKTrackRep::extrapolate ( const GFDetPlane pl,
TMatrixT< Double_t > &  statePred,
TMatrixT< Double_t > &  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 genf::GFAbsTrackRep.

Definition at line 468 of file RKTrackRep.cxx.

References E, Extrap(), fCachePlane, fCacheSpu, genf::GFAbsTrackRep::fCov, genf::GFAbsTrackRep::fRefPlane, fSpu, genf::GFAbsTrackRep::fState, genf::GFDetPlane::getNormal(), genf::GFDetPlane::getO(), genf::GFDetPlane::getU(), genf::GFDetPlane::getV(), genf::PrintROOTobject(), rescaleCovOffDiags(), w, X, Y, and Z.

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

471 {
472 
473  TMatrixT<Double_t> cov7x7(7, 7);
474  TMatrixT<Double_t> J_pM(7, 5);
475 
476  TVector3 o = fRefPlane.getO();
477  TVector3 u = fRefPlane.getU();
478  TVector3 v = fRefPlane.getV();
479  TVector3 w = u.Cross(v);
480  std::ostream* pOut = nullptr; // &std::cout if you really really want
481 
482  J_pM[0][3] = u.X();
483  J_pM[0][4] = v.X(); // dx/du
484  J_pM[1][3] = u.Y();
485  J_pM[1][4] = v.Y();
486  J_pM[2][3] = u.Z();
487  J_pM[2][4] = v.Z();
488 
489  TVector3 pTilde = fSpu * (w + fState[1][0] * u + fState[2][0] * v);
490  double pTildeMag = pTilde.Mag();
491 
492  //J_pM matrix is d(x,y,z,ax,ay,az,q/p) / d(q/p,u',v',u,v)
493 
494  // da_x/du'
495  J_pM[3][1] = fSpu / pTildeMag * (u.X() - pTilde.X() / (pTildeMag * pTildeMag) * u * pTilde);
496  J_pM[4][1] = fSpu / pTildeMag * (u.Y() - pTilde.Y() / (pTildeMag * pTildeMag) * u * pTilde);
497  J_pM[5][1] = fSpu / pTildeMag * (u.Z() - pTilde.Z() / (pTildeMag * pTildeMag) * u * pTilde);
498  // da_x/dv'
499  J_pM[3][2] = fSpu / pTildeMag * (v.X() - pTilde.X() / (pTildeMag * pTildeMag) * v * pTilde);
500  J_pM[4][2] = fSpu / pTildeMag * (v.Y() - pTilde.Y() / (pTildeMag * pTildeMag) * v * pTilde);
501  J_pM[5][2] = fSpu / pTildeMag * (v.Z() - pTilde.Z() / (pTildeMag * pTildeMag) * v * pTilde);
502  // dqOp/dqOp
503  J_pM[6][0] = 1.;
504 
505  TMatrixT<Double_t> J_pM_transp(J_pM);
506  J_pM_transp.T();
507 
508  cov7x7 = J_pM * (fCov * J_pM_transp);
509  if (cov7x7[0][0] >= 1000. || cov7x7[0][0] < 1.E-50) {
510  if (pOut) {
511  (*pOut) << "RKTrackRep::extrapolate(): cov7x7[0][0] is crazy. Rescale off-diags. Try again. "
512  "fCov, cov7x7 were: "
513  << std::endl;
514  PrintROOTobject(*pOut, fCov);
515  PrintROOTobject(*pOut, cov7x7);
516  }
518  cov7x7 = J_pM * (fCov * J_pM_transp);
519  if (pOut) {
520  (*pOut) << "New cov7x7 and fCov are ... " << std::endl;
521  PrintROOTobject(*pOut, cov7x7);
522  PrintROOTobject(*pOut, fCov);
523  }
524  }
525 
526  TVector3 pos = o + fState[3][0] * u + fState[4][0] * v;
527  TMatrixT<Double_t> state7(7, 1);
528  state7[0][0] = pos.X();
529  state7[1][0] = pos.Y();
530  state7[2][0] = pos.Z();
531  state7[3][0] = pTilde.X() / pTildeMag;
532  ;
533  state7[4][0] = pTilde.Y() / pTildeMag;
534  ;
535  state7[5][0] = pTilde.Z() / pTildeMag;
536  ;
537  state7[6][0] = fState[0][0];
538 
539  double coveredDistance = this->Extrap(pl, &state7, &cov7x7);
540 
541  TVector3 O = pl.getO();
542  TVector3 U = pl.getU();
543  TVector3 V = pl.getV();
544  TVector3 W = pl.getNormal();
545 
546  double X = state7[0][0];
547  double Y = state7[1][0];
548  double Z = state7[2][0];
549  double AX = state7[3][0];
550  double AY = state7[4][0];
551  double AZ = state7[5][0];
552  double QOP = state7[6][0];
553  TVector3 A(AX, AY, AZ);
554  TVector3 Point(X, Y, Z);
555  TMatrixT<Double_t> J_Mp(5, 7);
556 
557  // J_Mp matrix is d(q/p,u',v',u,v) / d(x,y,z,ax,ay,az,q/p)
558  J_Mp[0][6] = 1.;
559  //du'/da_x
560  double AtW = A * W;
561  J_Mp[1][3] = (U.X() * (AtW)-W.X() * (A * U)) / (AtW * AtW);
562  J_Mp[1][4] = (U.Y() * (AtW)-W.Y() * (A * U)) / (AtW * AtW);
563  J_Mp[1][5] = (U.Z() * (AtW)-W.Z() * (A * U)) / (AtW * AtW);
564  //dv'/da_x
565  J_Mp[2][3] = (V.X() * (AtW)-W.X() * (A * V)) / (AtW * AtW);
566  J_Mp[2][4] = (V.Y() * (AtW)-W.Y() * (A * V)) / (AtW * AtW);
567  J_Mp[2][5] = (V.Z() * (AtW)-W.Z() * (A * V)) / (AtW * AtW);
568  //du/dx
569  J_Mp[3][0] = U.X();
570  J_Mp[3][1] = U.Y();
571  J_Mp[3][2] = U.Z();
572  //dv/dx
573  J_Mp[4][0] = V.X();
574  J_Mp[4][1] = V.Y();
575  J_Mp[4][2] = V.Z();
576 
577  TMatrixT<Double_t> J_Mp_transp(J_Mp);
578  J_Mp_transp.T();
579 
580  covPred.ResizeTo(5, 5);
581  covPred = J_Mp * (cov7x7 * J_Mp_transp);
582 
583  statePred.ResizeTo(5, 1);
584  statePred[0][0] = QOP;
585  statePred[1][0] = (A * U) / (A * W);
586  statePred[2][0] = (A * V) / (A * W);
587  statePred[3][0] = (Point - O) * U;
588  statePred[4][0] = (Point - O) * V;
589 
590  fCachePlane = pl;
591  fCacheSpu = (A * W) / fabs(A * W);
592 
593  return coveredDistance;
594 }
Float_t Y
Definition: plot.C:37
GFDetPlane fCachePlane
Definition: RKTrackRep.h:160
Float_t Z
Definition: plot.C:37
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:92
TVector3 getO() const
Definition: GFDetPlane.h:77
Float_t E
Definition: plot.C:20
void rescaleCovOffDiags()
TVector3 getU() const
Definition: GFDetPlane.h:78
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:42
double Extrap(const GFDetPlane &plane, TMatrixT< Double_t > *state, TMatrixT< Double_t > *cov=NULL) const
Handles propagation and material effects.
void PrintROOTobject(std::ostream &, const ROOTOBJ &)
Small utility functions which print some ROOT objects into an output stream.
Definition: GFException.h:127
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
TVector3 getV() const
Definition: GFDetPlane.h:79
Float_t X
Definition: plot.C:37
Float_t w
Definition: plot.C:20
double genf::RKTrackRep::extrapolate ( const GFDetPlane pl,
TMatrixT< Double_t > &  statePred 
)
virtual

returns the tracklength spanned in this extrapolation

Reimplemented from genf::GFAbsTrackRep.

Definition at line 596 of file RKTrackRep.cxx.

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

597 {
598 
599  TVector3 o = fRefPlane.getO();
600  TVector3 u = fRefPlane.getU();
601  TVector3 v = fRefPlane.getV();
602  TVector3 w = u.Cross(v);
603 
604  TVector3 pTilde = fSpu * (w + fState[1][0] * u + fState[2][0] * v);
605  double pTildeMag = pTilde.Mag();
606 
607  TVector3 pos = o + fState[3][0] * u + fState[4][0] * v;
608 
609  TMatrixT<Double_t> state7(7, 1);
610  state7[0][0] = pos.X();
611  state7[1][0] = pos.Y();
612  state7[2][0] = pos.Z();
613  state7[3][0] = pTilde.X() / pTildeMag;
614  state7[4][0] = pTilde.Y() / pTildeMag;
615  state7[5][0] = pTilde.Z() / pTildeMag;
616  state7[6][0] = fState[0][0];
617 
618  TVector3 O = pl.getO();
619  TVector3 U = pl.getU();
620  TVector3 V = pl.getV();
621  TVector3 W = pl.getNormal();
622 
623  double coveredDistance = this->Extrap(pl, &state7);
624 
625  double X = state7[0][0];
626  double Y = state7[1][0];
627  double Z = state7[2][0];
628  double AX = state7[3][0];
629  double AY = state7[4][0];
630  double AZ = state7[5][0];
631  double QOP = state7[6][0];
632  TVector3 A(AX, AY, AZ);
633  TVector3 Point(X, Y, Z);
634 
635  statePred.ResizeTo(5, 1);
636  statePred[0][0] = QOP;
637  statePred[1][0] = (A * U) / (A * W);
638  statePred[2][0] = (A * V) / (A * W);
639  statePred[3][0] = (Point - O) * U;
640  statePred[4][0] = (Point - O) * V;
641 
642  return coveredDistance;
643 }
Float_t Y
Definition: plot.C:37
Float_t Z
Definition: plot.C:37
TVector3 getO() const
Definition: GFDetPlane.h:77
TVector3 getU() const
Definition: GFDetPlane.h:78
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:42
double Extrap(const GFDetPlane &plane, TMatrixT< Double_t > *state, TMatrixT< Double_t > *cov=NULL) const
Handles propagation and material effects.
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
TVector3 getV() const
Definition: GFDetPlane.h:79
Float_t X
Definition: plot.C:37
Float_t w
Definition: plot.C:20
double genf::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 53 of file GFAbsTrackRep.cxx.

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

54 {
55  TMatrixT<Double_t> statePred(fDimension, 1);
56  TMatrixT<Double_t> covPred(fDimension, fDimension);
57  double retVal = extrapolate(plane, statePred, covPred);
58  setData(statePred, plane, &covPred);
59  return retVal;
60 }
unsigned int fDimension
Dimensionality of track representation.
Definition: GFAbsTrackRep.h:86
virtual void setData(const TMatrixT< Double_t > &st, const GFDetPlane &pl, const TMatrixT< Double_t > *cov=NULL)
virtual double extrapolate(const GFDetPlane &plane, TMatrixT< Double_t > &statePred)
returns the tracklength spanned in this extrapolation
void genf::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 genf::GFAbsTrackRep.

Definition at line 415 of file RKTrackRep.cxx.

References Extrap(), genf::GFAbsTrackRep::fRefPlane, fSpu, genf::GFAbsTrackRep::fState, genf::GFDetPlane::getO(), genf::GFDetPlane::getU(), genf::GFDetPlane::getV(), MINSTEP, poca2Line(), GFException::setFatal(), genf::GFDetPlane::setO(), genf::GFDetPlane::setU(), genf::GFDetPlane::setV(), and w.

Referenced by prototype().

420 {
421  static const int maxIt(30);
422 
423  TVector3 o = fRefPlane.getO();
424  TVector3 u = fRefPlane.getU();
425  TVector3 v = fRefPlane.getV();
426  TVector3 w = u.Cross(v);
427 
428  TVector3 pTilde = fSpu * (w + fState[1][0] * u + fState[2][0] * v);
429  pTilde.SetMag(1.);
430 
431  TVector3 point = o + fState[3][0] * u + fState[4][0] * v;
432 
433  TMatrixT<Double_t> state7(7, 1);
434  state7[0][0] = point.X();
435  state7[1][0] = point.Y();
436  state7[2][0] = point.Z();
437  state7[3][0] = pTilde.X();
438  state7[4][0] = pTilde.Y();
439  state7[5][0] = pTilde.Z();
440  state7[6][0] = fState[0][0];
441 
442  double coveredDistance(0.);
443 
444  GFDetPlane pl;
445  int iterations(0);
446 
447  while (true) {
448  pl.setO(point1);
449  TVector3 currentDir(state7[3][0], state7[4][0], state7[5][0]);
450  pl.setU(currentDir.Cross(point2 - point1));
451  pl.setV(point2 - point1);
452  coveredDistance = this->Extrap(pl, &state7);
453 
454  if (fabs(coveredDistance) < MINSTEP) break;
455  if (++iterations == maxIt) {
456  throw GFException(
457  "RKTrackRep extrapolation to point failed, maximum number of iterations reached",
458  __LINE__,
459  __FILE__)
460  .setFatal();
461  }
462  }
463  poca.SetXYZ(state7[0][0], state7[1][0], state7[2][0]);
464  dirInPoca.SetXYZ(state7[3][0], state7[4][0], state7[5][0]);
465  poca_onwire = poca2Line(point1, point2, poca);
466 }
TVector3 getO() const
Definition: GFDetPlane.h:77
#define MINSTEP
Definition: RKTrackRep.cxx:39
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:47
TVector3 getU() const
Definition: GFDetPlane.h:78
double Extrap(const GFDetPlane &plane, TMatrixT< Double_t > *state, TMatrixT< Double_t > *cov=NULL) const
Handles propagation and material effects.
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:75
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
TVector3 getV() const
Definition: GFDetPlane.h:79
Float_t w
Definition: plot.C:20
TVector3 poca2Line(const TVector3 &extr1, const TVector3 &extr2, const TVector3 &point) const
Definition: RKTrackRep.cxx:398
void genf::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 genf::GFAbsTrackRep.

Definition at line 352 of file RKTrackRep.cxx.

References Extrap(), genf::GFAbsTrackRep::fRefPlane, fSpu, genf::GFAbsTrackRep::fState, genf::GFDetPlane::getO(), genf::GFDetPlane::getU(), genf::GFDetPlane::getV(), MINSTEP, GFException::setFatal(), genf::GFDetPlane::setON(), and w.

Referenced by prototype().

353 {
354 
355  static const int maxIt(30);
356 
357  TVector3 o = fRefPlane.getO();
358  TVector3 u = fRefPlane.getU();
359  TVector3 v = fRefPlane.getV();
360  TVector3 w = u.Cross(v);
361 
362  TVector3 pTilde = fSpu * (w + fState[1][0] * u + fState[2][0] * v);
363  pTilde.SetMag(1.);
364 
365  TVector3 point = o + fState[3][0] * u + fState[4][0] * v;
366 
367  TMatrixT<Double_t> state7(7, 1);
368  state7[0][0] = point.X();
369  state7[1][0] = point.Y();
370  state7[2][0] = point.Z();
371  state7[3][0] = pTilde.X();
372  state7[4][0] = pTilde.Y();
373  state7[5][0] = pTilde.Z();
374  state7[6][0] = fState[0][0];
375 
376  double coveredDistance(0.);
377 
378  GFDetPlane pl;
379  int iterations(0);
380 
381  while (true) {
382  pl.setON(pos, TVector3(state7[3][0], state7[4][0], state7[5][0]));
383  coveredDistance = this->Extrap(pl, &state7);
384 
385  if (fabs(coveredDistance) < MINSTEP) break;
386  if (++iterations == maxIt) {
387  throw GFException("RKTrackRep::extrapolateToPoint==> extrapolation to point failed, maximum "
388  "number of iterations reached",
389  __LINE__,
390  __FILE__)
391  .setFatal();
392  }
393  }
394  poca.SetXYZ(state7[0][0], state7[1][0], state7[2][0]);
395  dirInPoca.SetXYZ(state7[3][0], state7[4][0], state7[5][0]);
396 }
TVector3 getO() const
Definition: GFDetPlane.h:77
#define MINSTEP
Definition: RKTrackRep.cxx:39
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:47
TVector3 getU() const
Definition: GFDetPlane.h:78
double Extrap(const GFDetPlane &plane, TMatrixT< Double_t > *state, TMatrixT< Double_t > *cov=NULL) const
Handles propagation and material effects.
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:75
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
TVector3 getV() const
Definition: GFDetPlane.h:79
Float_t w
Definition: plot.C:20
const TMatrixT< double > * genf::RKTrackRep::getAuxInfo ( const GFDetPlane pl)

Definition at line 61 of file RKTrackRep.cxx.

References fAuxInfo, fCachePlane, fCacheSpu, and GFException::setFatal().

Referenced by switchDirection().

62 {
63 
64  if (pl != fCachePlane) {
65  throw GFException("RKTrackRep::getAuxInfo() trying to get auxillary information with planes "
66  "mismatch (Information returned does not belong to requested plane)",
67  __LINE__,
68  __FILE__)
69  .setFatal();
70  }
71  fAuxInfo.ResizeTo(1, 1);
72  fAuxInfo(0, 0) = fCacheSpu;
73  return &fAuxInfo;
74 }
GFDetPlane fCachePlane
Definition: RKTrackRep.h:160
TMatrixT< double > fAuxInfo
Definition: RKTrackRep.h:163
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:47
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:75
double genf::RKTrackRep::getCharge ( ) const
inlinevirtual

Returns charge.

Implements genf::GFAbsTrackRep.

Definition at line 133 of file RKTrackRep.h.

References fCharge.

133 { return fCharge; }
double fCharge
Charge.
Definition: RKTrackRep.h:174
double genf::GFAbsTrackRep::getChiSqu ( ) const
inlineinherited
const TMatrixT<Double_t>& genf::GFAbsTrackRep::getCov ( ) const
inlineinherited
double genf::GFAbsTrackRep::getCovElem ( int  i,
int  j 
) const
inlineinherited
unsigned int genf::GFAbsTrackRep::getDim ( ) const
inlineinherited

returns dimension of state vector

Definition at line 182 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fDimension, and genf::GFAbsTrackRep::Print().

Referenced by genf::GFKalman::getChi2Hit(), genf::GFAbsTrackRep::getNDF(), genf::GFTrack::getResiduals(), genf::GFKalman::processHit(), and genf::GFDaf::processTrack().

182 { return fDimension; }
unsigned int fDimension
Dimensionality of track representation.
Definition: GFAbsTrackRep.h:86
TMatrixT<Double_t> genf::GFAbsTrackRep::getFirstCov ( ) const
inlineinherited

Definition at line 218 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fFirstCov.

218 { return fFirstCov; }
TMatrixT< Double_t > fFirstCov
GFDetPlane genf::GFAbsTrackRep::getFirstPlane ( ) const
inlineinherited

Definition at line 219 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fFirstPlane.

219 { return fFirstPlane; }
GFDetPlane fFirstPlane
TMatrixT<Double_t> genf::GFAbsTrackRep::getFirstState ( ) const
inlineinherited

Definition at line 217 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fFirstState.

217 { return fFirstState; }
TMatrixT< Double_t > fFirstState
state, cov and plane for first and last point in fit
TMatrixT<Double_t> genf::GFAbsTrackRep::getLastCov ( ) const
inlineinherited

Definition at line 221 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fLastCov.

221 { return fLastCov; }
TMatrixT< Double_t > fLastCov
GFDetPlane genf::GFAbsTrackRep::getLastPlane ( ) const
inlineinherited

Definition at line 222 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fLastPlane.

Referenced by trkf::Track3DKalman::produce(), and trkf::Track3DKalmanSPS::produce().

222 { return fLastPlane; }
TMatrixT<Double_t> genf::GFAbsTrackRep::getLastState ( ) const
inlineinherited

Definition at line 220 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fLastState.

220 { return fLastState; }
TMatrixT< Double_t > fLastState
TVector3 genf::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 genf::GFAbsTrackRep.

Definition at line 312 of file RKTrackRep.cxx.

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

313 {
314  TMatrixT<Double_t> statePred(fState);
315  TVector3 retmom;
316  if (pl != fRefPlane) {
317  extrapolate(pl, statePred);
318  retmom =
319  fCacheSpu * (pl.getNormal() + statePred[1][0] * pl.getU() + statePred[2][0] * pl.getV());
320  }
321  else {
322  retmom = fSpu * (pl.getNormal() + statePred[1][0] * pl.getU() + statePred[2][0] * pl.getV());
323  }
324  retmom.SetMag(1. / fabs(statePred[0][0]));
325  return retmom;
326 }
double extrapolate(const GFDetPlane &, TMatrixT< Double_t > &statePred, TMatrixT< Double_t > &covPred)
returns the tracklength spanned in this extrapolation
Definition: RKTrackRep.cxx:468
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
TVector3 genf::GFAbsTrackRep::getMom ( )
inlineinherited
TVector3 genf::RKTrackRep::getMomLast ( const GFDetPlane pl)

Definition at line 328 of file RKTrackRep.cxx.

References genf::GFAbsTrackRep::fLastState, fSpu, genf::GFDetPlane::getNormal(), genf::GFDetPlane::getU(), and genf::GFDetPlane::getV().

Referenced by prototype().

329 {
330  TMatrixT<Double_t> statePred(fLastState);
331  TVector3 retmom;
332  retmom = fSpu * (pl.getNormal() + statePred[1][0] * pl.getU() + statePred[2][0] * pl.getV());
333 
334  retmom.SetMag(1. / fabs(statePred[0][0]));
335  return retmom;
336 }
TMatrixT< Double_t > fLastState
unsigned int genf::GFAbsTrackRep::getNDF ( ) const
inlineinherited

Definition at line 230 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::getDim().

Referenced by genf::GFTrack::getNDF(), genf::GFAbsTrackRep::getRedChiSqu(), trkf::Track3DKalman::produce(), and trkf::Track3DKalmanSPS::produce().

231  {
232  if (fNdf > getDim()) return fNdf - getDim();
233  return 0;
234  }
unsigned int fNdf
Definition: GFAbsTrackRep.h:96
unsigned int getDim() const
returns dimension of state vector
int genf::RKTrackRep::getPDG ( )

Definition at line 297 of file RKTrackRep.cxx.

References fPdg.

Referenced by switchDirection().

298 {
299  return fPdg;
300 }
int fPdg
PDG particle code.
Definition: RKTrackRep.h:170
TVector3 genf::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 genf::GFAbsTrackRep.

Definition at line 302 of file RKTrackRep.cxx.

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

303 {
304  if (pl != fRefPlane) {
305  TMatrixT<Double_t> s(5, 1);
306  extrapolate(pl, s);
307  return pl.getO() + s[3][0] * pl.getU() + s[4][0] * pl.getV();
308  }
309  return fRefPlane.getO() + fState[3][0] * fRefPlane.getU() + fState[4][0] * fRefPlane.getV();
310 }
double extrapolate(const GFDetPlane &, TMatrixT< Double_t > &statePred, TMatrixT< Double_t > &covPred)
returns the tracklength spanned in this extrapolation
Definition: RKTrackRep.cxx:468
TVector3 getO() const
Definition: GFDetPlane.h:77
TVector3 getU() const
Definition: GFDetPlane.h:78
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
TVector3 getV() const
Definition: GFDetPlane.h:79
TVector3 genf::GFAbsTrackRep::getPos ( )
inlineinherited
void genf::RKTrackRep::getPosMom ( const GFDetPlane pl,
TVector3 &  pos,
TVector3 &  mom 
)
virtual

Gets position and momentum in the plane by exrapolating or not.

If #GFDetPlane equals the reference plane fRefPlane, it gets current position and momentum; otherwise for getMom() it extrapolates the track to the plane and gets the position and momentum. GetMomLast() does no extrapn. EC.

Implements genf::GFAbsTrackRep.

Definition at line 338 of file RKTrackRep.cxx.

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

Referenced by prototype().

339 {
340  TMatrixT<Double_t> statePred(fState);
341  if (pl != fRefPlane) {
342  extrapolate(pl, statePred);
343  mom = fCacheSpu * (pl.getNormal() + statePred[1][0] * pl.getU() + statePred[2][0] * pl.getV());
344  }
345  else {
346  mom = fSpu * (pl.getNormal() + statePred[1][0] * pl.getU() + statePred[2][0] * pl.getV());
347  }
348  mom.SetMag(1. / fabs(statePred[0][0]));
349  pos = pl.getO() + (statePred[3][0] * pl.getU()) + (statePred[4][0] * pl.getV());
350 }
double extrapolate(const GFDetPlane &, TMatrixT< Double_t > &statePred, TMatrixT< Double_t > &covPred)
returns the tracklength spanned in this extrapolation
Definition: RKTrackRep.cxx:468
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
void genf::GFAbsTrackRep::getPosMomCov ( const GFDetPlane pl,
TVector3 &  pos,
TVector3 &  mom,
TMatrixT< Double_t > &  cov 
)
virtualinherited

method which gets position, momentum and 6x6 covariance matrix

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

Definition at line 101 of file GFAbsTrackRep.cxx.

References genf::GFAbsTrackRep::Abort().

Referenced by genf::GFAbsTrackRep::getCovElem(), genf::GFAbsTrackRep::getPosMomCov(), and genf::GFTrack::getPosMomCov().

105 {
106  Abort("getPosMomCov()");
107 }
void Abort(std::string method)
void genf::GFAbsTrackRep::getPosMomCov ( TVector3 &  pos,
TVector3 &  mom,
TMatrixT< Double_t > &  c 
)
inlineinherited

Definition at line 212 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::getPosMomCov().

213  {
214  getPosMomCov(fRefPlane, pos, mom, c);
215  }
virtual void getPosMomCov(const GFDetPlane &pl, TVector3 &pos, TVector3 &mom, TMatrixT< Double_t > &cov)
method which gets position, momentum and 6x6 covariance matrix
double genf::GFAbsTrackRep::getRedChiSqu ( ) const
inlineinherited

returns chi2/ndf

Definition at line 225 of file GFAbsTrackRep.h.

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

Referenced by genf::GFTrack::getRedChiSqu().

226  {
227  if (getNDF() > 0) return getChiSqu() / getNDF();
228  return 0;
229  }
unsigned int getNDF() const
double getChiSqu() const
const GFDetPlane& genf::GFAbsTrackRep::getReferencePlane ( ) const
inlineinherited
const TMatrixT<Double_t>& genf::GFAbsTrackRep::getState ( ) const
inlineinherited
double genf::GFAbsTrackRep::getStateElem ( int  i) const
inlineinherited

Definition at line 189 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fState.

189 { return fState(i, 0); }
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
bool genf::GFAbsTrackRep::getStatusFlag ( )
inlineinherited
bool genf::RKTrackRep::hasAuxInfo ( )
inline

Definition at line 157 of file RKTrackRep.h.

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

Definition at line 165 of file RKTrackRep.h.

165 { return *this; }
TVector3 genf::RKTrackRep::poca2Line ( const TVector3 &  extr1,
const TVector3 &  extr2,
const TVector3 &  point 
) const
private

Definition at line 398 of file RKTrackRep.cxx.

References GFException::setFatal().

Referenced by extrapolateToLine().

401 {
402 
403  TVector3 theWire = extr2 - extr1;
404  if (theWire.Mag() < 1.E-8) {
405  throw GFException("RKTrackRep::poca2Line(): try to find poca between line and point, but the "
406  "line is really just a point",
407  __LINE__,
408  __FILE__)
409  .setFatal();
410  }
411  double t = 1. / (theWire * theWire) * (point * theWire + extr1 * extr1 - extr1 * extr2);
412  return (extr1 + t * theWire);
413 }
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:47
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:75
void genf::GFAbsTrackRep::Print ( std::ostream &  out = std::cout) const
virtualinherited

Definition at line 122 of file GFAbsTrackRep.cxx.

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

Referenced by genf::GFAbsTrackRep::getDim(), and genf::GFTrack::Print().

123 {
124  out << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
125  out << "GFAbsTrackRep::Parameters at reference plane ";
126  fRefPlane.Print(out);
127  out << "GFAbsTrackRep::State" << std::endl;
128  PrintROOTmatrix(out, fState);
129  out << "GFAbsTrackRep::Covariances" << std::endl;
130  PrintROOTmatrix(out, fCov);
131  out << "GFAbsTrackRep::chi^2" << std::endl;
132  out << fChiSqu << std::endl;
133  out << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
134 }
double fChiSqu
chiSqu of the track fit
Definition: GFAbsTrackRep.h:95
void PrintROOTmatrix(std::ostream &out, const TMatrixT< T > &m)
Small utility functions which print some ROOT objects into an output stream.
Definition: GFException.h:133
void Print(std::ostream &out=std::cout) const
Definition: GFDetPlane.cxx:228
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:92
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
virtual GFAbsTrackRep* genf::RKTrackRep::prototype ( ) const
inlinevirtual
void genf::RKTrackRep::rescaleCovOffDiags ( )

Definition at line 1313 of file RKTrackRep.cxx.

References genf::GFAbsTrackRep::fCov.

Referenced by extrapolate(), and switchDirection().

1314 {
1315 
1316  for (int i = 0; i < fCov.GetNrows(); ++i) {
1317  for (int j = 0; j < fCov.GetNcols(); ++j) {
1318  if (i != j) { //off diagonal
1319  fCov[i][j] = 0.5 * fCov[i][j];
1320  }
1321  else { //diagonal. Do not let diag cov element be <=0.
1322  if (fCov[i][j] <= 0.0) fCov[i][j] = 0.01;
1323  }
1324  }
1325  }
1326 }
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:92
void genf::GFAbsTrackRep::reset ( )
virtualinherited

Definition at line 109 of file GFAbsTrackRep.cxx.

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

Referenced by genf::GFAbsTrackRep::getStatusFlag().

110 {
111  std::cout << "GFAbsTrackRep::reset" << std::endl;
112  TVector3 nullVec(0., 0., 0.);
113  fRefPlane.set(nullVec, nullVec, nullVec);
114  fState.Zero();
115  fCov.Zero();
116  fFirstState.Zero();
117  fFirstCov.Zero();
118  fLastState.Zero();
119  fLastCov.Zero();
120 }
void set(const TVector3 &o, const TVector3 &u, const TVector3 &v)
Definition: GFDetPlane.cxx:86
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:92
TMatrixT< Double_t > fFirstCov
TMatrixT< Double_t > fLastState
TMatrixT< Double_t > fLastCov
TMatrixT< Double_t > fFirstState
state, cov and plane for first and last point in fit
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
bool genf::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

Contains all material effects.

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 682 of file RKTrackRep.cxx.

References genf::GFDetPlane::distance(), fCharge, fPdg, genf::GFFieldManager::getFieldVal(), genf::GFMaterialEffects::getInstance(), genf::GFDetPlane::getNormal(), genf::GFDetPlane::getO(), genf::GFDetPlane::inActive(), MINSTEP, r, GFException::setFatal(), and genf::GFMaterialEffects::stepper().

Referenced by Extrap().

689 {
690 
691  static const double EC = .000149896229; // c/(2*10^12) resp. c/2Tera
692  static const double DLT = .0002; // max. deviation for approximation-quality test
693  static const double DLT32 = DLT / 32.; //
694  static const double P3 = 1. / 3.; // 1/3
695  static const double Smax = 100.0; // max. step allowed 100->.4 EC, 6-Jan-2-11, 100->1.0
696  static const double Wmax = 3000.; // max. way allowed.
697  //static const double Pmin = 4.E-3; // minimum momentum for propagation [GeV]
698  static const double Pmin = 1.E-11; // minimum momentum for propagation [GeV]
699 
700  static const int ND = 56; // number of variables for derivatives calculation
701  static const int ND1 = ND - 7; // = 49
702  double* R = &P[0]; // Start coordinates in cm ( x, y, z)
703  double* A = &P[3]; // Start directions (ax, ay, az); ax^2+ay^2+az^2=1
704  double SA[3] = {0., 0., 0.}; // Start directions derivatives
705  double Pinv = P[6] * EC; // P[6] is charge/momentum in e/(Gev/c)
706  double Way = 0.; // Total way of the trajectory
707  double Way2 = 0.; // Total way of the trajectory with correct signs
708  bool error = false; // Error of propogation
709  bool stopBecauseOfMaterial =
710  false; // does not go through main loop again when stepsize is reduced by stepper
711 
712  points.clear();
713  pointPaths.clear();
714 
715  if (fabs(fCharge / P[6]) < Pmin) {
716  std::cerr << "RKTrackRep::RKutta ==> momentum too low: " << fabs(fCharge / P[6]) * 1000.
717  << " MeV" << std::endl;
718  return (false);
719  }
720 
721  double SU[4];
722  TVector3 O = plane.getO();
723  TVector3 W = plane.getNormal();
724  if (W * O > 0) { // make SU vector point away from origin
725  SU[0] = W.X();
726  SU[1] = W.Y();
727  SU[2] = W.Z();
728  }
729  else {
730  SU[0] = -1. * W.X();
731  SU[1] = -1. * W.Y();
732  SU[2] = -1. * W.Z();
733  }
734  SU[3] = plane.distance(0., 0., 0.);
735 
736  //
737  // Step estimation until surface
738  //
739  double Step, An, Dist, Dis, S, Sl = 0;
740 
741  points.push_back(TVector3(R[0], R[1], R[2]));
742  pointPaths.push_back(0.);
743 
744  An = A[0] * SU[0] + A[1] * SU[1] + A[2] * SU[2]; // An = A * N; component of A normal to surface
745 
746  if (An == 0) {
747  // std::cerr<<"PaAlgo::RKutta ==> Component Normal to surface is 0!: "<<Step<<" cm !"<<std::endl;
748  std::cerr << "RKTrackRep::RKutta ==> cannot propagate perpendicular to plane " << std::endl;
749  return false; // no propagation if A perpendicular to surface}
750  }
751  if (plane.inActive(
752  TVector3(R[0], R[1], R[2]),
753  TVector3(A[0], A[1], A[2]))) { // if direction is not pointing to active part of surface
754  Dist = SU[3] - R[0] * SU[0] - R[1] * SU[1] -
755  R[2] * SU[2]; // Distance between start coordinates and surface
756  Step = Dist / An;
757  }
758  else { // make sure to extrapolate towards surface
759  if ((O.X() - R[0]) * A[0] + (O.Y() - R[1]) * A[1] + (O.Z() - R[2]) * A[2] >
760  0) { // if direction A pointing from start coordinates R towards surface
761  Dist =
762  sqrt((R[0] - O.X()) *
763  (R[0] - O.X()) + // |R-O|; Distance between start coordinates and origin of surface
764  (R[1] - O.Y()) * (R[1] - O.Y()) +
765  (R[2] - O.Z()) * (R[2] - O.Z()));
766  }
767  else { // if direction pointing away from surface
768  Dist = -1. * sqrt((R[0] - O.X()) * (R[0] - O.X()) + (R[1] - O.Y()) * (R[1] - O.Y()) +
769  (R[2] - O.Z()) * (R[2] - O.Z()));
770  }
771  Step = Dist;
772  }
773 
774  if (fabs(Step) > Wmax) {
775  // std::cerr<<"PaAlgo::RKutta ==> Too long extrapolation requested : "<<Step<<" cm !"<<std::endl;
776  std::cerr << "RKTrackRep::RKutta ==> Too long extrapolation requested : " << Step << " cm !"
777  << std::endl;
778  std::cerr << "X = " << R[0] << " Y = " << R[1] << " Z = " << R[2] << " COSx = " << A[0]
779  << " COSy = " << A[1] << " COSz = " << A[2] << std::endl;
780  std::cout << "Destination X = " << SU[0] * SU[3] << std::endl;
781 
782  mf::LogInfo("RKTrackRep::RKutta(): ") << "Throw cet exception here, ... ";
783  throw GFException("RKTrackRep::RKutta(): Runge Kutta propagation failed", __LINE__, __FILE__)
784  .setFatal();
785 
786  return (false);
787 
788  //std::cout << "RKUtta.EC: But keep plowing on, lowering Step"<< std::endl;
789  }
790 
791  // reduce maximum stepsize S to Smax
792  Step > Smax ? S = Smax : Step < -Smax ? S = -Smax : S = Step;
793 
794  //
795  // Main cycle of Runge-Kutta method
796  //
797 
798  //for saving points only when the direction didnt change
799  int Ssign = 1;
800  if (S < 0) Ssign = -1;
801 
802  while (fabs(Step) > MINSTEP && !stopBecauseOfMaterial) {
803 
804  // call stepper and reduce stepsize
805  double stepperLen;
806  //std::cout<< "RKTrackRep: About to enter fEffect->stepper()" << std::endl;
807  /*
808  stepperLen = fEffect->stepper(fabs(S),
809  R[0],R[1],R[2],
810  Ssign*A[0],Ssign*A[1],Ssign*A[2],
811  fabs(fCharge/P[6]),
812  fPdg);
813  */
814  // From Genfit svn, 27-Sep-2011.
815  stepperLen = GFMaterialEffects::getInstance()->stepper(fabs(S),
816  R[0],
817  R[1],
818  R[2],
819  Ssign * A[0],
820  Ssign * A[1],
821  Ssign * A[2],
822  fabs(fCharge / P[6]),
823  fPdg);
824 
825  //std::cout<< "RKTrackRep: S,R[0],R[1],R[2], A[0],A[1],A[2], stepperLen is " << S <<", "<< R[0] <<", "<< R[1]<<", " << R[2] <<", "<< A[0]<<", " << A[1]<<", " << A[2]<<", " << stepperLen << std::endl;
826  if (stepperLen < MINSTEP)
827  stepperLen = MINSTEP; // prevents tiny stepsizes that can slow down the program
828  if (S > stepperLen) {
829  S = stepperLen;
830  stopBecauseOfMaterial = true;
831  }
832  else if (S < -stepperLen) {
833  S = -stepperLen;
834  stopBecauseOfMaterial = true;
835  }
836 
837  double H0[12], H1[12], H2[12], r[3];
838  double S3 = P3 * S, S4 = .25 * S, PS2 = Pinv * S;
839 
840  //
841  // First point
842  //
843  r[0] = R[0];
844  r[1] = R[1];
845  r[2] = R[2];
846  TVector3 pos(r[0], r[1], r[2]); // vector of start coordinates R0 (x, y, z)
847  TVector3 H0vect = GFFieldManager::getFieldVal(pos); // magnetic field in 10^-4 T = kGauss
848  H0[0] = PS2 * H0vect.X();
849  H0[1] = PS2 * H0vect.Y();
850  H0[2] = PS2 * H0vect.Z(); // H0 is PS2*(Hx, Hy, Hz) @ R0
851  double A0 = A[1] * H0[2] - A[2] * H0[1], B0 = A[2] * H0[0] - A[0] * H0[2],
852  C0 = A[0] * H0[1] - A[1] * H0[0]; // (ax, ay, az) x H0
853  double A2 = A[0] + A0, B2 = A[1] + B0, C2 = A[2] + C0; // (A0, B0, C0) + (ax, ay, az)
854  double A1 = A2 + A[0], B1 = B2 + A[1], C1 = C2 + A[2]; // (A0, B0, C0) + 2*(ax, ay, az)
855 
856  //
857  // Second point
858  //
859  r[0] += A1 * S4;
860  r[1] += B1 * S4;
861  r[2] += C1 * S4; //setup.Field(r,H1);
862  pos.SetXYZ(r[0], r[1], r[2]);
863  TVector3 H1vect = GFFieldManager::getFieldVal(pos);
864  H1[0] = H1vect.X() * PS2;
865  H1[1] = H1vect.Y() * PS2;
866  H1[2] = H1vect.Z() *
867  PS2; // H1 is PS2*(Hx, Hy, Hz) @ (x, y, z) + 0.25*S * [(A0, B0, C0) + 2*(ax, ay, az)]
868  double A3, B3, C3, A4, B4, C4, A5, B5, C5;
869  A3 = B2 * H1[2] - C2 * H1[1] + A[0];
870  B3 = C2 * H1[0] - A2 * H1[2] + A[1];
871  C3 = A2 * H1[1] - B2 * H1[0] + A[2]; // (A2, B2, C2) x H1 + (ax, ay, az)
872  A4 = B3 * H1[2] - C3 * H1[1] + A[0];
873  B4 = C3 * H1[0] - A3 * H1[2] + A[1];
874  C4 = A3 * H1[1] - B3 * H1[0] + A[2]; // (A3, B3, C3) x H1 + (ax, ay, az)
875  A5 = A4 - A[0] + A4;
876  B5 = B4 - A[1] + B4;
877  C5 = C4 - A[2] + C4; // 2*(A4, B4, C4) - (ax, ay, az)
878 
879  //
880  // Last point
881  //
882  r[0] = R[0] + S * A4;
883  r[1] = R[1] + S * B4;
884  r[2] = R[2] + S * C4; //setup.Field(r,H2);
885  pos.SetXYZ(r[0], r[1], r[2]);
886  TVector3 H2vect = GFFieldManager::getFieldVal(pos);
887  H2[0] = H2vect.X() * PS2;
888  H2[1] = H2vect.Y() * PS2;
889  H2[2] = H2vect.Z() * PS2; // H2 is PS2*(Hx, Hy, Hz) @ (x, y, z) + 0.25*S * (A4, B4, C4)
890  double A6 = B5 * H2[2] - C5 * H2[1], B6 = C5 * H2[0] - A5 * H2[2],
891  C6 = A5 * H2[1] - B5 * H2[0]; // (A5, B5, C5) x H2
892 
893  //
894  // Test approximation quality on given step and possible step reduction
895  //
896  double EST =
897  fabs((A1 + A6) - (A3 + A4)) + fabs((B1 + B6) - (B3 + B4)) +
898  fabs(
899  (C1 + C6) -
900  (C3 +
901  C4)); // EST = ||(ABC1+ABC6)-(ABC3+ABC4)||_1 = ||(axzy x H0 + ABC5 x H2) - (ABC2 x H1 + ABC3 x H1)||_1
902  if (EST > DLT) {
903  S *= 0.5;
904  stopBecauseOfMaterial = false;
905  continue;
906  }
907 
908  //
909  // Derivatives of track parameters in last point
910  //
911  if (calcCov) {
912  for (int i = 7; i != ND;
913  i += 7) { // i = 7, 14, 21, 28, 35, 42, 49; ND = 56; ND1 = 49; rows of Jacobian
914 
915  double* dR = &P[i]; // dR = (dX/dpN, dY/dpN, dZ/dpN)
916  double* dA = &P[i + 3]; // dA = (dAx/dpN, dAy/dpN, dAz/dpN); N = X,Y,Z,Ax,Ay,Az,q/p
917 
918  //first point
919  double dA0 = H0[2] * dA[1] - H0[1] * dA[2]; // dA0/dp }
920  double dB0 = H0[0] * dA[2] - H0[2] * dA[0]; // dB0/dp } = dA x H0
921  double dC0 = H0[1] * dA[0] - H0[0] * dA[1]; // dC0/dp }
922 
923  if (i == ND1) {
924  dA0 += A0;
925  dB0 += B0;
926  dC0 += C0;
927  } // if last row: (dA0, dB0, dC0) := (dA0, dB0, dC0) + (A0, B0, C0)
928 
929  double dA2 = dA0 + dA[0]; // }
930  double dB2 = dB0 + dA[1]; // } = (dA0, dB0, dC0) + dA
931  double dC2 = dC0 + dA[2]; // }
932 
933  //second point
934  double dA3 = dA[0] + dB2 * H1[2] - dC2 * H1[1]; // dA3/dp }
935  double dB3 =
936  dA[1] + dC2 * H1[0] - dA2 * H1[2]; // dB3/dp } = dA + (dA2, dB2, dC2) x H1
937  double dC3 = dA[2] + dA2 * H1[1] - dB2 * H1[0]; // dC3/dp }
938 
939  if (i == ND1) {
940  dA3 += A3 - A[0];
941  dB3 += B3 - A[1];
942  dC3 += C3 - A[2];
943  } // if last row: (dA3, dB3, dC3) := (dA3, dB3, dC3) + (A3, B3, C3) - (ax, ay, az)
944 
945  double dA4 = dA[0] + dB3 * H1[2] - dC3 * H1[1]; // dA4/dp }
946  double dB4 =
947  dA[1] + dC3 * H1[0] - dA3 * H1[2]; // dB4/dp } = dA + (dA3, dB3, dC3) x H1
948  double dC4 = dA[2] + dA3 * H1[1] - dB3 * H1[0]; // dC4/dp }
949 
950  if (i == ND1) {
951  dA4 += A4 - A[0];
952  dB4 += B4 - A[1];
953  dC4 += C4 - A[2];
954  } // if last row: (dA4, dB4, dC4) := (dA4, dB4, dC4) + (A4, B4, C4) - (ax, ay, az)
955 
956  //last point
957  double dA5 = dA4 + dA4 - dA[0]; // }
958  double dB5 = dB4 + dB4 - dA[1]; // } = 2*(dA4, dB4, dC4) - dA
959  double dC5 = dC4 + dC4 - dA[2]; // }
960 
961  double dA6 = dB5 * H2[2] - dC5 * H2[1]; // dA6/dp }
962  double dB6 = dC5 * H2[0] - dA5 * H2[2]; // dB6/dp } = (dA5, dB5, dC5) x H2
963  double dC6 = dA5 * H2[1] - dB5 * H2[0]; // dC6/dp }
964 
965  if (i == ND1) {
966  dA6 += A6;
967  dB6 += B6;
968  dC6 += C6;
969  } // if last row: (dA6, dB6, dC6) := (dA6, dB6, dC6) + (A6, B6, C6)
970 
971  dR[0] += (dA2 + dA3 + dA4) * S3;
972  dA[0] = (dA0 + dA3 + dA3 + dA5 + dA6) *
973  P3; // dR := dR + S3*[(dA2, dB2, dC2) + (dA3, dB3, dC3) + (dA4, dB4, dC4)]
974  dR[1] += (dB2 + dB3 + dB4) * S3;
975  dA[1] =
976  (dB0 + dB3 + dB3 + dB5 + dB6) *
977  P3; // dA := 1/3*[(dA0, dB0, dC0) + 2*(dA3, dB3, dC3) + (dA5, dB5, dC5) + (dA6, dB6, dC6)]
978  dR[2] += (dC2 + dC3 + dC4) * S3;
979  dA[2] = (dC0 + dC3 + dC3 + dC5 + dC6) * P3;
980  }
981  }
982 
983  Way2 += S; // add stepsize to way (signed)
984  if ((Way += fabs(S)) > Wmax) {
985  std::cerr << "PaAlgo::RKutta ==> Trajectory is longer than length limit : " << Way << " cm !"
986  << " p/q = " << 1. / P[6] << " GeV" << std::endl;
987  return (false);
988  }
989 
990  //
991  // Track parameters in last point
992  //
993  R[0] += (A2 + A3 + A4) * S3;
994  A[0] += (SA[0] = (A0 + A3 + A3 + A5 + A6) * P3 -
995  A[0]); // R = R0 + S3*[(A2, B2, C2) + (A3, B3, C3) + (A4, B4, C4)]
996  R[1] += (B2 + B3 + B4) * S3;
997  A[1] +=
998  (SA[1] = (B0 + B3 + B3 + B5 + B6) * P3 -
999  A[1]); // A = 1/3*[(A0, B0, C0) + 2*(A3, B3, C3) + (A5, B5, C5) + (A6, B6, C6)]
1000  R[2] += (C2 + C3 + C4) * S3;
1001  A[2] += (SA[2] = (C0 + C3 + C3 + C5 + C6) * P3 - A[2]); // SA = A_new - A_old
1002  Sl = S; // last S used
1003 
1004  // if extrapolation has changed direction, delete the last point, because it is
1005  // not a consecutive point to be used for material estimations
1006  if (Ssign * S < 0.) {
1007  pointPaths.at(pointPaths.size() - 1) += S;
1008  points.erase(points.end());
1009  }
1010  else {
1011  pointPaths.push_back(S);
1012  }
1013 
1014  points.push_back(TVector3(R[0], R[1], R[2]));
1015 
1016  double CBA = 1. / sqrt(A[0] * A[0] + A[1] * A[1] + A[2] * A[2]); // 1/|A|
1017  A[0] *= CBA;
1018  A[1] *= CBA;
1019  A[2] *= CBA; // normalize A
1020 
1021  // Step estimation until surface and test conditions for stop of propogation
1022  if (fabs(Way2) > Wmax) {
1023  Dis = 0.;
1024  Dist = 0.;
1025  S = 0;
1026  Step = 0.;
1027  break;
1028  }
1029 
1030  An = A[0] * SU[0] + A[1] * SU[1] + A[2] * SU[2];
1031 
1032  if (plane.inActive(TVector3(R[0], R[1], R[2]), TVector3(A[0], A[1], A[2]))) {
1033  Dis = SU[3] - R[0] * SU[0] - R[1] * SU[1] - R[2] * SU[2];
1034  Step = Dis / An;
1035  }
1036  else {
1037  if ((O.X() - R[0]) * A[0] + (O.Y() - R[1]) * A[1] + (O.Z() - R[2]) * A[2] > 0) {
1038  Dis = sqrt((R[0] - O.X()) * (R[0] - O.X()) + (R[1] - O.Y()) * (R[1] - O.Y()) +
1039  (R[2] - O.Z()) * (R[2] - O.Z()));
1040  }
1041  else {
1042  Dis = -1. * sqrt((R[0] - O.X()) * (R[0] - O.X()) + (R[1] - O.Y()) * (R[1] - O.Y()) +
1043  (R[2] - O.Z()) * (R[2] - O.Z()));
1044  }
1045  Step = Dis; // signed distance to surface
1046  }
1047 
1048  // if (An==0 || (Dis*Dist>0 && fabs(Dis)>fabs(Dist))) { // did not get closer to surface
1049  if (Dis * Dist > 0 && fabs(Dis) > fabs(Dist)) { // did not get closer to surface
1050  error = true;
1051  Step = 0;
1052  break;
1053  }
1054  Dist = Dis;
1055 
1056  //
1057  // reset & check step size
1058  //
1059  // reset S to Step if extrapolation too long or in wrong direction
1060  if (S * Step < 0. || fabs(S) > fabs(Step))
1061  S = Step;
1062  else if (EST < DLT32 && fabs(2. * S) <= Smax)
1063  S *= 2.;
1064 
1065  } //end of main loop
1066 
1067  //
1068  // Output information preparation for main track parameteres
1069  //
1070 
1071  if (!stopBecauseOfMaterial) { // linear extrapolation to surface
1072  if (Sl != 0) Sl = 1. / Sl; // Sl = inverted last Stepsize Sl
1073  A[0] += (SA[0] *= Sl) * Step; // Step = distance to surface
1074  A[1] +=
1075  (SA[1] *= Sl) *
1076  Step; // SA*Sl = delta A / delta way; local derivative of A with respect to the length of the way
1077  A[2] += (SA[2] *= Sl) * Step; // A = A + Step * SA*Sl
1078 
1079  P[0] =
1080  R[0] +
1081  Step * (A[0] -
1082  .5 * Step *
1083  SA[0]); // P = R + Step*(A - 1/2*Step*SA); approximation for final point on surface
1084  P[1] = R[1] + Step * (A[1] - .5 * Step * SA[1]);
1085  P[2] = R[2] + Step * (A[2] - .5 * Step * SA[2]);
1086 
1087  points.push_back(TVector3(P[0], P[1], P[2]));
1088  pointPaths.push_back(Step);
1089  }
1090 
1091  double CBA = 1. / sqrt(A[0] * A[0] + A[1] * A[1] + A[2] * A[2]);
1092 
1093  P[3] = A[0] * CBA; // normalize A
1094  P[4] = A[1] * CBA;
1095  P[5] = A[2] * CBA;
1096 
1097  //
1098  // Output derivatives of track parameters preparation
1099  //
1100  An = A[0] * SU[0] + A[1] * SU[1] + A[2] * SU[2];
1101  // An!=0 ? An=1./An : An = 0; // 1/A_normal
1102  fabs(An) < 1.E-6 ? An = 1. / An : An = 0; // 1/A_normal
1103 
1104  if (calcCov && !stopBecauseOfMaterial) {
1105  for (int i = 7; i != ND; i += 7) {
1106  double* dR = &P[i];
1107  double* dA = &P[i + 3];
1108  S = (dR[0] * SU[0] + dR[1] * SU[1] + dR[2] * SU[2]) * An; // dR_normal / A_normal
1109  dR[0] -= S * A[0];
1110  dR[1] -= S * A[1];
1111  dR[2] -= S * A[2];
1112  dA[0] -= S * SA[0];
1113  dA[1] -= S * SA[1];
1114  dA[2] -= S * SA[2];
1115  }
1116  }
1117 
1118  if (error) {
1119  // std::cerr << "PaAlgo::RKutta ==> Do not get closer. Path = " << Way << " cm" << " p/q = " << 1./P[6] << " GeV" << std::endl;
1120  std::cerr << "RKTrackRep::RKutta ==> Do not get closer. Path = " << Way << " cm"
1121  << " p/q = " << 1. / P[6] << " GeV" << std::endl;
1122  return (false);
1123  }
1124 
1125  // calculate covered distance
1126  if (!stopBecauseOfMaterial)
1127  coveredDistance = Way2 + Step;
1128  else
1129  coveredDistance = Way2;
1130 
1131  return (true);
1132 }
TRandom r
Definition: spectrum.C:23
int fPdg
PDG particle code.
Definition: RKTrackRep.h:170
static GFMaterialEffects * getInstance()
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
#define MINSTEP
Definition: RKTrackRep.cxx:39
static TVector3 getFieldVal(const TVector3 &x)
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:47
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:75
double fCharge
Charge.
Definition: RKTrackRep.h:174
Definition: Step.hh:40
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.
void genf::GFAbsTrackRep::setChiSqu ( double  aChiSqu)
inlineinherited

Definition at line 262 of file GFAbsTrackRep.h.

Referenced by genf::GFKalman::fittingPass().

262 { fChiSqu = aChiSqu; }
double fChiSqu
chiSqu of the track fit
Definition: GFAbsTrackRep.h:95
void genf::GFAbsTrackRep::setCov ( const TMatrixT< Double_t > &  aCov)
inlineinherited

Definition at line 244 of file GFAbsTrackRep.h.

Referenced by genf::GFDaf::blowUpCovs(), genf::GFKalman::blowUpCovs(), and genf::GFKalman::blowUpCovsDiag().

244 { fCov = aCov; }
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:92
void genf::RKTrackRep::setData ( const TMatrixT< double > &  st,
const GFDetPlane pl,
const TMatrixT< double > *  cov = NULL,
const TMatrixT< double > *  aux = NULL 
)

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.

Definition at line 41 of file RKTrackRep.cxx.

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

Referenced by switchDirection().

45 {
46  if (aux != NULL) { fCacheSpu = (*aux)(0, 0); }
47  else {
48  if (pl != fCachePlane) {
49  throw GFException("RKTrackRep::setData() called with a reference plane not the same as the "
50  "one the last extrapolate(plane,state,cov) was made",
51  __LINE__,
52  __FILE__)
53  .setFatal();
54  }
55  }
56  GFAbsTrackRep::setData(st, pl, cov);
57  if (fCharge * fState[0][0] < 0) fCharge *= -1; // set charge accordingly! (fState[0][0] = q/p)
58  fSpu = fCacheSpu;
59 }
GFDetPlane fCachePlane
Definition: RKTrackRep.h:160
virtual void setData(const TMatrixT< Double_t > &st, const GFDetPlane &pl, const TMatrixT< Double_t > *cov=NULL)
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:47
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:75
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
double fCharge
Charge.
Definition: RKTrackRep.h:174
virtual void genf::GFAbsTrackRep::setData ( const TMatrixT< Double_t > &  st,
const GFDetPlane pl,
const TMatrixT< Double_t > *  cov = NULL 
)
inlinevirtualinherited

Definition at line 236 of file GFAbsTrackRep.h.

Referenced by genf::GFAbsTrackRep::extrapolate(), genf::GFKalman::processHit(), genf::GFDaf::processTrack(), setData(), and switchDirection().

239  {
240  fState = st;
241  fRefPlane = pl;
242  if (cov != NULL) fCov = *cov;
243  }
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:92
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:89
void genf::GFAbsTrackRep::setFirstCov ( const TMatrixT< Double_t > &  aCov)
inlineinherited

Definition at line 246 of file GFAbsTrackRep.h.

Referenced by genf::GFKalman::processTrack().

246 { fFirstCov = aCov; }
TMatrixT< Double_t > fFirstCov
void genf::GFAbsTrackRep::setFirstPlane ( const GFDetPlane aPlane)
inlineinherited

Definition at line 247 of file GFAbsTrackRep.h.

Referenced by genf::GFKalman::processTrack().

248  {
249  fFirstPlane = aPlane;
250  ;
251  }
GFDetPlane fFirstPlane
void genf::GFAbsTrackRep::setFirstState ( const TMatrixT< Double_t > &  aState)
inlineinherited

Definition at line 245 of file GFAbsTrackRep.h.

Referenced by genf::GFKalman::processTrack().

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

Deprecated. Should be removed soon.

Definition at line 271 of file GFAbsTrackRep.h.

References f.

272  {
273  fInverted = f;
274  return true;
275  }
bool fInverted
specifies the direction of flight of the particle
TFile f
Definition: plotHisto.C:6
void genf::GFAbsTrackRep::setLastCov ( const TMatrixT< Double_t > &  aCov)
inlineinherited

Definition at line 253 of file GFAbsTrackRep.h.

Referenced by genf::GFKalman::processTrack().

253 { fLastCov = aCov; }
TMatrixT< Double_t > fLastCov
void genf::GFAbsTrackRep::setLastPlane ( const GFDetPlane aPlane)
inlineinherited

Definition at line 254 of file GFAbsTrackRep.h.

Referenced by genf::GFKalman::processTrack().

255  {
256  fLastPlane = aPlane;
257  ;
258  }
void genf::GFAbsTrackRep::setLastState ( const TMatrixT< Double_t > &  aState)
inlineinherited

Definition at line 252 of file GFAbsTrackRep.h.

Referenced by genf::GFKalman::processTrack().

252 { fLastState = aState; }
TMatrixT< Double_t > fLastState
void genf::GFAbsTrackRep::setNDF ( unsigned int  n)
inlineinherited

Definition at line 263 of file GFAbsTrackRep.h.

References n.

Referenced by genf::GFKalman::fittingPass().

263 { fNdf = n; }
unsigned int fNdf
Definition: GFAbsTrackRep.h:96
Char_t n[5]
void genf::RKTrackRep::setPDG ( int  i)

Set PDG particle code.

Definition at line 284 of file RKTrackRep.cxx.

References fCharge, fMass, fPdg, and part.

Referenced by RKTrackRep(), and switchDirection().

285 {
286  fPdg = i;
287  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(fPdg);
288  if (part == 0) {
289  std::cerr << "RKTrackRep::setPDG particle " << i << " not known to TDatabasePDG -> abort"
290  << std::endl;
291  exit(1);
292  }
293  fMass = part->Mass();
294  fCharge = part->Charge() / (3.);
295 }
int fPdg
PDG particle code.
Definition: RKTrackRep.h:170
double fMass
Mass (in GeV)
Definition: RKTrackRep.h:172
TString part[npart]
Definition: Style.C:32
double fCharge
Charge.
Definition: RKTrackRep.h:174
void genf::GFAbsTrackRep::setStatusFlag ( int  _val)
inlineinherited

Definition at line 266 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::switchDirection().

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

266 { fStatusFlag = _val; }
int fStatusFlag
status of track representation: 0 means everything&#39;s OK
Definition: GFAbsTrackRep.h:99
void genf::GFAbsTrackRep::stepalong ( double  h)
virtualinherited

make step of h cm along the track

There is an emply implementation in GFAbsTrackRep.cxx which will abort (see one of the extrapolate methods above). This can be overwritten, if this feature is needed.

Definition at line 96 of file GFAbsTrackRep.cxx.

References genf::GFAbsTrackRep::Abort().

97 {
98  Abort("stepalong()");
99 }
void Abort(std::string method)
void genf::RKTrackRep::switchDirection ( )
inlinevirtual

Member Data Documentation

TMatrixT<double> genf::RKTrackRep::fAuxInfo
private

Definition at line 163 of file RKTrackRep.h.

Referenced by getAuxInfo().

GFDetPlane genf::RKTrackRep::fCachePlane
private

Definition at line 160 of file RKTrackRep.h.

Referenced by extrapolate(), getAuxInfo(), and setData().

double genf::RKTrackRep::fCacheSpu
private

Definition at line 161 of file RKTrackRep.h.

Referenced by extrapolate(), getAuxInfo(), getMom(), getPosMom(), and setData().

double genf::RKTrackRep::fCharge
private

Charge.

Definition at line 174 of file RKTrackRep.h.

Referenced by Extrap(), getCharge(), RKTrackRep(), RKutta(), setData(), and setPDG().

double genf::GFAbsTrackRep::fChiSqu
protectedinherited

chiSqu of the track fit

Definition at line 95 of file GFAbsTrackRep.h.

Referenced by genf::GFAbsTrackRep::getChiSqu(), and genf::GFAbsTrackRep::Print().

unsigned int genf::GFAbsTrackRep::fDimension
protectedinherited

Dimensionality of track representation.

Definition at line 86 of file GFAbsTrackRep.h.

Referenced by genf::SlTrackRep::extrapolate(), genf::GFAbsTrackRep::extrapolate(), and genf::GFAbsTrackRep::getDim().

bool genf::RKTrackRep::fDirection
private

Definition at line 167 of file RKTrackRep.h.

Referenced by switchDirection().

TMatrixT<Double_t> genf::GFAbsTrackRep::fFirstCov
protectedinherited
GFDetPlane genf::GFAbsTrackRep::fFirstPlane
protectedinherited

Definition at line 109 of file GFAbsTrackRep.h.

Referenced by genf::GFAbsTrackRep::getFirstPlane().

TMatrixT<Double_t> genf::GFAbsTrackRep::fFirstState
protectedinherited

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

Definition at line 104 of file GFAbsTrackRep.h.

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

bool genf::GFAbsTrackRep::fInverted
protectedinherited

specifies the direction of flight of the particle

Definition at line 101 of file GFAbsTrackRep.h.

TMatrixT<Double_t> genf::GFAbsTrackRep::fLastCov
protectedinherited

Definition at line 108 of file GFAbsTrackRep.h.

Referenced by genf::GFAbsTrackRep::getLastCov(), and genf::GFAbsTrackRep::reset().

GFDetPlane genf::GFAbsTrackRep::fLastPlane
protectedinherited

Definition at line 110 of file GFAbsTrackRep.h.

Referenced by genf::GFAbsTrackRep::getLastPlane().

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

Mass (in GeV)

Definition at line 172 of file RKTrackRep.h.

Referenced by setPDG().

unsigned int genf::GFAbsTrackRep::fNdf
protectedinherited

Definition at line 96 of file GFAbsTrackRep.h.

int genf::RKTrackRep::fPdg
private

PDG particle code.

Definition at line 170 of file RKTrackRep.h.

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

double genf::RKTrackRep::fSpu
private
int genf::GFAbsTrackRep::fStatusFlag
protectedinherited

status of track representation: 0 means everything's OK

Definition at line 99 of file GFAbsTrackRep.h.

Referenced by genf::GFAbsTrackRep::getStatusFlag().


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