LArSoft  v07_13_02
Liquid Argon Software toolkit - http://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 52 of file RKTrackRep.h.

Constructor & Destructor Documentation

genf::RKTrackRep::RKTrackRep ( )

Definition at line 72 of file RKTrackRep.cxx.

Referenced by clone(), and prototype().

72  : GFAbsTrackRep(5),
73  fCachePlane(),
74  fCacheSpu(1),
75  fSpu(1),
76  fAuxInfo(1,1),
77  fDirection(true),
78  fPdg(0),
79  fMass(0.),
80  fCharge(-1)
81 { }
int fPdg
PDG particle code.
Definition: RKTrackRep.h:177
double fMass
Mass (in GeV)
Definition: RKTrackRep.h:179
GFDetPlane fCachePlane
Definition: RKTrackRep.h:166
TMatrixT< double > fAuxInfo
Definition: RKTrackRep.h:169
double fCharge
Charge.
Definition: RKTrackRep.h:181
genf::RKTrackRep::RKTrackRep ( const TVector3 &  pos,
const TVector3 &  mom,
const TVector3 &  poserr,
const TVector3 &  momerr,
const int &  PDGCode 
)

Definition at line 84 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.

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

Definition at line 141 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.

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

Definition at line 192 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.

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

Definition at line 246 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.

248  :
249  GFAbsTrackRep(5),
250  fCachePlane(), fCacheSpu(1), fAuxInfo(1,1),
251  fDirection(true)
252 {
253  setPDG(PDGCode); // also sets charge and mass
254 
255  //fEffect = new GFMaterialEffects();
256 
257  fRefPlane = pl;
258  TVector3 o=fRefPlane.getO();
259  TVector3 u=fRefPlane.getU();
260  TVector3 v=fRefPlane.getV();
261  TVector3 w=u.Cross(v);
262 
263  double pu=mom*u;
264  double pv=mom*v;
265  double pw=mom*w;
266 
267  fState[0][0]=fCharge/mom.Mag();
268  //u' and v'
269  fState[1][0]=pu/pw;
270  fState[2][0]=pv/pw;
271  //u and v
272  fState[3][0]=0.0;
273  fState[4][0]=0.0;
274  //spu
275  fSpu=1.;
276 
277  static const TVector3 stdPosErr2(1.,1.,1.);
278  static const TVector3 stdMomErr2(1.,1.,1.);
279 
280  fCov[3][3] = stdPosErr2.X()*stdPosErr2.X() * u.X()*u.X() +
281  stdPosErr2.Y()*stdPosErr2.Y() * u.Y()*u.Y() +
282  stdPosErr2.Z()*stdPosErr2.Z() * u.Z()*u.Z();
283  fCov[4][4] = stdPosErr2.X()*stdPosErr2.X() * v.X()*v.X() +
284  stdPosErr2.Y()*stdPosErr2.Y() * v.Y()*v.Y() +
285  stdPosErr2.Z()*stdPosErr2.Z() * v.Z()*v.Z();
286  fCov[0][0] = fCharge*fCharge/pow(mom.Mag(),6.) *
287  (mom.X()*mom.X() * stdMomErr2.X()*stdMomErr2.X()+
288  mom.Y()*mom.Y() * stdMomErr2.Y()*stdMomErr2.Y()+
289  mom.Z()*mom.Z() * stdMomErr2.Z()*stdMomErr2.Z());
290  fCov[1][1] = pow((u.X()/pw - w.X()*pu/(pw*pw)),2.)*stdMomErr2.X()*stdMomErr2.X() +
291  pow((u.Y()/pw - w.Y()*pu/(pw*pw)),2.)*stdMomErr2.Y()*stdMomErr2.Y() +
292  pow((u.Z()/pw - w.Z()*pu/(pw*pw)),2.)*stdMomErr2.Z()*stdMomErr2.Z();
293  fCov[2][2] = pow((v.X()/pw - w.X()*pv/(pw*pw)),2.)*stdMomErr2.X()*stdMomErr2.X() +
294  pow((v.Y()/pw - w.Y()*pv/(pw*pw)),2.)*stdMomErr2.Y()*stdMomErr2.Y() +
295  pow((v.Z()/pw - w.Z()*pv/(pw*pw)),2.)*stdMomErr2.Z()*stdMomErr2.Z();
296 
297 }
GFDetPlane fCachePlane
Definition: RKTrackRep.h:166
TMatrixT< double > fAuxInfo
Definition: RKTrackRep.h:169
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:96
TVector3 getO() const
Definition: GFDetPlane.h:76
void setPDG(int)
Set PDG particle code.
Definition: RKTrackRep.cxx:300
TVector3 getU() const
Definition: GFDetPlane.h:77
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:93
TVector3 getV() const
Definition: GFDetPlane.h:78
double fCharge
Charge.
Definition: RKTrackRep.h:181
Float_t w
Definition: plot.C:23
genf::RKTrackRep::~RKTrackRep ( )
virtual

