StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StiKalmanTrackNode.h
1 #ifndef StiKalmanTrackNode_H
2 #define StiKalmanTrackNode_H 1
3 #define STI_NODE_DEBUG
4 
5 #include <Stiostream.h>
6 #include <stdlib.h>
7 #include <stdexcept>
8 #include <math.h>
9 #include "StiTrackNode.h"
10 #include "StThreeVector.hh"
11 #include "StThreeVectorF.hh"
12 #include "StDetectorDbMaker/StiKalmanTrackFinderParameters.h"
13 #include "StiShape.h"
14 #include "StiPlanarShape.h"
15 #include "StiCylindricalShape.h"
16 #include "StiPlacement.h"
17 #include "StiHit.h"
18 #include "StiMaterial.h"
19 class StiDetector;
20 class StiMaterial;
21 //class StiElossCalculator;
22 
23 typedef enum {
24  kFailed = -1, // could not find intersection
25  kHit,
26  kEdgePhiPlus, kEdgeZminus, kEdgePhiMinus, kEdgeZplus,
27  kMissPhiPlus, kMissZminus, kMissPhiMinus, kMissZplus,
28  kEnded
29 } StiIntersection;
30 
31 class StiNodeStat {
32 public:
33  StiNodeStat(){reset();}
34 void reset(){memset(this,0,sizeof(StiNodeStat));}
35  double x1; // local x position of current node
36  double y1; // local y position of current node
37  double cosCA1,sinCA1; // cos & sin of local pfi angle
38  double x2,y2; // x & y of track in the next node but in current local
39  double cosCA2,sinCA2; // cos & sin in this point
40  double dx, dy; // x2-x1 & y2-y1
41  double dl0; // linear distance
42  double dl; // curved distance
43  double sumSin, sumCos;
44 };
45 
46 class StiNodeExt {
47 public:
48 void reset(){mPP.reset();mPE.reset();mMtx.reset();}
49 void unset(){;}
50 public:
51  StiNodePars mPP; //Predicted params
52  StiNodeMtx mMtx;
53  StiNodeErrs mPE; //Predicted params errs
54 };
55 
56 class StiNodeInf {
57 public:
58 void reset(){mPP.reset();mPE.reset();mHrr.reset();}
59 void unset(){;}
60 public:
61  StiNodePars mPP;
62  StiNodeErrs mPE;
63  StiHitErrs mHrr;
64 };
65 
66 
67 
68 
69 
82 {
83 friend class StiTrackNodeHelper;
84 public:
85  StiKalmanTrackNode(){_inf=0; _ext=0; reset();}
87  virtual ~StiKalmanTrackNode(){reduce();mId=-1;};
88  const StiKalmanTrackNode& operator=(const StiKalmanTrackNode &node);
89 
90  static double mcs2(double relRadThickness, double beta2, double p2) {return 14.1*14.1*relRadThickness/(beta2*p2*1e6);}
92  void reset();
93  void unset();
95  void resetError(double fak=0);
97  void initialize(StiHit*h);
98  void initialize(StiDetector *d);
99  void initialize(const double dir[3]); //dir[0]**2+dir[1]**2 ==1
100 //void initialize(StiHit*h,double alpha, double eta, double curvature, double tanl);
101 
103  void setState(const StiKalmanTrackNode * node);
105  void get(double& alpha, double& xRef, double x[kNPars], double cc[kNErrs], double& chi2);
106 
108  void getGlobalRadial(double x[6],double e[15]);
109 
111  void getGlobalTpt (float x[6],float e[15]);
112 
114  int getCharge() const {return (mFP.ptin() > 0) ? -1 : 1;}
120  StThreeVector<double> getMomentum() const;
121  StThreeVector<double> getGlobalMomentum() const;
124  void getMomentum(double p[3], double e[6]=0) const;
126  double getCurvature() const {return mFP.curv();}
127  void setCurvature(double curvature) {mFP.curv()=curvature;}
128  double getDipAngle() const {return atan(mFP.tanl());}
129  double getTanL() const {return mFP.tanl();}
131  double getPt() const;
133  double getP() const {return (getPt()*::sqrt(1.+mFP.tanl()*mFP.tanl()));}
136  double getHz() const;
137  double getField() const {return getHz();}
138  double x_g() const;
139  double y_g() const;
140  double z_g() const;
141  double getX() const { return mFP.x();}
142  double getY() const { return mFP.y();}
143  double getZ() const { return mFP.z();}
144  double x() const { return mFP.x();}
145  double y() const { return mFP.y();}
146  double z() const { return mFP.z();}
147  double getRxy() const { return sqrt(mFP.x()*mFP.x()+mFP.y()*mFP.y());}
148 
149  double getEta () const {return mFP.eta(); }
150  double getSin () const {return mFP._sinCA;}
151  double getCos () const {return mFP._cosCA;}
152  double getAlpha() const {return _alpha; }
153  const double *hitErrs() const {return mHrr.G(); }
154  double getEyy() const {return mHrr.hYY;}
155  double getEzz() const {return mHrr.hZZ;}
156  double getCyy() const {return mFE._cYY;}
157  double getCzz() const {return mFE._cZZ;}
158  double const *getPars()const {return mFP.P;}
159  int getHitCount () const {return hitCount;}
160  int getNullCount() const {return nullCount;}
161  int getContigHitCount () const {return contiguousHitCount ;}
162  int getContigNullCount() const {return contiguousNullCount;}
163  //const char &getHitCount () const {return hitCount;}
164  //const char &getNullCount() const {return nullCount;}
165  //const char &getContigHitCount () const {return contiguousHitCount ;}
166  //const char &getContigNullCount() const {return contiguousNullCount;}
167  int incHitCount () {return ++hitCount;}
168  int incNullCount() {return ++nullCount;}
169  int incContigHitCount () {return ++contiguousHitCount ;}
170  int incContigNullCount() {return ++contiguousNullCount;}
171  void setHitCount (char c=0) { hitCount=c;}
172  void setNullCount(char c=0) { nullCount=c;}
173  void setContigHitCount (char c=0) { contiguousHitCount=c ;}
174  void setContigNullCount(char c=0) { contiguousNullCount=c;}
175  double getTime() const;
176 
177  void setHitCand(int nhits) {mHitCand = nhits;}
178  void setIHitCand(int ihit) {mIHitCand = ihit;}
179  int getHitCand()const {return mHitCand;}
180  int getIHitCand()const {return mIHitCand;}
181  StiELoss *getELoss() {return mELoss;}
182 const StiELoss *getELoss()const {return mELoss;}
183  static void PrintStep();
184  StThreeVector<double>getPoint() const;
185  StThreeVector<double>getGlobalPoint() const;
187  void getGlobalMomentum(double p[3], double e[6]=0) const;
188  int isEnded() const;
189  int isDca() const;
190 
192  int propagate(StiKalmanTrackNode *p, const StiDetector * tDet, int dir); //throw (Exception);
193 
195  bool propagate(const StiKalmanTrackNode *p, StiHit * vertex, int dir);
196 
197  bool propagateToBeam(const StiKalmanTrackNode *p, int dir);
198  int propagateToRadius(StiKalmanTrackNode *pNode, double radius,int dir);
199  int locate();
200  int propagate(double x,int option,int dir);
201  void propagateMtx();
202  void propagateError();
203  void saveInfo(int kase=1);
204 const StiNodeInf *getInfo() const {return _inf;}
205  int testError(double *emx,int begend);
206  void numeDeriv(double val,int kind,int shape=0,int dir=0);
207  int testDeriv(double *der);
208  void propagateMCS(StiKalmanTrackNode * previousNode, const StiDetector * tDet);
209  int nudge(StiHit *hit=0);
210  double evaluateChi2 (const StiHit *hit);
211  double evaluateChi2Info(const StiHit *hit) const;
212  int updateNode();
213  int rotate(double alpha);
214  int getHelicity()const {return (mFP.curv() < 0) ? -1 : 1;}
215  double getPhase() const;
216  double getPsi() const;
217  double getWindowY();
218  double getWindowZ();
219  double pitchAngle() const {return atan(mFP.tanl());}
220  double crossAngle() const {return asin(mFP._sinCA);}
221  double sinCrossAngle() const {return mFP._sinCA;}
222  double pathlength() const;
223  double pathLToNode(const StiKalmanTrackNode * const oNode);
224  StThreeVectorD* getLengths(StiKalmanTrackNode *nextNode);
225 
226  double length(const StThreeVector<double>& delta, double curv);
227  double getDedx() const;
228  static double nice(double angle);
231  void setHitErrors(const StiHit *hit=0);
232  StiHitErrs getGlobalHitErrs(const StiHit *hit) const;
233  friend ostream& operator<<(ostream& os, const StiKalmanTrackNode& n);
234 
235  double getX0() const;
236  double getGasX0() const;
237  double getDensity() const;
238  double getGasDensity() const;
239 
240  void extend();
241  void reduce();
242  static Int_t debug() {return _debug;}
243  static void setDebug(Int_t m) {_debug = m;}
244  static void SetLaser(Int_t m) {_laser = m;}
245  static Int_t IsLaser() {return _laser;}
246  void PrintpT(const Char_t *opt="") const ;
247  int getFlipFlop() const {return mFlipFlop;}
248  static void ResetComment(const Char_t *m = "") {comment = m; commentdEdx = "";}
249  static const Char_t *Comment() {return comment.Data();}
251  int print(const char *opt) const;
252 
253  private:
254  void extinf(); //add inf block
255  void static saveStatics(double *sav);
256  void static backStatics(double *sav);
257  static StiNodeExt *nodeExtInstance();
258  static StiNodeInf *nodeInfInstance();
259  void propagateCurv(const StiKalmanTrackNode *parent);
260 
261 // Extended members
262  public:
263 
264  const StiNodePars &fitPars() const {return mFP; }
265  StiNodePars &fitPars() {return mFP; }
266  const StiNodeErrs &fitErrs() const {return mFE; }
267  StiNodeErrs &fitErrs() {return mFE; }
268  const StiNodePars &mPP() const {return _ext->mPP; }
269  StiNodePars &mPP() {if (!_ext) extend(); return _ext->mPP; }
270  const StiNodeErrs &mPE() const {return _ext->mPE; }
271  StiNodeErrs &mPE() {if (!_ext) extend(); return _ext->mPE; }
272  const StiNodeMtx &mMtx()const {return _ext->mMtx;}
273  StiNodeMtx &mMtx() {if (!_ext) extend();return _ext->mMtx;}
274  const StiNode2Pars &unTouched() const {return mUnTouch;}
275  void setUntouched();
276  protected:
277 
278  char _beg[1];
279  double _alpha;
281  mutable double mHz;
282  StiNodePars mFP;
285  StiNode2Pars mUnTouch;
286  StiHitErrs mHrr;
287  StiELoss mELoss[2];
288  char hitCount;
289  char nullCount;
290  char contiguousHitCount;
291  char contiguousNullCount;
292  char mFlipFlop;
293  char mHitCand;
294  char mIHitCand;
295  char _end[1];
296  StiNodeExt *_ext;
297  StiNodeInf *_inf;
298  static StiNodeStat mgP;
299  static bool useCalculatedHitError;
300 // debug variables
301  static int fDerivTestOn;
302  static double fDerivTest[kNPars][kNPars];
303  static int _debug;
304  static TString comment;
305  static TString commentdEdx;
306  static int _laser;
307 public:
308  int mId; //for debug only
309 };
310 
311 
312 inline double StiKalmanTrackNode::nice(double angle)
313 {
314  if (angle <= -M_PI) angle += 2*M_PI;
315  if (angle > M_PI) angle -= 2*M_PI;
316  return angle;
317 }
318 
319 inline StThreeVector<double> StiKalmanTrackNode::getMomentum() const
320 {
321  double pt = getPt();
322  return StThreeVector<double>(pt*mFP._cosCA,pt*mFP._sinCA,pt*mFP.tanl());
323 }
324 
326 {
327  double pt = getPt();
328  return StThreeVectorF(pt*mFP._cosCA,pt*mFP._sinCA,pt*mFP.tanl());
329 }
330 
331 inline StThreeVector<double> StiKalmanTrackNode::getGlobalMomentum() const
332 {
333  StThreeVector<double> p = getMomentum();
334  p.rotateZ(_alpha);
335  return p;
336 }
337 
339 {
341  p.rotateZ(_alpha);
342  return p;
343 }
344 
345 
346 //stl helper functor
347 
349 {
350  bool operator()(const StiKalmanTrackNode& lhs, const StiKalmanTrackNode& rhs) const;
351 };
352 
353 struct StreamX
354 {
355  void operator()(const StiKalmanTrackNode& node)
356  {
357  cout <<node.getX()<<endl;
358  }
359 };
360 
363 inline double StiKalmanTrackNode::getX0() const
364 {
365  const StiDetector * det = getDetector();
366  if (!det)
367  return 0.;
368  return det->getMaterial()->getX0();
369 }
370 
373 inline double StiKalmanTrackNode::getGasX0() const
374 {
375  const StiDetector * det = getDetector();
376  if (!det)
377  return 0.;
378  return det->getGas()->getX0();
379 }
380 
381 inline double StiKalmanTrackNode::getDensity() const
382 {
383  const StiDetector * det = getDetector();
384  if (!det)
385  return 0.;
386  return det->getMaterial()->getDensity();
387 }
388 
389 inline double StiKalmanTrackNode::getGasDensity() const
390 {
391  const StiDetector * det = getDetector();
392  if (!det)
393  return 0.;
394  return det->getGas()->getDensity();
395 }
396 
397 
398 inline StThreeVectorD* StiKalmanTrackNode::getLengths(StiKalmanTrackNode* nextNode)
399 {
400  double x1=pathlength()/2.;
401  double x3=nextNode->pathlength()/2.;
402  double x2=pathLToNode(nextNode);
403  if (x2> (x1+x3)) x2=x2-x1-x3;
404  else x2=0;
405 
406  return new StThreeVectorD(x1/getX0(),
407  x2/getDetector()->getMaterial()->getX0(),
408  x3/nextNode->getX0());
409 }
410 
411 inline double StiKalmanTrackNode::getDedx() const
412 {
413 
414  StiHit *hit = getHit();
415  if (!hit) return -1;
416  double de=hit->getEloss();
417  double dx=pathlength();
418  if(dx>0 && de>0) return de/dx;
419  return -1;
420 }
421 
422 #endif
423 
void setState(const StiKalmanTrackNode *node)
Sets the Kalman state of this node equal to that of the given node.
int propagate(StiKalmanTrackNode *p, const StiDetector *tDet, int dir)
Propagates a track encapsulated by the given node &quot;p&quot; to the given detector &quot;tDet&quot;.
double getX0() const
double getP() const
Calculates and returns the momentum of the track at this node.
int propagateToRadius(StiKalmanTrackNode *pNode, double radius, int dir)
void getGlobalRadial(double x[6], double e[15])
Extract state information from this node in Radial representation.
bool propagateToBeam(const StiKalmanTrackNode *p, int dir)
Definition: StiHit.h:51
StThreeVectorF getGlobalMomentumF() const
void reset()
Resets the node to a &quot;null&quot; un-used state.
double pathLToNode(const StiKalmanTrackNode *const oNode)
int rotate(double alpha)
double getGasX0() const
int print(const char *opt) const
rotation angle of local coordinates wrt global coordinates
StThreeVector< double > getHelixCenter() const
Return center of helix circle in global coordinates.
double getPt() const
Calculates and returns the transverse momentum of the track at this node.
friend ostream & operator<<(ostream &os, const StiKalmanTrackNode &n)
StThreeVectorF getMomentumF() const
Convenience Method that returns the track momentum at this node.
Float_t getEloss()
Return the energy deposition associated with this point.
Definition: StiHit.cxx:253
double _cosCA
sine and cosine of cross angle
Definition: StiNodePars.h:51
void resetError(double fak=0)
Resets errors for refit.
double evaluateChi2(const StiHit *hit)
double getCurvature() const
Calculates and returns the tangent of the track pitch angle at this node.
void getGlobalTpt(float x[6], float e[15])
Extract state information from this node in TPT representation.
double getX0() const
Get the radiation length in centimeter.
Definition: StiMaterial.h:37
double mHz
Z mag field in units PGev = Hz*Rcm.
void propagateMCS(StiKalmanTrackNode *previousNode, const StiDetector *tDet)
StiNodeErrs mFE
covariance matrix of the track parameters
int getCharge() const
Get the charge (sign) of the track at this node.
double evaluateChi2Info(const StiHit *hit) const
double getDensity() const
Get the material density in grams/cubic centimeter.
Definition: StiMaterial.h:35
void initialize(StiHit *h)
Initialize this node with the given hit information.