00001 #ifndef StiKalmanTrackNode_H
00002 #define StiKalmanTrackNode_H 1
00003 #define STI_NODE_DEBUG
00004
00005 #include <Stiostream.h>
00006 #include <stdlib.h>
00007 #include <stdexcept>
00008 #include <math.h>
00009 #include "StiTrackNode.h"
00010 #include "StThreeVector.hh"
00011 #include "StThreeVectorF.hh"
00012 #include "StDetectorDbMaker/StiKalmanTrackFinderParameters.h"
00013 #include "StiShape.h"
00014 #include "StiPlanarShape.h"
00015 #include "StiCylindricalShape.h"
00016 #include "StiPlacement.h"
00017 #include "StiHit.h"
00018 #include "StiMaterial.h"
00019 class StiDetector;
00020 class StiMaterial;
00021 class StiElossCalculator;
00022
00023 typedef enum {
00024 kFailed = -1,
00025 kHit,
00026 kEdgePhiPlus, kEdgeZminus, kEdgePhiMinus, kEdgeZplus,
00027 kMissPhiPlus, kMissZminus, kMissPhiMinus, kMissZplus,
00028 kEnded
00029 } StiIntersection;
00030
00031 class StiNodeStat {
00032 public:
00033 StiNodeStat(){reset();}
00034 void reset(){memset(this,0,sizeof(StiNodeStat));}
00035 double dx, dy, dl0, dl;
00036 double x1,y1,cosCA1,sinCA1;
00037 double x2,y2,cosCA2,sinCA2;
00038 double sumSin, sumCos;
00039 };
00040
00041 class StiNodeExt {
00042 public:
00043 void reset(){mPP.reset();mPE.reset();mMtx.reset();}
00044 void unset(){;}
00045 public:
00046 StiNodePars mPP;
00047 StiNodeMtx mMtx;
00048 StiNodeErrs mPE;
00049 };
00050
00051 class StiNodeInf {
00052 public:
00053 void reset(){mPP.reset();mPE.reset();mHrr.reset();}
00054 void unset(){;}
00055 public:
00056 StiNodePars mPP;
00057 StiNodeErrs mPE;
00058 StiHitErrs mHrr;
00059 };
00060
00061
00062
00063
00064
00076 class StiKalmanTrackNode : public StiTrackNode
00077 {
00078 friend class StiTrackNodeHelper;
00079 public:
00080 StiKalmanTrackNode(){_inf=0; _ext=0; reset();}
00081 StiKalmanTrackNode(const StiKalmanTrackNode &node);
00082 virtual ~StiKalmanTrackNode(){reduce();mId=-1;};
00083 const StiKalmanTrackNode& operator=(const StiKalmanTrackNode &node);
00084
00085 static double mcs2(double relRadThickness, double beta2, double p2) {return 14.1*14.1*relRadThickness/(beta2*p2*1e6);}
00087 void reset();
00088 void unset();
00090 void resetError(double fak=0);
00092 void initialize(StiHit*h);
00093 void initialize(StiDetector *d);
00094
00095
00097 void setState(const StiKalmanTrackNode * node);
00099 void get(double& alpha, double& xRef, double x[kNPars], double cc[kNErrs], double& chi2);
00100
00102 void getGlobalRadial(double x[6],double e[15]);
00103
00105 void getGlobalTpt (float x[6],float e[15]);
00106
00108 int getCharge() const {return (mFP.ptin() > 0) ? -1 : 1;}
00110 StThreeVectorF getMomentumF() const;
00113 StThreeVectorF getGlobalMomentumF() const;
00114 StThreeVector<double> getMomentum() const;
00115 StThreeVector<double> getGlobalMomentum() const;
00118 void getMomentum(double p[3], double e[6]=0) const;
00120 double getCurvature() const {return mFP.curv();}
00121 void setCurvature(double curvature) {mFP.curv()=curvature;}
00122 double getDipAngle() const {return atan(mFP.tanl());}
00123 double getTanL() const {return mFP.tanl();}
00125 double getPt() const;
00127 double getP() const {return (getPt()*::sqrt(1.+mFP.tanl()*mFP.tanl()));}
00130 double getHz() const;
00131 double getField() const {return getHz();}
00132 double x_g() const;
00133 double y_g() const;
00134 double z_g() const;
00135 void getXYZ_g(double *xyz) const;
00136 double getX() const { return mFP.x();}
00137 double getY() const { return mFP.y();}
00138 double getZ() const { return mFP.z();}
00139 double x() const { return mFP.x();}
00140 double y() const { return mFP.y();}
00141 double z() const { return mFP.z();}
00142 double getRxy() const { return sqrt(mFP.x()*mFP.x()+mFP.y()*mFP.y());}
00143
00144 double getEta () const {return mFP.eta(); }
00145 double getSin () const {return mFP._sinCA;}
00146 double getCos () const {return mFP._cosCA;}
00147 double getAlpha() const {return _alpha; }
00148 const double *hitErrs() const {return mHrr.A; }
00149 double getEyy() const {return mHrr.hYY;}
00150 double getEzz() const {return mHrr.hZZ;}
00151 double getCyy() const {return mFE._cYY;}
00152 double getCzz() const {return mFE._cZZ;}
00153 double const *getPars()const {return mFP.P;}
00154 double getDiag(int idx)const {return mFE.A[(idx*(idx+3))/2];}
00155 int getHitCount () const {return hitCount;}
00156 int getNullCount() const {return nullCount;}
00157 int getContigHitCount () const {return contiguousHitCount ;}
00158 int getContigNullCount() const {return contiguousNullCount;}
00159
00160
00161
00162
00163 int incHitCount () {return ++hitCount;}
00164 int incNullCount() {return ++nullCount;}
00165 int incContigHitCount () {return ++contiguousHitCount ;}
00166 int incContigNullCount() {return ++contiguousNullCount;}
00167 void setHitCount (char c=0) { hitCount=c;}
00168 void setNullCount(char c=0) { nullCount=c;}
00169 void setContigHitCount (char c=0) { contiguousHitCount=c ;}
00170 void setContigNullCount(char c=0) { contiguousNullCount=c;}
00171 double getTime() const;
00172
00173 void setHitCand(int nhits) {mHitCand = nhits;}
00174 void setIHitCand(int ihit) {mIHitCand = ihit;}
00175 int getHitCand()const {return mHitCand;}
00176 int getIHitCand()const {return mIHitCand;}
00177 static void Break(int kase);
00178 static void PrintStep();
00179 StThreeVector<double>getPoint() const;
00180 StThreeVector<double>getGlobalPoint() const;
00182 void getGlobalMomentum(double p[3], double e[6]=0) const;
00183 int isEnded() const;
00184 int isDca() const;
00185
00187 int propagate(StiKalmanTrackNode *p, const StiDetector * tDet, int dir);
00188
00190 bool propagate(const StiKalmanTrackNode *p, StiHit * vertex, int dir);
00191
00192 bool propagateToBeam(const StiKalmanTrackNode *p, int dir);
00193 int propagateToRadius(StiKalmanTrackNode *pNode, double radius,int dir);
00194
00201 double evaluateDedx();
00202
00203 int locate();
00204 int propagate(double x,int option,int dir);
00205 void propagateMtx();
00206 void propagateError();
00207 void saveInfo(int kase=1);
00208 const StiNodeInf *getInfo() const {return _inf;}
00209 int testError(double *emx,int begend);
00210 void numeDeriv(double val,int kind,int shape=0,int dir=0);
00211 int testDeriv(double *der);
00212 void propagateMCS(StiKalmanTrackNode * previousNode, const StiDetector * tDet);
00213
00216 StThreeVector<double> getPointAt(double xk) const;
00217
00218 int nudge(StiHit *hit=0);
00219 double evaluateChi2(const StiHit *hit);
00220 int updateNode();
00221 int rotate(double alpha);
00222 int getHelicity()const {return (mFP.curv() < 0) ? -1 : 1;}
00223 double getPhase() const;
00224 double getPsi() const;
00225 double getWindowY();
00226 double getWindowZ();
00227 double pitchAngle() const {return atan(mFP.tanl());}
00228 double crossAngle() const {return asin(mFP._sinCA);}
00229 double sinCrossAngle() const {return mFP._sinCA;}
00230 double pathlength() const;
00231 double pathLToNode(const StiKalmanTrackNode * const oNode);
00232 StThreeVectorD* getLengths(StiKalmanTrackNode *nextNode);
00233
00234 double length(const StThreeVector<double>& delta, double curv);
00235 double getDedx() const;
00236 static double nice(double angle);
00238 StThreeVector<double> getHelixCenter() const;
00239 void setHitErrors(const StiHit *hit=0);
00240 StiHitErrs getGlobalHitErrs(const StiHit *hit) const;
00241 friend ostream& operator<<(ostream& os, const StiKalmanTrackNode& n);
00242
00243 double getX0() const;
00244 double getGasX0() const;
00245 double getDensity() const;
00246 double getGasDensity() const;
00247
00248 void extend();
00249 void reduce();
00250 static Int_t debug() {return _debug;}
00251 static void setDebug(Int_t m) {_debug = m;}
00252 static void SetLaser(Int_t m) {_laser = m;}
00253 static Int_t IsLaser() {return _laser;}
00254 void PrintpT(const Char_t *opt="") const ;
00255 int getFlipFlop() const {return mFlipFlop;}
00256 static void ResetComment(const Char_t *m = "") {comment = m; commentdEdx = "";}
00257 static const Char_t *Comment() {return comment.Data();}
00259 int print(const char *opt) const;
00260
00261 private:
00262 void extinf();
00263 void static saveStatics(double *sav);
00264 void static backStatics(double *sav);
00265 static StiNodeExt *nodeExtInstance();
00266 static StiNodeInf *nodeInfInstance();
00267 void propagateCurv(const StiKalmanTrackNode *parent);
00268
00269
00270 public:
00271
00272 const StiNodePars &fitPars() const {return mFP; }
00273 StiNodePars &fitPars() {return mFP; }
00274 const StiNodeErrs &fitErrs() const {return mFE; }
00275 StiNodeErrs &fitErrs() {return mFE; }
00276 const StiNodePars &mPP() const {return _ext->mPP; }
00277 StiNodePars &mPP() {if (!_ext) extend(); return _ext->mPP; }
00278 const StiNodeErrs &mPE() const {return _ext->mPE; }
00279 StiNodeErrs &mPE() {if (!_ext) extend(); return _ext->mPE; }
00280 const StiNodeMtx &mMtx()const {return _ext->mMtx;}
00281 StiNodeMtx &mMtx() {if (!_ext) extend();return _ext->mMtx;}
00282 const StiNode2Pars &unTouched() const {return mUnTouch;}
00283 void setUntouched();
00284 protected:
00285
00286 char _beg[1];
00287 double _alpha;
00289 mutable double mHz;
00290 StiNodePars mFP;
00292 StiNodeErrs mFE;
00293 StiNode2Pars mUnTouch;
00294 StiHitErrs mHrr;
00295 char hitCount;
00296 char nullCount;
00297 char contiguousHitCount;
00298 char contiguousNullCount;
00299 char mFlipFlop;
00300 char mHitCand;
00301 char mIHitCand;
00302 char _end[1];
00303 StiNodeExt *_ext;
00304 StiNodeInf *_inf;
00305 static StiNodeStat mgP;
00306 static bool useCalculatedHitError;
00307
00308 static int fDerivTestOn;
00309 static double fDerivTest[kNPars][kNPars];
00310 static int _debug;
00311 static TString comment;
00312 static TString commentdEdx;
00313 static int _laser;
00314 public:
00315 int mId;
00316 };
00317
00318
00319 inline double StiKalmanTrackNode::nice(double angle)
00320 {
00321 if (angle <= -M_PI) angle += 2*M_PI;
00322 if (angle > M_PI) angle -= 2*M_PI;
00323 return angle;
00324 }
00325
00326 inline StThreeVector<double> StiKalmanTrackNode::getMomentum() const
00327 {
00328 double pt = getPt();
00329 return StThreeVector<double>(pt*mFP._cosCA,pt*mFP._sinCA,pt*mFP.tanl());
00330 }
00331
00332 inline StThreeVectorF StiKalmanTrackNode::getMomentumF() const
00333 {
00334 double pt = getPt();
00335 return StThreeVectorF(pt*mFP._cosCA,pt*mFP._sinCA,pt*mFP.tanl());
00336 }
00337
00338 inline StThreeVector<double> StiKalmanTrackNode::getGlobalMomentum() const
00339 {
00340 StThreeVector<double> p = getMomentum();
00341 p.rotateZ(_alpha);
00342 return p;
00343 }
00344
00345 inline StThreeVectorF StiKalmanTrackNode::getGlobalMomentumF() const
00346 {
00347 StThreeVectorF p = getMomentumF();
00348 p.rotateZ(_alpha);
00349 return p;
00350 }
00351
00352
00353
00354
00355 struct StiKTNXLessThan
00356 {
00357 bool operator()(const StiKalmanTrackNode& lhs, const StiKalmanTrackNode& rhs) const;
00358 };
00359
00360 struct StreamX
00361 {
00362 void operator()(const StiKalmanTrackNode& node)
00363 {
00364 cout <<node.getX()<<endl;
00365 }
00366 };
00367
00371 inline double StiKalmanTrackNode::pathlength() const
00372 {
00373 const StiDetector * det = getDetector();
00374 if (!det) return 0.;
00375 double thickness = det->getShape()->getThickness();
00376 return (thickness*::sqrt(1.+mFP.tanl()*mFP.tanl())) / mFP._cosCA;
00377 }
00378
00381 inline double StiKalmanTrackNode::getX0() const
00382 {
00383 const StiDetector * det = getDetector();
00384 if (!det)
00385 return 0.;
00386 return det->getMaterial()->getX0();
00387 }
00388
00391 inline double StiKalmanTrackNode::getGasX0() const
00392 {
00393 const StiDetector * det = getDetector();
00394 if (!det)
00395 return 0.;
00396 return det->getGas()->getX0();
00397 }
00398
00399 inline double StiKalmanTrackNode::getDensity() const
00400 {
00401 const StiDetector * det = getDetector();
00402 if (!det)
00403 return 0.;
00404 return det->getMaterial()->getDensity();
00405 }
00406
00407 inline double StiKalmanTrackNode::getGasDensity() const
00408 {
00409 const StiDetector * det = getDetector();
00410 if (!det)
00411 return 0.;
00412 return det->getGas()->getDensity();
00413 }
00414
00415
00416 inline StThreeVectorD* StiKalmanTrackNode::getLengths(StiKalmanTrackNode* nextNode)
00417 {
00418 double x1=pathlength()/2.;
00419 double x3=nextNode->pathlength()/2.;
00420 double x2=pathLToNode(nextNode);
00421 if (x2> (x1+x3)) x2=x2-x1-x3;
00422 else x2=0;
00423
00424 return new StThreeVectorD(x1/getX0(),
00425 x2/getDetector()->getMaterial()->getX0(),
00426 x3/nextNode->getX0());
00427 }
00428
00429 inline double StiKalmanTrackNode::getDedx() const
00430 {
00431
00432 StiHit *hit = getHit();
00433 if (!hit) return -1;
00434 double de=hit->getEloss();
00435 double dx=pathlength();
00436 if(dx>0 && de>0) return de/dx;
00437 return -1;
00438 }
00439
00440 #endif
00441
