All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
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 GFAbsTrackRep * clone () const
 
virtual GFAbsTrackRep * prototype () 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 ()
 
- Public Member Functions inherited from genf::GFAbsTrackRep
 GFAbsTrackRep ()
 
 GFAbsTrackRep (int)
 
virtual ~GFAbsTrackRep ()
 
virtual void stepalong (double h)
 make step of h cm along the track More...
 
double extrapolate (const GFDetPlane &plane)
 This changes the state and cov and plane of the rep. More...
 
unsigned int getDim () const
 returns dimension of state vector More...
 
virtual void Print (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
 
virtual void getPosMomCov (const GFDetPlane &pl, TVector3 &pos, TVector3 &mom, TMatrixT< Double_t > &cov)
 method which gets position, momentum and 6x6 covariance matrix More...
 
TVector3 getPos ()
 
TVector3 getMom ()
 
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 GFDetPlane & getReferencePlane () 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 ()
 

Private Member Functions

RKTrackRep & operator= (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...
 

Additional Inherited Members

- Protected Attributes inherited from genf::GFAbsTrackRep
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
 

Detailed Description

Definition at line 53 of file RKTrackRep.h.

Constructor & Destructor Documentation

genf::RKTrackRep::RKTrackRep ( )

Definition at line 71 of file RKTrackRep.cxx.

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

Definition at line 83 of file RKTrackRep.cxx.

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

Definition at line 140 of file RKTrackRep.cxx.

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

Definition at line 191 of file RKTrackRep.cxx.

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

Definition at line 245 of file RKTrackRep.cxx.

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

Definition at line 65 of file RKTrackRep.cxx.

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

Definition at line 174 of file RKTrackRep.h.

174 : GFAbsTrackRep() {}

Member Function Documentation

virtual GFAbsTrackRep* genf::RKTrackRep::clone ( ) const
inlinevirtual

Implements genf::GFAbsTrackRep.

Definition at line 78 of file RKTrackRep.h.

78 {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 1058 of file RKTrackRep.cxx.

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

477  {
478 
479  TMatrixT<Double_t> cov7x7(7,7);
480  TMatrixT<Double_t> J_pM(7,5);
481 
482  TVector3 o=fRefPlane.getO();
483  TVector3 u=fRefPlane.getU();
484  TVector3 v=fRefPlane.getV();
485  TVector3 w=u.Cross(v);
486  std::ostream* pOut = nullptr; // &std::cout if you really really want
487 
488  J_pM[0][3] = u.X();J_pM[0][4]=v.X(); // dx/du
489  J_pM[1][3] = u.Y();J_pM[1][4]=v.Y();
490  J_pM[2][3] = u.Z();J_pM[2][4]=v.Z();
491 
492  TVector3 pTilde = fSpu * (w + fState[1][0] * u + fState[2][0] * v);
493  double pTildeMag = pTilde.Mag();
494 
495  //J_pM matrix is d(x,y,z,ax,ay,az,q/p) / d(q/p,u',v',u,v)
496 
497  // da_x/du'
498  J_pM[3][1] = fSpu/pTildeMag*(u.X()-pTilde.X()/(pTildeMag*pTildeMag)*u*pTilde);
499  J_pM[4][1] = fSpu/pTildeMag*(u.Y()-pTilde.Y()/(pTildeMag*pTildeMag)*u*pTilde);
500  J_pM[5][1] = fSpu/pTildeMag*(u.Z()-pTilde.Z()/(pTildeMag*pTildeMag)*u*pTilde);
501  // da_x/dv'
502  J_pM[3][2] = fSpu/pTildeMag*(v.X()-pTilde.X()/(pTildeMag*pTildeMag)*v*pTilde);
503  J_pM[4][2] = fSpu/pTildeMag*(v.Y()-pTilde.Y()/(pTildeMag*pTildeMag)*v*pTilde);
504  J_pM[5][2] = fSpu/pTildeMag*(v.Z()-pTilde.Z()/(pTildeMag*pTildeMag)*v*pTilde);
505  // dqOp/dqOp
506  J_pM[6][0] = 1.;
507 
508  TMatrixT<Double_t> J_pM_transp(J_pM);
509  J_pM_transp.T();
510 
511  cov7x7 = J_pM*(fCov*J_pM_transp);
512  if (cov7x7[0][0]>=1000. || cov7x7[0][0]<1.E-50)
513  {
514  if (pOut) {
515  (*pOut) << "RKTrackRep::extrapolate(): cov7x7[0][0] is crazy. Rescale off-diags. Try again. fCov, cov7x7 were: " << std::endl;
516  PrintROOTobject(*pOut, fCov);
517  PrintROOTobject(*pOut, cov7x7);
518  }
519  rescaleCovOffDiags();
520  cov7x7 = J_pM*(fCov*J_pM_transp);
521  if (pOut) {
522  (*pOut) << "New cov7x7 and fCov are ... " << std::endl;
523  PrintROOTobject(*pOut, cov7x7);
524  PrintROOTobject(*pOut, fCov);
525  }
526  }
527 
528 
529  TVector3 pos = o + fState[3][0]*u + fState[4][0]*v;
530  TMatrixT<Double_t> state7(7,1);
531  state7[0][0] = pos.X();
532  state7[1][0] = pos.Y();
533  state7[2][0] = pos.Z();
534  state7[3][0] = pTilde.X()/pTildeMag;;
535  state7[4][0] = pTilde.Y()/pTildeMag;;
536  state7[5][0] = pTilde.Z()/pTildeMag;;
537  state7[6][0] = fState[0][0];
538 
539  double coveredDistance = this->Extrap(pl,&state7,&cov7x7);
540 
541 
542  TVector3 O = pl.getO();
543  TVector3 U = pl.getU();
544  TVector3 V = pl.getV();
545  TVector3 W = pl.getNormal();
546 
547  double X = state7[0][0];
548  double Y = state7[1][0];
549  double Z = state7[2][0];
550  double AX = state7[3][0];
551  double AY = state7[4][0];
552  double AZ = state7[5][0];
553  double QOP = state7[6][0];
554  TVector3 A(AX,AY,AZ);
555  TVector3 Point(X,Y,Z);
556  TMatrixT<Double_t> J_Mp(5,7);
557 
558  // J_Mp matrix is d(q/p,u',v',u,v) / d(x,y,z,ax,ay,az,q/p)
559  J_Mp[0][6] = 1.;
560  //du'/da_x
561  double AtW = A*W;
562  J_Mp[1][3] = (U.X()*(AtW)-W.X()*(A*U))/(AtW*AtW);
563  J_Mp[1][4] = (U.Y()*(AtW)-W.Y()*(A*U))/(AtW*AtW);
564  J_Mp[1][5] = (U.Z()*(AtW)-W.Z()*(A*U))/(AtW*AtW);
565  //dv'/da_x
566  J_Mp[2][3] = (V.X()*(AtW)-W.X()*(A*V))/(AtW*AtW);
567  J_Mp[2][4] = (V.Y()*(AtW)-W.Y()*(A*V))/(AtW*AtW);
568  J_Mp[2][5] = (V.Z()*(AtW)-W.Z()*(A*V))/(AtW*AtW);
569  //du/dx
570  J_Mp[3][0] = U.X();
571  J_Mp[3][1] = U.Y();
572  J_Mp[3][2] = U.Z();
573  //dv/dx
574  J_Mp[4][0] = V.X();
575  J_Mp[4][1] = V.Y();
576  J_Mp[4][2] = V.Z();
577 
578  TMatrixT<Double_t> J_Mp_transp(J_Mp);
579  J_Mp_transp.T();
580 
581  covPred.ResizeTo(5,5);
582  covPred = J_Mp*(cov7x7*J_Mp_transp);
583 
584 
585  statePred.ResizeTo(5,1);
586  statePred[0][0] = QOP;
587  statePred[1][0] = (A*U)/(A*W);
588  statePred[2][0] = (A*V)/(A*W);
589  statePred[3][0] = (Point-O)*U;
590  statePred[4][0] = (Point-O)*V;
591 
592  fCachePlane = pl;
593  fCacheSpu = (A*W)/fabs(A*W);
594 
595  return coveredDistance;
596 }
genf::RKTrackRep::fCachePlane
GFDetPlane fCachePlane
Definition: RKTrackRep.h:167
E
process_name E
Definition: rawtorootconvert.fcl:26
X
then echo echo For and will not be changed by echo further linking echo echo B echo The symbol is in the uninitialized data multiple common symbols may appear with the echo same name If the symbol is defined the common echo symbols are treated as undefined references For more echo details on common see the discussion of warn common echo in *Note Linker see the discussion of warn common echo in *Note Linker such as a global int variable echo as opposed to a large global array echo echo I echo The symbol is an indirect reference to another symbol This echo is a GNU extension to the a out object file format which is echo rarely used echo echo N echo The symbol is a debugging symbol echo echo R echo The symbol is in a read only data section echo echo S echo The symbol is in an uninitialized data section for small echo objects echo echo T echo The symbol is in the the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo W echo The symbol is a weak symbol that has not been specifically echo tagged as a weak object symbol When a weak defined symbol echo is linked with a normal defined the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo echo The symbol is a stabs symbol in an a out object file In echo this the next values printed are the stabs other echo the stabs desc and the stab type Stabs symbols are echo used to hold debugging information For more echo see *Note or object file format specific echo echo For Mac OS X
Definition: find_global_symbol.sh:107
genf::GFAbsTrackRep::fCov
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:94
V
BEGIN_PROLOG V
Definition: trackfinderalgorithms.fcl:26
genf::GFDetPlane::getO
TVector3 getO() const
Definition: GFDetPlane.h:81
dcel2d::Point
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:44
genf::RKTrackRep::rescaleCovOffDiags
void rescaleCovOffDiags()
Definition: RKTrackRep.cxx:1222
gen_crt_frags.w
tuple w
Definition: gen_crt_frags.py:1416
genf::GFDetPlane::getU
TVector3 getU() const
Definition: GFDetPlane.h:82
genf::RKTrackRep::Extrap
double Extrap(const GFDetPlane &plane, TMatrixT< Double_t > *state, TMatrixT< Double_t > *cov=NULL) const
Handles propagation and material effects.
Definition: RKTrackRep.cxx:1058
genf::PrintROOTobject
void PrintROOTobject(std::ostream &, const ROOTOBJ &)
Small utility functions which print some ROOT objects into an output stream.
Definition: GFException.h:127
genf::GFAbsTrackRep::fState
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:91
genf::GFDetPlane::getV
TVector3 getV() const
Definition: GFDetPlane.h:83
lib.dedx.A
float A
Definition: dedx.py:137
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 601 of file RKTrackRep.cxx.

602  {
603 
604  TVector3 o=fRefPlane.getO();
605  TVector3 u=fRefPlane.getU();
606  TVector3 v=fRefPlane.getV();
607  TVector3 w=u.Cross(v);
608 
609 
610  TVector3 pTilde = fSpu* (w + fState[1][0] * u + fState[2][0] * v);
611  double pTildeMag = pTilde.Mag();
612 
613  TVector3 pos = o + fState[3][0]*u + fState[4][0]*v;
614 
615  TMatrixT<Double_t> state7(7,1);
616  state7[0][0] = pos.X();
617  state7[1][0] = pos.Y();
618  state7[2][0] = pos.Z();
619  state7[3][0] = pTilde.X()/pTildeMag;
620  state7[4][0] = pTilde.Y()/pTildeMag;
621  state7[5][0] = pTilde.Z()/pTildeMag;
622  state7[6][0] = fState[0][0];
623 
624  TVector3 O = pl.getO();
625  TVector3 U = pl.getU();
626  TVector3 V = pl.getV();
627  TVector3 W = pl.getNormal();
628 
629  double coveredDistance = this->Extrap(pl,&state7);
630 
631  double X = state7[0][0];
632  double Y = state7[1][0];
633  double Z = state7[2][0];
634  double AX = state7[3][0];
635  double AY = state7[4][0];
636  double AZ = state7[5][0];
637  double QOP = state7[6][0];
638  TVector3 A(AX,AY,AZ);
639  TVector3 Point(X,Y,Z);
640 
641  statePred.ResizeTo(5,1);
642  statePred[0][0] = QOP;
643  statePred[1][0] = (A*U)/(A*W);
644  statePred[2][0] = (A*V)/(A*W);
645  statePred[3][0] = (Point-O)*U;
646  statePred[4][0] = (Point-O)*V;
647 
648  return coveredDistance;
649 }
X
then echo echo For and will not be changed by echo further linking echo echo B echo The symbol is in the uninitialized data multiple common symbols may appear with the echo same name If the symbol is defined the common echo symbols are treated as undefined references For more echo details on common see the discussion of warn common echo in *Note Linker see the discussion of warn common echo in *Note Linker such as a global int variable echo as opposed to a large global array echo echo I echo The symbol is an indirect reference to another symbol This echo is a GNU extension to the a out object file format which is echo rarely used echo echo N echo The symbol is a debugging symbol echo echo R echo The symbol is in a read only data section echo echo S echo The symbol is in an uninitialized data section for small echo objects echo echo T echo The symbol is in the the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo W echo The symbol is a weak symbol that has not been specifically echo tagged as a weak object symbol When a weak defined symbol echo is linked with a normal defined the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo echo The symbol is a stabs symbol in an a out object file In echo this the next values printed are the stabs other echo the stabs desc and the stab type Stabs symbols are echo used to hold debugging information For more echo see *Note or object file format specific echo echo For Mac OS X
Definition: find_global_symbol.sh:107
V
BEGIN_PROLOG V
Definition: trackfinderalgorithms.fcl:26
genf::GFDetPlane::getO
TVector3 getO() const
Definition: GFDetPlane.h:81
dcel2d::Point
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:44
gen_crt_frags.w
tuple w
Definition: gen_crt_frags.py:1416
genf::GFDetPlane::getU
TVector3 getU() const
Definition: GFDetPlane.h:82
genf::RKTrackRep::Extrap
double Extrap(const GFDetPlane &plane, TMatrixT< Double_t > *state, TMatrixT< Double_t > *cov=NULL) const
Handles propagation and material effects.
Definition: RKTrackRep.cxx:1058
genf::GFAbsTrackRep::fState
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:91
genf::GFDetPlane::getV
TVector3 getV() const
Definition: GFDetPlane.h:83
lib.dedx.A
float A
Definition: dedx.py:137
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 424 of file RKTrackRep.cxx.

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

367  {
368 
369  static const int maxIt(30);
370 
371  TVector3 o=fRefPlane.getO();
372  TVector3 u=fRefPlane.getU();
373  TVector3 v=fRefPlane.getV();
374  TVector3 w=u.Cross(v);
375 
376  TVector3 pTilde = fSpu* (w + fState[1][0] * u + fState[2][0] * v);
377  pTilde.SetMag(1.);
378 
379  TVector3 point = o + fState[3][0]*u + fState[4][0]*v;
380 
381  TMatrixT<Double_t> state7(7,1);
382  state7[0][0] = point.X();
383  state7[1][0] = point.Y();
384  state7[2][0] = point.Z();
385  state7[3][0] = pTilde.X();
386  state7[4][0] = pTilde.Y();
387  state7[5][0] = pTilde.Z();
388  state7[6][0] = fState[0][0];
389 
390  double coveredDistance(0.);
391 
392  GFDetPlane pl;
393  int iterations(0);
394 
395  while(true){
396  pl.setON(pos,TVector3(state7[3][0],state7[4][0],state7[5][0]));
397  coveredDistance = this->Extrap(pl,&state7);
398 
399  if(fabs(coveredDistance)<MINSTEP) break;
400  if(++iterations == maxIt) {
401  throw GFException("RKTrackRep::extrapolateToPoint==> extrapolation to point failed, maximum number of iterations reached",__LINE__,__FILE__).setFatal();
402  }
403  }
404  poca.SetXYZ(state7[0][0],state7[1][0],state7[2][0]);
405  dirInPoca.SetXYZ(state7[3][0],state7[4][0],state7[5][0]);
406 }
genf::GFDetPlane::getO
TVector3 getO() const
Definition: GFDetPlane.h:81
MINSTEP
#define MINSTEP
Definition: RKTrackRep.cxx:39
gen_crt_frags.w
tuple w
Definition: gen_crt_frags.py:1416
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:48
genf::GFDetPlane::getU
TVector3 getU() const
Definition: GFDetPlane.h:82
genf::RKTrackRep::Extrap
double Extrap(const GFDetPlane &plane, TMatrixT< Double_t > *state, TMatrixT< Double_t > *cov=NULL) const
Handles propagation and material effects.
Definition: RKTrackRep.cxx:1058
GFException::setFatal
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:78
genf::GFAbsTrackRep::fState
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:91
genf::GFDetPlane::getV
TVector3 getV() const
Definition: GFDetPlane.h:83
const TMatrixT< double > * genf::RKTrackRep::getAuxInfo ( const GFDetPlane &  pl)

Definition at line 55 of file RKTrackRep.cxx.

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

Returns charge.

Implements genf::GFAbsTrackRep.

Definition at line 142 of file RKTrackRep.h.

142 {return fCharge;}
genf::RKTrackRep::fCharge
double fCharge
Charge.
Definition: RKTrackRep.h:182
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 323 of file RKTrackRep.cxx.

323  {
324  TMatrixT<Double_t> statePred(fState);
325  TVector3 retmom;
326  if(pl!=fRefPlane) {
327  extrapolate(pl,statePred);
328  retmom = fCacheSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
329  }
330  else{
331  retmom = fSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
332  }
333  retmom.SetMag(1./fabs(statePred[0][0]));
334  return retmom;
335 }
genf::RKTrackRep::extrapolate
double extrapolate(const GFDetPlane &, TMatrixT< Double_t > &statePred, TMatrixT< Double_t > &covPred)
returns the tracklength spanned in this extrapolation
Definition: RKTrackRep.cxx:475
genf::GFAbsTrackRep::fState
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:91
TVector3 genf::RKTrackRep::getMomLast ( const GFDetPlane &  pl)

Definition at line 337 of file RKTrackRep.cxx.

337  {
338  TMatrixT<Double_t> statePred(fLastState);
339  TVector3 retmom;
340  retmom = fSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
341 
342  retmom.SetMag(1./fabs(statePred[0][0]));
343  return retmom;
344 }
genf::GFAbsTrackRep::fLastState
TMatrixT< Double_t > fLastState
Definition: GFAbsTrackRep.h:109
int genf::RKTrackRep::getPDG ( )

Definition at line 311 of file RKTrackRep.cxx.

311 {return fPdg;}
genf::RKTrackRep::fPdg
int fPdg
PDG particle code.
Definition: RKTrackRep.h:178
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 313 of file RKTrackRep.cxx.

313  {
314  if(pl!=fRefPlane){
315  TMatrixT<Double_t> s(5,1);
316  extrapolate(pl,s);
317  return pl.getO()+s[3][0]*pl.getU()+s[4][0]*pl.getV();
318  }
319  return fRefPlane.getO()+fState[3][0]*fRefPlane.getU()+fState[4][0]*fRefPlane.getV();
320 }
genf::RKTrackRep::extrapolate
double extrapolate(const GFDetPlane &, TMatrixT< Double_t > &statePred, TMatrixT< Double_t > &covPred)
returns the tracklength spanned in this extrapolation
Definition: RKTrackRep.cxx:475
genf::GFDetPlane::getO
TVector3 getO() const
Definition: GFDetPlane.h:81
genf::GFDetPlane::getU
TVector3 getU() const
Definition: GFDetPlane.h:82
s
then echo File list $list not found else cat $list while read file do echo $file sed s
Definition: file_to_url.sh:60
genf::GFAbsTrackRep::fState
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:91
genf::GFDetPlane::getV
TVector3 getV() const
Definition: GFDetPlane.h:83
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 347 of file RKTrackRep.cxx.

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

Definition at line 163 of file RKTrackRep.h.

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

Definition at line 173 of file RKTrackRep.h.

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

Definition at line 411 of file RKTrackRep.cxx.

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

Implements genf::GFAbsTrackRep.

Definition at line 79 of file RKTrackRep.h.

79 {return new RKTrackRep();}
void genf::RKTrackRep::rescaleCovOffDiags ( )

Definition at line 1222 of file RKTrackRep.cxx.

1223 {
1224 
1225  for(int i=0;i<fCov.GetNrows();++i){
1226  for(int j=0;j<fCov.GetNcols();++j){
1227  if(i!=j){//off diagonal
1228  fCov[i][j]=0.5*fCov[i][j];
1229  }
1230  else{//diagonal. Do not let diag cov element be <=0.
1231  if (fCov[i][j]<=0.0) fCov[i][j] = 0.01;
1232  }
1233  }
1234  }
1235 }
genf::GFAbsTrackRep::fCov
TMatrixT< Double_t > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:94
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 691 of file RKTrackRep.cxx.

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

42  {
43  if(aux != NULL) {
44  fCacheSpu = (*aux)(0,0);
45  } else {
46  if(pl!=fCachePlane){
47  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();
48  }
49  }
50  GFAbsTrackRep::setData(st,pl,cov);
51  if (fCharge*fState[0][0] < 0) fCharge *= -1; // set charge accordingly! (fState[0][0] = q/p)
52  fSpu=fCacheSpu;
53 }
genf::RKTrackRep::fCachePlane
GFDetPlane fCachePlane
Definition: RKTrackRep.h:167
genf::GFAbsTrackRep::setData
virtual void setData(const TMatrixT< Double_t > &st, const GFDetPlane &pl, const TMatrixT< Double_t > *cov=NULL)
Definition: GFAbsTrackRep.h:256
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:48
GFException::setFatal
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:78
genf::GFAbsTrackRep::fState
TMatrixT< Double_t > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:91
genf::RKTrackRep::fCharge
double fCharge
Charge.
Definition: RKTrackRep.h:182
void genf::RKTrackRep::setPDG ( int  i)

Set PDG particle code.

Definition at line 299 of file RKTrackRep.cxx.

299  {
300  fPdg = i;
301  TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(fPdg);
302  if(part == 0){
303  std::cerr << "RKTrackRep::setPDG particle " << i
304  << " not known to TDatabasePDG -> abort" << std::endl;
305  exit(1);
306  }
307  fMass = part->Mass();
308  fCharge = part->Charge()/(3.);
309 }
genf::RKTrackRep::fPdg
int fPdg
PDG particle code.
Definition: RKTrackRep.h:178
cerr
BEGIN_PROLOG could also be cerr
Definition: messageservice.fcl:10
exit
exit
Definition: generate_icarus_geometry_splitwire.pl:430
genf::RKTrackRep::fMass
double fMass
Mass (in GeV)
Definition: RKTrackRep.h:180
genf::RKTrackRep::fCharge
double fCharge
Charge.
Definition: RKTrackRep.h:182
void genf::RKTrackRep::switchDirection ( )
inlinevirtual

deprecated

Implements genf::GFAbsTrackRep.

Definition at line 144 of file RKTrackRep.h.

144 {fDirection = (!fDirection);}

Member Data Documentation

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

Definition at line 170 of file RKTrackRep.h.

GFDetPlane genf::RKTrackRep::fCachePlane
private

Definition at line 167 of file RKTrackRep.h.

double genf::RKTrackRep::fCacheSpu
private

Definition at line 168 of file RKTrackRep.h.

double genf::RKTrackRep::fCharge
private

Charge.

Definition at line 182 of file RKTrackRep.h.

bool genf::RKTrackRep::fDirection
private

Definition at line 175 of file RKTrackRep.h.

double genf::RKTrackRep::fMass
private

Mass (in GeV)

Definition at line 180 of file RKTrackRep.h.

int genf::RKTrackRep::fPdg
private

PDG particle code.

Definition at line 178 of file RKTrackRep.h.

double genf::RKTrackRep::fSpu
private

Definition at line 169 of file RKTrackRep.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.5