Definition at line 66 of file RKTrackRep.cxx.

66  {
67  // delete fEffect;
68 //GFMaterialEffects is now a Singleton
69 }
genf::RKTrackRep::RKTrackRep ( const RKTrackRep )
inlineprivate

Definition at line 173 of file RKTrackRep.h.

Member Function Documentation

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

Definition at line 293 of file GFAbsTrackRep.h.

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

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

Definition at line 296 of file GFAbsTrackRep.h.

References n.

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

296  {
297  fNdf += n;
298  }
Char_t n[5]
virtual GFAbsTrackRep* genf::RKTrackRep::clone ( ) const
inlinevirtual

Implements genf::GFAbsTrackRep.

Definition at line 77 of file RKTrackRep.h.

References RKTrackRep().

77 {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 1059 of file RKTrackRep.cxx.

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

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

1059  {
1060 
1061  static const int maxNumIt(2000);
1062  int numIt(0);
1063  bool calcCov(true);
1064  if(cov==NULL) calcCov=false;
1065 
1066  double P[56]; // only 7 used if no covariance matrix
1067  if(calcCov) std::fill(P, P + std::extent<decltype(P)>::value, 0);
1068 // else {} // not needed std::fill(P, P + 7, 0);
1069 
1070  for(int i=0;i<7;++i){
1071  P[i] = (*state)[i][0];
1072  }
1073 
1074  TMatrixT<Double_t> jac(7,7);
1075  TMatrixT<Double_t> jacT(7,7);
1076  TMatrixT<Double_t> oldCov(7,7);
1077  if(calcCov) oldCov=(*cov);
1078  double coveredDistance(0.);
1079  double sumDistance(0.);
1080 
1081  while(true){
1082  if(numIt++ > maxNumIt){
1083  throw GFException("RKTrackRep::Extrap ==> maximum number of iterations exceeded",
1084  __LINE__,__FILE__).setFatal();
1085  }
1086 
1087  if(calcCov){
1088  memset(&P[7],0x00,49*sizeof(double));
1089  for(int i=0; i<6; ++i){
1090  P[(i+1)*7+i] = 1.;
1091  }
1092  P[55] = (*state)[6][0];
1093  }
1094 
1095 // double dir(1.);
1096  {
1097  TVector3 Pvect(P[0],P[1],P[2]); //position
1098  TVector3 Avect(P[3],P[4],P[5]); //direction
1099  TVector3 dist = plane.dist(Pvect); //from point to plane
1100  // std::cout << "RKTrackRep::Extrap(): Position, Direction, Plane, dist " << std::endl;
1101  //Pvect.Print();
1102  //Avect.Print();
1103  //plane.Print();
1104  //dist.Print();
1105  // if(dist*Avect<0.) dir=-1.;
1106  }
1107 
1108  TVector3 directionBefore(P[3],P[4],P[5]); // direction before propagation
1109  directionBefore.SetMag(1.);
1110 
1111  // propagation
1112  std::vector<TVector3> points;
1113  std::vector<double> pointPaths;
1114  if( ! this->RKutta(plane,P,coveredDistance,points,pointPaths,-1.,calcCov) ) { // maxLen currently not used
1115  //GFException exc("RKTrackRep::Extrap ==> Runge Kutta propagation failed",__LINE__,__FILE__);
1116 
1117 
1118  //exc.setFatal(); // stops propagation; faster, but some hits will be lost
1119  mf::LogInfo("RKTrackRep::RKutta(): ") << "Throw cet exception here, ... ";
1120  throw GFException("RKTrackRep::RKutta(): Runge Kutta propagation failed",
1121  __LINE__,__FILE__).setFatal();
1122  }
1123 
1124  TVector3 directionAfter(P[3],P[4],P[5]); // direction after propagation
1125  directionAfter.SetMag(1.);
1126 
1127  sumDistance+=coveredDistance;
1128 
1129  // filter Points
1130  std::vector<TVector3> pointsFilt(1, points.at(0));
1131  std::vector<double> pointPathsFilt(1, 0.);
1132  // only if in right direction
1133  for(unsigned int i=1;i<points.size();++i){
1134  if (pointPaths.at(i) * coveredDistance > 0.) {
1135  pointsFilt.push_back(points.at(i));
1136  pointPathsFilt.push_back(pointPaths.at(i));
1137  }
1138  else {
1139  pointsFilt.back() = points.at(i);
1140  pointPathsFilt.back() += pointPaths.at(i);
1141  }
1142  // clean up tiny steps
1143  int position = pointsFilt.size()-1; // position starts with 0
1144  if (fabs(pointPathsFilt.back()) < MINSTEP && position > 1) {
1145  pointsFilt.at(position-1) = pointsFilt.at(position);
1146  pointsFilt.pop_back();
1147  pointPathsFilt.at(position-1) += pointPathsFilt.at(position);
1148  pointPathsFilt.pop_back();
1149  }
1150  }
1151 
1152  //consistency check
1153  double checkSum(0.);
1154  for(unsigned int i=0;i<pointPathsFilt.size();++i){
1155  checkSum+=pointPathsFilt.at(i);
1156  }
1157  //assert(fabs(checkSum-coveredDistance)<1.E-7);
1158  if(fabs(checkSum-coveredDistance)>1.E-7){
1159  throw GFException("RKTrackRep::Extrap ==> fabs(checkSum-coveredDistance)>1.E-7",__LINE__,__FILE__).setFatal();
1160  }
1161 
1162  if(calcCov){ //calculate Jacobian jac
1163  for(int i=0;i<7;++i){
1164  for(int j=0;j<7;++j){
1165  if(i<6) jac[i][j] = P[ (i+1)*7+j ];
1166  else jac[i][j] = P[ (i+1)*7+j ]/P[6];
1167  }
1168  }
1169  jacT = jac;
1170  jacT.T();
1171  }
1172 
1173  TMatrixT<Double_t> noise(7,7); // zero everywhere by default
1174 
1175  // call MatEffects
1176  double momLoss; // momLoss has a sign - negative loss means momentum gain
1177 
1178  /*
1179  momLoss = fEffect->effects(pointsFilt,
1180  pointPathsFilt,
1181  fabs(fCharge/P[6]), // momentum
1182  fPdg,
1183  calcCov,
1184  &noise,
1185  &jac,
1186  &directionBefore,
1187  &directionAfter);
1188  */
1189  momLoss = GFMaterialEffects::getInstance()->effects(pointsFilt,
1190  pointPathsFilt,
1191  fabs(fCharge/P[6]), // momentum
1192  fPdg,
1193  calcCov,
1194  &noise,
1195  &jac,
1196  &directionBefore,
1197  &directionAfter);
1198 
1199  if(fabs(P[6])>1.E-10){ // do momLoss only for defined 1/momentum .ne.0
1200  P[6] = fCharge/(fabs(fCharge/P[6])-momLoss);
1201  }
1202 
1203  if(calcCov){ //propagate cov and add noise
1204  oldCov = *cov;
1205  *cov = jacT*((oldCov)*jac)+noise;
1206  }
1207 
1208 
1209  //we arrived at the destination plane, if we point to the active area
1210  //of the plane (if it is finite), and the distance is below threshold
1211  if( plane.inActive(TVector3(P[0],P[1],P[2]),TVector3(P[3],P[4],P[5]))) {
1212  if(plane.distance(P[0],P[1],P[2])<MINSTEP) break;
1213  }
1214  }
1215  (*state)[0][0] = P[0]; (*state)[1][0] = P[1];
1216  (*state)[2][0] = P[2]; (*state)[3][0] = P[3];
1217  (*state)[4][0] = P[4]; (*state)[5][0] = P[5];
1218  (*state)[6][0] = P[6];
1219 
1220  return sumDistance;
1221 }
int fPdg
PDG particle code.
Definition: RKTrackRep.h:177
static GFMaterialEffects * getInstance()
Float_t E
Definition: plot.C:23
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
unsigned int noise()
Definition: chem4.cc:265
void fill(const art::PtrVector< recob::Hit > &hits, int only_plane)
#define MINSTEP
Definition: RKTrackRep.cxx:40
std::string value(boost::any const &)
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:50
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:80
double fCharge
Charge.
Definition: RKTrackRep.h:181
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:692
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 476 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().

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

603  {
604 
605  TVector3 o=fRefPlane.getO();
606  TVector3 u=fRefPlane.getU();
607  TVector3 v=fRefPlane.getV();
608  TVector3 w=u.Cross(v);
609 
610 
611  TVector3 pTilde = fSpu* (w + fState[1][0] * u + fState[2][0] * v);
612  double pTildeMag = pTilde.Mag();
613 
614  TVector3 pos = o + fState[3][0]*u + fState[4][0]*v;
615 
616  TMatrixT<Double_t> state7(7,1);
617  state7[0][0] = pos.X();
618  state7[1][0] = pos.Y();
619  state7[2][0] = pos.Z();
620  state7[3][0] = pTilde.X()/pTildeMag;
621  state7[4][0] = pTilde.Y()/pTildeMag;
622  state7[5][0] = pTilde.Z()/pTildeMag;
623  state7[6][0] = fState[0][0];
624 
625  TVector3 O = pl.getO();
626  TVector3 U = pl.getU();
627  TVector3 V = pl.getV();
628  TVector3 W = pl.getNormal();
629 
630  double coveredDistance = this->Extrap(pl,&state7);
631 
632  double X = state7[0][0];
633  double Y = state7[1][0];
634  double Z = state7[2][0];
635  double AX = state7[3][0];
636  double AY = state7[4][0];
637  double AZ = state7[5][0];
638  double QOP = state7[6][0];
639  TVector3 A(AX,AY,AZ);
640  TVector3 Point(X,Y,Z);
641 
642  statePred.ResizeTo(5,1);
643  statePred[0][0] = QOP;
644  statePred[1][0] = (A*U)/(A*W);
645  statePred[2][0] = (A*V)/(A*W);
646  statePred[3][0] = (Point-O)*U;
647  statePred[4][0] = (Point-O)*V;
648 
649  return coveredDistance;
650 }
Float_t Y
Definition: plot.C:39
Float_t Z
Definition: plot.C:39
TVector3 getO() const
Definition: GFDetPlane.h:76
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:34
TVector3 getU() const
Definition: GFDetPlane.h:77
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:93
TVector3 getV() const
Definition: GFDetPlane.h:78
Float_t X
Definition: plot.C:39
Float_t w
Definition: plot.C:23
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 33 of file GFAbsTrackRep.cxx.

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

33  {
34  TMatrixT<Double_t> statePred(fDimension,1);
35  TMatrixT<Double_t> covPred(fDimension,fDimension);
36  double retVal = extrapolate(plane,statePred,covPred);
37  setData(statePred,plane,&covPred);
38  return retVal;
39 }
unsigned int fDimension
Dimensionality of track representation.
Definition: GFAbsTrackRep.h:90
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 425 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().

429  {
430  static const int maxIt(30);
431 
432  TVector3 o=fRefPlane.getO();
433  TVector3 u=fRefPlane.getU();
434  TVector3 v=fRefPlane.getV();
435  TVector3 w=u.Cross(v);
436 
437  TVector3 pTilde = fSpu* (w + fState[1][0] * u + fState[2][0] * v);
438  pTilde.SetMag(1.);
439 
440  TVector3 point = o + fState[3][0]*u + fState[4][0]*v;
441 
442  TMatrixT<Double_t> state7(7,1);
443  state7[0][0] = point.X();
444  state7[1][0] = point.Y();
445  state7[2][0] = point.Z();
446  state7[3][0] = pTilde.X();
447  state7[4][0] = pTilde.Y();
448  state7[5][0] = pTilde.Z();
449  state7[6][0] = fState[0][0];
450 
451  double coveredDistance(0.);
452 
453  GFDetPlane pl;
454  int iterations(0);
455 
456  while(true){
457  pl.setO(point1);
458  TVector3 currentDir(state7[3][0],state7[4][0],state7[5][0]);
459  pl.setU(currentDir.Cross(point2-point1));
460  pl.setV(point2-point1);
461  coveredDistance = this->Extrap(pl,&state7);
462 
463  if(fabs(coveredDistance)<MINSTEP) break;
464  if(++iterations == maxIt) {
465  throw GFException("RKTrackRep extrapolation to point failed, maximum number of iterations reached",__LINE__,__FILE__).setFatal();
466  }
467  }
468  poca.SetXYZ(state7[0][0],state7[1][0],state7[2][0]);
469  dirInPoca.SetXYZ(state7[3][0],state7[4][0],state7[5][0]);
470  poca_onwire = poca2Line(point1,point2,poca);
471 }
TVector3 getO() const
Definition: GFDetPlane.h:76
#define MINSTEP
Definition: RKTrackRep.cxx:40
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:50
TVector3 getU() const
Definition: GFDetPlane.h:77
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:80
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:93
TVector3 getV() const
Definition: GFDetPlane.h:78
Float_t w
Definition: plot.C:23
TVector3 poca2Line(const TVector3 &extr1, const TVector3 &extr2, const TVector3 &point) const
Definition: RKTrackRep.cxx:412
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 366 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().

368  {
369 
370  static const int maxIt(30);
371 
372  TVector3 o=fRefPlane.getO();
373  TVector3 u=fRefPlane.getU();
374  TVector3 v=fRefPlane.getV();
375  TVector3 w=u.Cross(v);
376 
377  TVector3 pTilde = fSpu* (w + fState[1][0] * u + fState[2][0] * v);
378  pTilde.SetMag(1.);
379 
380  TVector3 point = o + fState[3][0]*u + fState[4][0]*v;
381 
382  TMatrixT<Double_t> state7(7,1);
383  state7[0][0] = point.X();
384  state7[1][0] = point.Y();
385  state7[2][0] = point.Z();
386  state7[3][0] = pTilde.X();
387  state7[4][0] = pTilde.Y();
388  state7[5][0] = pTilde.Z();
389  state7[6][0] = fState[0][0];
390 
391  double coveredDistance(0.);
392 
393  GFDetPlane pl;
394  int iterations(0);
395 
396  while(true){
397  pl.setON(pos,TVector3(state7[3][0],state7[4][0],state7[5][0]));
398  coveredDistance = this->Extrap(pl,&state7);
399 
400  if(fabs(coveredDistance)<MINSTEP) break;
401  if(++iterations == maxIt) {
402  throw GFException("RKTrackRep::extrapolateToPoint==> extrapolation to point failed, maximum number of iterations reached",__LINE__,__FILE__).setFatal();
403  }
404  }
405  poca.SetXYZ(state7[0][0],state7[1][0],state7[2][0]);
406  dirInPoca.SetXYZ(state7[3][0],state7[4][0],state7[5][0]);
407 }
TVector3 getO() const
Definition: GFDetPlane.h:76
#define MINSTEP
Definition: RKTrackRep.cxx:40
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:50
TVector3 getU() const
Definition: GFDetPlane.h:77
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:80
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:93
TVector3 getV() const
Definition: GFDetPlane.h:78
Float_t w
Definition: plot.C:23
const TMatrixT< double > * genf::RKTrackRep::getAuxInfo ( const GFDetPlane pl)

Definition at line 56 of file RKTrackRep.cxx.

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

Referenced by switchDirection().

56  {
57 
58  if(pl!=fCachePlane) {
59  throw GFException("RKTrackRep::getAuxInfo() trying to get auxillary information with planes mismatch (Information returned does not belong to requested plane)", __LINE__, __FILE__).setFatal();
60  }
61  fAuxInfo.ResizeTo(1,1);
62  fAuxInfo(0,0) = fCacheSpu;
63  return &fAuxInfo;
64 
65 }
GFDetPlane fCachePlane
Definition: RKTrackRep.h:166
TMatrixT< double > fAuxInfo
Definition: RKTrackRep.h:169
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:50
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:80
double genf::RKTrackRep::getCharge ( ) const
inlinevirtual

Returns charge.

Implements genf::GFAbsTrackRep.

Definition at line 141 of file RKTrackRep.h.

References fCharge.

141 {return fCharge;}
double fCharge
Charge.
Definition: RKTrackRep.h:181
double genf::GFAbsTrackRep::getChiSqu ( ) const
inlineinherited

Definition at line 245 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fChiSqu.

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

245  {
246  return fChiSqu;
247  }
double fChiSqu
chiSqu of the track fit
Definition: GFAbsTrackRep.h:99
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 194 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().

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

Definition at line 230 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fFirstCov.

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

Definition at line 233 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fFirstPlane.

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

Definition at line 227 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fFirstState.

227  {
228  return fFirstState;
229  }
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 239 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fLastCov.

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

Definition at line 242 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fLastPlane.

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

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

Definition at line 236 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fLastState.

236  {
237  return fLastState;
238  }
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 324 of file RKTrackRep.cxx.

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

324  {
325  TMatrixT<Double_t> statePred(fState);
326  TVector3 retmom;
327  if(pl!=fRefPlane) {
328  extrapolate(pl,statePred);
329  retmom = fCacheSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
330  }
331  else{
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 }
double extrapolate(const GFDetPlane &, TMatrixT< Double_t > &statePred, TMatrixT< Double_t > &covPred)
returns the tracklength spanned in this extrapolation
Definition: RKTrackRep.cxx:476
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:93
TVector3 genf::GFAbsTrackRep::getMom ( )
inlineinherited
TVector3 genf::RKTrackRep::getMomLast ( const GFDetPlane pl)

Definition at line 338 of file RKTrackRep.cxx.

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

Referenced by prototype().

338  {
339  TMatrixT<Double_t> statePred(fLastState);
340  TVector3 retmom;
341  retmom = fSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
342 
343  retmom.SetMag(1./fabs(statePred[0][0]));
344  return retmom;
345 }
TMatrixT< Double_t > fLastState
unsigned int genf::GFAbsTrackRep::getNDF ( ) const
inlineinherited

Definition at line 253 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::getDim().

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

253  {
254  if(fNdf>getDim()) return fNdf-getDim();
255  return 0;
256  }
unsigned int getDim() const
returns dimension of state vector
int genf::RKTrackRep::getPDG ( )

Definition at line 312 of file RKTrackRep.cxx.

References fPdg.

Referenced by switchDirection().

312 {return fPdg;}
int fPdg
PDG particle code.
Definition: RKTrackRep.h:177
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 314 of file RKTrackRep.cxx.

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

314  {
315  if(pl!=fRefPlane){
316  TMatrixT<Double_t> s(5,1);
317  extrapolate(pl,s);
318  return pl.getO()+s[3][0]*pl.getU()+s[4][0]*pl.getV();
319  }
320  return fRefPlane.getO()+fState[3][0]*fRefPlane.getU()+fState[4][0]*fRefPlane.getV();
321 }
Float_t s
Definition: plot.C:23
double extrapolate(const GFDetPlane &, TMatrixT< Double_t > &statePred, TMatrixT< Double_t > &covPred)
returns the tracklength spanned in this extrapolation
Definition: RKTrackRep.cxx:476
TVector3 getO() const
Definition: GFDetPlane.h:76
TVector3 getU() const
Definition: GFDetPlane.h:77
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:93
TVector3 getV() const
Definition: GFDetPlane.h:78
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 348 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().

349  {
350  TMatrixT<Double_t> statePred(fState);
351  if(pl!=fRefPlane) {
352  extrapolate(pl,statePred);
353  mom = fCacheSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
354  }
355  else{
356  mom = fSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
357  }
358  mom.SetMag(1./fabs(statePred[0][0]));
359  pos = pl.getO()+(statePred[3][0]*pl.getU())+(statePred[4][0]*pl.getV());
360 }
double extrapolate(const GFDetPlane &, TMatrixT< Double_t > &statePred, TMatrixT< Double_t > &covPred)
returns the tracklength spanned in this extrapolation
Definition: RKTrackRep.cxx:476
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:93
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 75 of file GFAbsTrackRep.cxx.

References genf::GFAbsTrackRep::Abort().

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

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

Definition at line 223 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::getPosMomCov().

223  {
224  getPosMomCov(fRefPlane,pos,mom,c);
225  }
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 249 of file GFAbsTrackRep.h.

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

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

249  {
250  if(getNDF()>0) return getChiSqu()/getNDF();
251  return 0;
252  }
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 201 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::fState.

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

Definition at line 162 of file RKTrackRep.h.

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

Definition at line 172 of file RKTrackRep.h.

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

Definition at line 412 of file RKTrackRep.cxx.

References GFException::setFatal().

Referenced by extrapolateToLine().

412  {
413 
414  TVector3 theWire = extr2-extr1;
415  if(theWire.Mag()<1.E-8){
416  throw GFException("RKTrackRep::poca2Line(): try to find poca between line and point, but the line is really just a point",__LINE__,__FILE__).setFatal();
417  }
418  double t = 1./(theWire*theWire)*(point*theWire+extr1*extr1-extr1*extr2);
419  return (extr1+t*theWire);
420 }
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:50
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:80
void genf::GFAbsTrackRep::Print ( std::ostream &  out = std::cout) const
virtualinherited

Definition at line 93 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().

93  {
94  out << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
95  out <<"GFAbsTrackRep::Parameters at reference plane ";
96  fRefPlane.Print(out);
97  out <<"GFAbsTrackRep::State"<<std::endl;
98  PrintROOTmatrix(out, fState);
99  out <<"GFAbsTrackRep::Covariances"<<std::endl;
100  PrintROOTmatrix(out, fCov);
101  out <<"GFAbsTrackRep::chi^2"<<std::endl;
102  out <<fChiSqu<<std::endl;
103  out << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
104 }
double fChiSqu
chiSqu of the track fit
Definition: GFAbsTrackRep.h:99
void PrintROOTmatrix(std::ostream &out, const TMatrixT< T > &m)
Small utility functions which print some ROOT objects into an output stream.
Definition: GFException.h:132
void Print(std::ostream &out=std::cout) const
Definition: GFDetPlane.cxx:242
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:96
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:93
virtual GFAbsTrackRep* genf::RKTrackRep::prototype ( ) const
inlinevirtual
void genf::RKTrackRep::rescaleCovOffDiags ( )

Definition at line 1223 of file RKTrackRep.cxx.

References genf::GFAbsTrackRep::fCov.

Referenced by extrapolate(), and switchDirection().

1224 {
1225 
1226  for(int i=0;i<fCov.GetNrows();++i){
1227  for(int j=0;j<fCov.GetNcols();++j){
1228  if(i!=j){//off diagonal
1229  fCov[i][j]=0.5*fCov[i][j];
1230  }
1231  else{//diagonal. Do not let diag cov element be <=0.
1232  if (fCov[i][j]<=0.0) fCov[i][j] = 0.01;
1233  }
1234  }
1235  }
1236 }
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:96
void genf::GFAbsTrackRep::reset ( )
virtualinherited

Definition at line 80 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().

80  {
81  std::cout<<"GFAbsTrackRep::reset"<<std::endl;
82  TVector3 nullVec(0.,0.,0.);
83  fRefPlane.set(nullVec,nullVec,nullVec);
84  fState.Zero();
85  fCov.Zero();
86  fFirstState.Zero();
87  fFirstCov.Zero();
88  fLastState.Zero();
89  fLastCov.Zero();
90 }
void set(const TVector3 &o, const TVector3 &u, const TVector3 &v)
Definition: GFDetPlane.cxx:85
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:96
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:93
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 692 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().

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

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

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

Definition at line 263 of file GFAbsTrackRep.h.

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

263  {
264  fCov=aCov;
265  }
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:96
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 43 of file RKTrackRep.cxx.

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

Referenced by switchDirection().

43  {
44  if(aux != NULL) {
45  fCacheSpu = (*aux)(0,0);
46  } else {
47  if(pl!=fCachePlane){
48  throw GFException("RKTrackRep::setData() called with a reference plane not the same as the one the last extrapolate(plane,state,cov) was made", __LINE__, __FILE__).setFatal();
49  }
50  }
51  GFAbsTrackRep::setData(st,pl,cov);
52  if (fCharge*fState[0][0] < 0) fCharge *= -1; // set charge accordingly! (fState[0][0] = q/p)
54 }
GFDetPlane fCachePlane
Definition: RKTrackRep.h:166
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:50
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:80
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:93
double fCharge
Charge.
Definition: RKTrackRep.h:181
virtual void genf::GFAbsTrackRep::setData ( const TMatrixT< Double_t > &  st,
const GFDetPlane pl,
const TMatrixT< Double_t > *  cov = NULL 
)
inlinevirtualinherited

Definition at line 258 of file GFAbsTrackRep.h.

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

258  {
259  fState=st;
260  fRefPlane=pl;
261  if(cov!=NULL) fCov=*cov;
262  }
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:96
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:93
void genf::GFAbsTrackRep::setFirstCov ( const TMatrixT< Double_t > &  aCov)
inlineinherited

Definition at line 269 of file GFAbsTrackRep.h.

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

Definition at line 272 of file GFAbsTrackRep.h.

272  {
273  fFirstPlane = aPlane;;
274  }
GFDetPlane fFirstPlane
void genf::GFAbsTrackRep::setFirstState ( const TMatrixT< Double_t > &  aState)
inlineinherited

Definition at line 266 of file GFAbsTrackRep.h.

266  {
267  fFirstState = aState;
268  }
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 306 of file GFAbsTrackRep.h.

References f.

306 {fInverted=f; return true;}
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 278 of file GFAbsTrackRep.h.

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

Definition at line 281 of file GFAbsTrackRep.h.

281  {
282  fLastPlane = aPlane;;
283  }
void genf::GFAbsTrackRep::setLastState ( const TMatrixT< Double_t > &  aState)
inlineinherited

Definition at line 275 of file GFAbsTrackRep.h.

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

Definition at line 290 of file GFAbsTrackRep.h.

References n.

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

290  {
291  fNdf = n;
292  }
Char_t n[5]
void genf::RKTrackRep::setPDG ( int  i)

Set PDG particle code.

Definition at line 300 of file RKTrackRep.cxx.

References fCharge, fMass, fPdg, and part.

Referenced by RKTrackRep(), and switchDirection().

300  {
301  fPdg = i;
302  TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(fPdg);
303  if(part == 0){
304  std::cerr << "RKTrackRep::setPDG particle " << i
305  << " not known to TDatabasePDG -> abort" << std::endl;
306  exit(1);
307  }
308  fMass = part->Mass();
309  fCharge = part->Charge()/(3.);
310 }
int fPdg
PDG particle code.
Definition: RKTrackRep.h:177
double fMass
Mass (in GeV)
Definition: RKTrackRep.h:179
TString part[npart]
Definition: Style.C:32
double fCharge
Charge.
Definition: RKTrackRep.h:181
void genf::GFAbsTrackRep::setStatusFlag ( int  _val)
inlineinherited

Definition at line 299 of file GFAbsTrackRep.h.

References genf::GFAbsTrackRep::switchDirection().

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

299  {
300  fStatusFlag = _val;
301  }
int fStatusFlag
status of track representation: 0 means everything&#39;s OK
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 71 of file GFAbsTrackRep.cxx.

References genf::GFAbsTrackRep::Abort().

71  {
72  Abort("stepalong()");
73 }
void Abort(std::string method)
void genf::RKTrackRep::switchDirection ( )
inlinevirtual

Member Data Documentation

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

Definition at line 169 of file RKTrackRep.h.

Referenced by getAuxInfo().

GFDetPlane genf::RKTrackRep::fCachePlane
private

Definition at line 166 of file RKTrackRep.h.

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

double genf::RKTrackRep::fCacheSpu
private

Definition at line 167 of file RKTrackRep.h.

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

double genf::RKTrackRep::fCharge
private

Charge.

Definition at line 181 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 99 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 90 of file GFAbsTrackRep.h.

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

bool genf::RKTrackRep::fDirection
private

Definition at line 174 of file RKTrackRep.h.

Referenced by switchDirection().

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

Definition at line 113 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 108 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 105 of file GFAbsTrackRep.h.

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

Definition at line 112 of file GFAbsTrackRep.h.

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

GFDetPlane genf::GFAbsTrackRep::fLastPlane
protectedinherited

Definition at line 114 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 179 of file RKTrackRep.h.

Referenced by setPDG().

unsigned int genf::GFAbsTrackRep::fNdf
protectedinherited

Definition at line 100 of file GFAbsTrackRep.h.

int genf::RKTrackRep::fPdg
private

PDG particle code.

Definition at line 177 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 103 of file GFAbsTrackRep.h.

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


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