StRoot: Sti/StiTrackNodeHelper.cxx Source File

00001 #include <stdio.h> 00002 #include <stdlib.h> 00003 #include "StiUtilities/StiDebug.h" 00004 #include "StiTrackNodeHelper.h" 00005 #include "StiElossCalculator.h" 00006 #include "StiHitErrorCalculator.h" 00007 #include "StMessMgr.h" 00008 #include "TArrayD.h" 00009 #if ROOT_VERSION_CODE < 331013 00010 #include "TCL.h" 00011 #else 00012 #include "TCernLib.h" 00013 #endif 00014 00015 #define NICE(a) ( ((a) <= -M_PI)? ((a)+2*M_PI) :\ 00016 ((a) > M_PI)? ((a)-2*M_PI) : (a)) 00017 00018 #define sinX(a) StiTrackNode::sinX(a) 00019 static const double kMaxEta = 1.5; 00020 static const double kMaxCur = 0.2; 00021 00022 static const double DY=0.9,DZ=0.9,DEta=0.1,DPti=3,DTan=0.1; 00023 static const double MAXSTEP[]={0,DY,DZ,DEta,DPti,DTan}; 00024 static const double ERROR_FACTOR = 2.; 00025 int StiTrackNodeHelper::_debug = 0; 00026 int StiTrackNodeHelper::mgCutStep=0; 00027 //______________________________________________________________________________ 00028 int errTest(StiNodePars &predP,StiNodeErrs &predE, 00029 const StiHit *hit,StiHitErrs &hitErr, 00030 StiNodePars &fitdP,StiNodeErrs &fitdE,double chi2); 00031 00032 //______________________________________________________________________________ 00033 void StiTrackNodeHelper::set(double chi2Max,double chi2Vtx,double errConfidence,int iter) 00034 { 00035 reset(); 00036 mChi2Max = chi2Max; 00037 mChi2Vtx = chi2Vtx; 00038 mNodeErrFactor = 10; 00039 mHitsErrFactor = 1.; 00040 if (errConfidence>0.1) { 00041 mNodeErrFactor = (1./(errConfidence)); 00042 mHitsErrFactor = (1./(1. - errConfidence)); 00043 } 00044 mIter = iter; 00045 if (!mIter) mFlipFlopNode = 0; 00046 } 00047 //______________________________________________________________________________ 00048 void StiTrackNodeHelper::reset() 00049 { 00050 memset(mBeg,'A',mEnd-mBeg+1); 00051 mWorstNode =0; 00052 mVertexNode=0; 00053 mUsed = 0; 00054 mCurvQa.reset(); 00055 mTanlQa.reset(); 00056 } 00057 00058 //______________________________________________________________________________ 00059 void StiTrackNodeHelper::set(StiKalmanTrackNode *pNode,StiKalmanTrackNode *sNode) 00060 { 00061 if(!pNode) reset(); 00062 mParentNode = pNode; 00063 mTargetNode = sNode; 00064 mTargetHz = mTargetNode->getHz(); 00065 mParentHz = mTargetHz; 00066 if (mParentNode) { 00067 mParentHz = mParentNode->getHz(); 00068 assert(fabs(mParentHz-mParentNode->mFP._hz)<1e-10); 00069 mParentNode->mFP.check("2StiTrackNodeHelper::set"); 00070 } 00071 if (mTargetNode->isValid()) { 00072 mTargetNode->mFP.check("1StiTrackNodeHelper::set"); 00073 assert(fabs(mTargetHz-mTargetNode->mFP._hz)<1e-10); 00074 } 00075 00076 mDetector = mTargetNode->getDetector(); 00077 if (!mDetector) mVertexNode = mTargetNode; 00078 mHit = mTargetNode->getHit(); 00079 if (mHit) { 00080 double time = 0; 00081 if (mHit->vy() || mHit->vz()) time = mTargetNode->getTime(); 00082 mHitPars[0]=mHit->x(); 00083 mHitPars[1]=mHit->y(time); 00084 mHitPars[2]=mHit->z(time); 00085 } 00086 00087 if (!mIter) mTargetNode->mFlipFlop=0; 00088 } 00089 //______________________________________________________________________________ 00090 int StiTrackNodeHelper::propagatePars(const StiNodePars &parPars 00091 , StiNodePars &rotPars 00092 , StiNodePars &proPars) 00093 { 00094 int ierr = 0; 00095 alpha = mTargetNode->_alpha - mParentNode->_alpha; 00096 ca=1;sa=0; 00097 parPars.check("1propagatePars"); 00098 rotPars = parPars; 00099 if (fabs(alpha) > 1.e-6) { //rotation part 00100 00101 double xt1=parPars._x; 00102 double yt1=parPars._y; 00103 double cosCA0 = parPars._cosCA; 00104 double sinCA0 = parPars._sinCA; 00105 00106 ca = cos(alpha); 00107 sa = sin(alpha); 00108 00109 rotPars._x = xt1*ca + yt1*sa; 00110 rotPars._y= -xt1*sa + yt1*ca; 00111 rotPars._cosCA = cosCA0*ca+sinCA0*sa; 00112 rotPars._sinCA = -cosCA0*sa+sinCA0*ca; 00113 double nor = 0.5*(rotPars._sinCA*rotPars._sinCA+rotPars._cosCA*rotPars._cosCA +1); 00114 rotPars._cosCA /= nor; 00115 rotPars._sinCA /= nor; 00116 rotPars._eta= NICE(parPars._eta-alpha); 00117 }// end of rotation part 00118 ierr = rotPars.check(); 00119 if (ierr) return 1; 00120 00121 // Propagation 00122 x1 = rotPars._x; 00123 x2 = (mDetector)? mDetector->getPlacement()->getNormalRadius():mHitPars[0]; 00124 dx = x2-x1; 00125 rho = 0.5*(mTargetHz*rotPars._ptin+rotPars._curv); 00126 dsin = rho*dx; 00127 sinCA2=rotPars._sinCA + dsin; 00128 if (sinCA2> 0.95) sinCA2= 0.95; 00129 if (sinCA2<-0.95) sinCA2=-0.95; 00130 cosCA2 = ::sqrt((1.-sinCA2)*(1.+sinCA2)); 00131 sumSin = rotPars._sinCA+sinCA2; 00132 sumCos = rotPars._cosCA+cosCA2; 00133 dy = dx*(sumSin/sumCos); 00134 y2 = rotPars._y+dy; 00135 dl0 = rotPars._cosCA*dx+rotPars._sinCA*dy; 00136 sind = dl0*rho; 00137 if (fabs(dsin) < 0.02 && rotPars._cosCA >0) { //tiny angle 00138 dl = dl0*(1.+sind*sind/6); 00139 } else { 00140 double cosd = cosCA2*rotPars._cosCA+sinCA2*rotPars._sinCA; 00141 dl = atan2(sind,cosd)/rho; 00142 } 00143 00144 proPars._x = x2; 00145 proPars._y = y2; 00146 proPars._z= rotPars._z + dl*rotPars._tanl; 00147 proPars._eta = (rotPars._eta+rho*dl); 00148 proPars._eta = NICE(proPars._eta); 00149 proPars._ptin = rotPars._ptin; 00150 proPars._hz = mTargetHz; 00151 proPars._curv = proPars._ptin*mTargetHz; 00152 proPars._tanl = rotPars._tanl; 00153 proPars._sinCA = sinCA2; 00154 proPars._cosCA = cosCA2; 00155 ierr = proPars.check(); 00156 if (ierr) return 2; 00157 return 0; 00158 } 00159 //______________________________________________________________________________ 00160 int StiTrackNodeHelper::propagateFitd() 00161 { 00162 // account eloss in matrix of transofrmation 00163 // remember that in matrix diagonal and correction factors 00164 // 1. subtructed. 00165 int ierr = 0; 00166 mBestPars._ptin *= (1+mMcs._ptinCorr); 00167 mBestPars._curv *= (1+mMcs._ptinCorr); 00168 mMtx.A[4][4] = (mMtx.A[4][4]+1)*(1+mMcs._ptinCorr) -1; 00169 00170 StiNodePars rotPars; 00171 ierr = propagatePars(mFitdParentPars,rotPars,mPredPars); 00172 if (ierr) return 1; 00173 // cutStep(&mBestPars,&mPredPars); 00174 mPredPars._hz = mTargetHz; 00175 mPredPars._ptin *= (1+mMcs._ptinCorr); 00176 mPredPars._curv *= (1+mMcs._ptinCorr); 00177 mPredPars.ready(); 00178 ierr = mPredPars.check(); 00179 if (ierr) return 2; 00180 return 0; 00181 } 00182 00183 00184 00185 //______________________________________________________________________________ 00186 int StiTrackNodeHelper::propagateMtx() 00187 { 00188 // fYE == dY/dEta 00189 double fYE= dx*(1.+mBestParentRotPars._cosCA*cosCA2+mBestParentRotPars._sinCA*sinCA2)/(sumCos*cosCA2); 00190 // fZE == dZ/dEta 00191 double dLdEta = dy/cosCA2; 00192 double fZE = mBestPars._tanl*dLdEta; 00193 // fZT == dZ/dTanL; 00194 double fZT= dl; 00195 00196 00197 // fEC == dEta/dRho 00198 double fEC = dx/cosCA2; 00199 // fYC == dY/dRho 00200 double fYC=(dl0)/sumCos*fEC; 00201 // fZC == dZ/dRho 00202 double dang = dl*rho; 00203 double C2LDX = dl*dl*( 00204 0.5*sinCA2*pow((1+pow(dang/2,2)*sinX(dang/2)),2) + 00205 cosCA2*dang*sinX(dang)); 00206 00207 double fZC = mBestPars._tanl*C2LDX/cosCA2; 00208 00209 fEC*=mTargetHz; fYC*=mTargetHz;fZC*=mTargetHz; 00210 00211 00212 mMtx.reset(); 00213 // X related derivatives 00214 mMtx.A[0][0] = -1; 00215 mMtx.A[1][0] = -sinCA2/cosCA2; 00216 mMtx.A[2][0] = -mBestPars._tanl/cosCA2 ; 00217 mMtx.A[3][0] = -mBestPars._curv/cosCA2 ; ; 00218 00219 mMtx.A[1][3]=fYE; mMtx.A[1][4]=fYC; mMtx.A[2][3]=fZE; 00220 mMtx.A[2][4]=fZC; mMtx.A[2][5]=fZT; mMtx.A[3][4]=fEC; 00221 double fYX = mMtx.A[1][0]; 00222 mMtx.A[1][0] = fYX*ca-sa; 00223 mMtx.A[1][1] = fYX*sa+ca-1; 00224 return 0; 00225 } 00226 //______________________________________________________________________________ 00227 int StiTrackNodeHelper::propagateError() 00228 { 00229 mPredErrs = mFitdParentErrs; 00230 StiTrackNode::errPropag6(mPredErrs.A,mMtx.A,kNPars); 00231 mPredErrs.recov(); 00232 mPredErrs._cEE+=mMcs._cEE; //add err to <eta*eta> eta crossing angle//add err to <eta*eta> eta crossing angle 00233 mPredErrs._cPP+=mMcs._cPP; //add err to <curv*curv> //add err to <curv*curv> 00234 mPredErrs._cTP+=mMcs._cTP; //add err to <tanL*curv> //add err to <tanL*curv> 00235 mPredErrs._cTT+=mMcs._cTT; //add err to <tanL*tanL> //add err to <tanL*tanL> 00236 int ierr = mPredErrs.check(); 00237 if (ierr) return 1; 00238 return 0; 00239 } 00240 00241 //______________________________________________________________________________ 00242 int StiTrackNodeHelper::makeFit(int smooth) 00243 { 00244 static int nCall=0; 00245 nCall++; 00246 StiDebug::Break(nCall); 00247 int ierr=0; 00248 mState = 0; 00249 mChi2 = 1e13; 00250 if (mParentNode) { 00251 // Select the best params to make matrix 00252 00253 // Get best params for derivatives 00254 if (smooth) {//joined pars always the best 00255 mBestParentPars = mJoinPars; 00256 mBestParentErrs = mJoinErrs; 00257 mBestDelta = mJoinErrs.getDelta(); 00258 } else { 00259 double delta = mFitdErrs.getDelta(); 00260 double er1 = delta*delta; 00261 double er2 = mSavdDelta*mSavdDelta; 00262 double wt = er1/(er1+er2); 00263 mBestParentPars = mFitdPars; 00264 mBestParentPars.merge(wt,mSavdParentPars); 00265 mBestDelta = sqrt(er1*er2/(er1+er2)); 00266 00267 } 00268 mBestParentPars.check("1makeFit"); 00269 mFitdParentPars = mFitdPars; 00270 mFitdParentErrs = mFitdErrs; 00271 mFitdParentPars.check("2makeFit"); 00272 mFitdParentErrs.check("3makeFit"); 00273 00274 ierr = propagatePars(mBestParentPars,mBestParentRotPars,mBestPars); 00275 if(ierr) return 1; 00276 ierr = propagateMtx(); if(ierr) return 2; 00277 ierr = propagateMCS(); if(ierr) return 3; 00278 ierr = propagateFitd(); if(ierr) return 4; 00279 ierr = propagateError(); if(ierr) return 5; 00280 00281 } 00282 mState = StiTrackNode::kTNReady; 00283 00284 // Save parameters for future, 00285 // when target node will became parent one 00286 mSavdParentPars = mTargetNode->mFP; 00287 mSavdDelta = (mTargetNode->isValid())? mTargetNode->mFE.getDelta():3e33; 00288 00289 if (!mParentNode) { 00290 if (!smooth) mgCutStep = 0; 00291 mPredErrs = mTargetNode->mFE; 00292 ierr = mPredErrs.check(); if (ierr) return 11; 00293 mPredPars = mTargetNode->mFP; 00294 ierr = mPredPars.check(); if (ierr) return 12; 00295 mBestPars = mPredPars; 00296 mBestDelta = mPredErrs.getDelta(); 00297 mJoinPars = mPredPars; 00298 resetError(mNodeErrFactor); 00299 } 00300 // Set fitted pars to predicted for the absent hit case 00301 mFitdPars = mPredPars; 00302 mFitdErrs = mPredErrs; 00303 00304 int ians = 1; 00305 mChi2 =0; 00306 do {//technical (fake) loop 00307 if (!mHit) break; 00308 setHitErrs(); 00309 if (nudge()) return 13; 00310 mChi2 = 3e33; 00311 double chi2 = evalChi2(); 00312 if (mTargetNode == mVertexNode) { 00313 if (chi2>mChi2Vtx) return 14; 00314 } else { 00315 if (chi2>mChi2Max) break; 00316 } 00317 mChi2 = chi2; if (mChi2>999) mChi2=999; 00318 ians = updateNode(); 00319 if (!ians) break; 00320 if (debug() & 8) {cout << Form("%5d ",ians); StiKalmanTrackNode::PrintStep();} 00321 if (mTargetNode == mVertexNode) return 15; 00322 mState = StiTrackNode::kTNReady; 00323 mFitdPars = mPredPars; 00324 mFitdErrs = mPredErrs; 00325 mChi2 = 3e33; 00326 }while(0); 00327 00328 ierr = (!smooth)? save():join(); 00329 if (ierr) return 16; 00330 00331 00332 do { //fake loop 00333 if (!smooth) break; 00334 if (!mHit) break; 00335 if (mDetector && (!mFlipFlopNode || mTargetNode->mFlipFlop > mFlipFlopNode->mFlipFlop)) 00336 {mFlipFlopNode=mTargetNode;} 00337 if(mState!=StiTrackNode::kTNFitEnd) break; 00338 mUsed++; 00339 if (mDetector && (!mWorstNode || mChi2>mWorstNode->getChi2())) 00340 {mWorstNode=mTargetNode;} 00341 if (!mParentNode) break; 00342 double accu,errr; 00343 accu = mJoinPars._ptin - mBestParentPars._ptin*(1+mMcs._ptinCorr); 00344 errr = sqrt(0.5*(fabs(mJoinErrs._cPP)+fabs(mBestParentErrs._cPP))); 00345 mCurvQa.add(accu/errr); 00346 accu = mJoinPars._tanl - mBestParentPars._tanl; 00347 errr = sqrt(0.5*(fabs(mJoinErrs._cTT)+fabs(mBestParentErrs._cTT))); 00348 mTanlQa.add(accu/errr); 00349 }while(0); 00350 00351 return 0; 00352 } 00353 00354 //______________________________________________________________________________ 00355 int StiTrackNodeHelper::join() 00356 { 00357 enum {kOLdValid=1,kNewFitd=2,kJoiUnFit=4}; 00358 // if (!mParentNode) return 0; 00359 00360 00361 int ierr = 0; 00362 double chi2; 00363 00364 StiDebug::Break(mTargetNode->mId); 00365 int kase = mTargetNode->isValid(); 00366 if (mState==StiTrackNode::kTNFitEnd) kase |=kNewFitd; 00367 static int oldJoinPrim = StiDebug::iFlag("StiOldJoinPrim"); 00368 if (!oldJoinPrim) { 00369 if (mTargetNode==mVertexNode) kase = kNewFitd; //ignore old info for primVtx 00370 //Hack to accoont specific 00371 //fit to primVtx 00372 } 00373 do { 00374 switch(kase) { 00375 case 0: // Old invalid & New UnFitd 00376 mChi2 = (mHit)? 3e33:0; 00377 mState = StiTrackNode::kTNReady; 00378 00379 case kNewFitd: // Old invalid & New Fitd 00380 mJoinPars = mFitdPars; 00381 mJoinErrs = mFitdErrs; 00382 kase = -1; break; 00383 00384 case kOLdValid|kNewFitd|kJoiUnFit: // Old valid & New Fitd & Join UnFit 00385 mChi2 = 3e33; 00386 mFitdPars = mPredPars; 00387 mFitdErrs = mPredErrs; 00388 mState = StiTrackNode::kTNReady; 00389 00390 case kOLdValid:; // Old valid & New UnFitd 00391 case kOLdValid|kNewFitd:; // Old valid & New Fitd 00392 00393 chi2 = joinTwo(kNPars,mTargetNode->mPP().P,mTargetNode->mPE().A 00394 ,kNPars,mFitdPars.P ,mFitdErrs.A 00395 ,mJoinPars.P ,mJoinErrs.A); 00396 mJoinPars._hz = mTargetHz; 00397 mJoinErrs.recov(); 00398 00399 00400 if (kase == (kOLdValid|kNewFitd)) { //Check errors improvements 00401 assert(mHrr.hYY > mJoinErrs._cYY); 00402 assert(mHrr.hZZ > mJoinErrs._cZZ); 00403 if (mHrr.hYY <= mJoinErrs._cYY) { 00404 LOG_DEBUG << Form("StiTrackNodeHelper::updateNode() WRONG hYY(%g) < nYY(%g)" 00405 ,mHrr.hYY,mFitdErrs._cYY)<< endm; 00406 return -13; 00407 } 00408 if (mHrr.hZZ <= mJoinErrs._cZZ) { 00409 LOG_DEBUG << Form("StiTrackNodeHelper::updateNode() WRONG hZZ(%g) < nZZ(%g)" 00410 ,mHrr.hZZ,mFitdErrs._cZZ) << endm; 00411 return -14; 00412 } 00413 } //End Check errors improvements 00414 mJoinPars.ready(); 00415 ierr = mJoinPars.check(); if (ierr) return 2; 00416 if (kase!=(kOLdValid|kNewFitd)) {kase = -1; break;} 00417 00418 mChi2 = 3e33; 00419 // chi2 = joinChi2(); //Test join Chi2 00420 chi2 = recvChi2(); //Test join Chi2 00421 mChi2 = 3e33; 00422 if (chi2>mChi2Max && mTargetNode!=mVertexNode) 00423 { kase |=kJoiUnFit; break;} //join Chi2 too big 00424 mChi2 = (chi2>999)? 999:chi2; 00425 mState = StiTrackNode::kTNFitEnd; 00426 kase = -1; break; 00427 00428 default: assert(0); 00429 }//end Switch 00430 } while(kase>=0); 00431 00432 assert(fabs(mJoinPars._hz-mTargetHz)<=1e-10); 00433 assert(fabs(mTargetNode->getHz()-mTargetHz)<=1e-10); 00434 00435 00436 mTargetNode->mFE = mJoinErrs; 00437 mTargetNode->mPE() = mJoinErrs; 00438 mTargetNode->mFP = mJoinPars; 00439 mTargetNode->mPP() = mJoinPars; 00440 mTargetNode->mHrr = mHrr; 00441 if (mTargetNode!=mVertexNode) mTargetNode->mUnTouch = mUnTouch; 00442 mTargetNode->_state = mState; 00443 if (mHit && ((mChi2<1000) != (mTargetNode->getChi2()<1000))) mTargetNode->mFlipFlop++; 00444 if (mTargetNode!=mVertexNode) mTargetNode->setChi2(mChi2); 00445 00446 // Sanity check 00447 double myMaxChi2 = (mTargetNode==mVertexNode)? 1000:mChi2Max; 00448 if (mState == StiTrackNode::kTNFitEnd) { 00449 assert(mChi2 <myMaxChi2); 00450 } else { 00451 assert(mChi2 <=0 || mChi2 >=myMaxChi2); 00452 } 00453 00454 return 0; 00455 } 00456 //______________________________________________________________________________ 00457 double StiTrackNodeHelper::joinTwo(int nP1,const double *P1,const double *E1 00458 ,int nP2,const double *P2,const double *E2 00459 , double *PJ, double *EJ) 00460 { 00461 00462 assert(nP1<=nP2); 00463 int nE1 = nP1*(nP1+1)/2; 00464 int nE2 = nP2*(nP2+1)/2; 00465 TArrayD ard(nE2*6); 00466 double *a = ard.GetArray(); 00467 double *sumE = (a); 00468 double *sumEI = (a+=nE2); 00469 double *e1sumEIe1 = (a+=nE2); 00470 double *subP = (a+=nE2); 00471 double *sumEIsubP = (a+=nE2); 00472 double chi2=3e33,p,q; 00473 00474 // Choose the smalest errors 00475 const double *p1 = P1, *p2 = P2, *e1 = E1, *e2 = E2, *t; 00476 double choice = (nP1==nP2)? 0:1; 00477 if (!choice ) { 00478 for (int i=0,n=1;i<nE2;i+=(++n)) { 00479 p=fabs(e1[i]);q=fabs(e2[i]);choice += (p-q)/(p+q+1e-10); 00480 }} 00481 if ( choice >0) {t = p2; p2 = p1; p1 = t; t = e2; e2 = e1; e1 = t;} 00482 00483 do {//empty loop 00484 // Join errors 00485 TCL::vadd(e1,e2,sumE,nE1); 00486 int negati = sumE[2]<0; 00487 if (negati) TCL::vcopyn(sumE,sumE,nE1); 00488 int ign0re = sumE[0]<=0; 00489 if (ign0re) sumE[0] = 1; 00490 TCL::trsinv(sumE,sumEI,nP1); 00491 if (ign0re) {sumE[0] = 0; sumEI[0] = 0;} 00492 if (negati) {TCL::vcopyn(sumE,sumE,nE1);TCL::vcopyn(sumEI,sumEI,nE1);} 00493 TCL::vsub(p2 ,p1 ,subP ,nP1); 00494 TCL::trasat(subP,sumEI,&chi2,1,nP1); 00495 if (!EJ) break; 00496 TCL::trqsq (e1 ,sumEI,e1sumEIe1,nP2); 00497 TCL::vsub(e1,e1sumEIe1,EJ,nE2); 00498 } while(0); 00499 // Join params 00500 if (PJ) { 00501 TCL::tras(subP ,sumEI,sumEIsubP,1,nP1); 00502 TCL::tras(sumEIsubP,e1 ,PJ ,1,nP2); 00503 TCL::vadd(PJ ,p1 ,PJ ,nP2); 00504 } 00505 return chi2; 00506 } 00507 #if 0 00508 //______________________________________________________________________________ 00509 double StiTrackNodeHelper::joinVtx(const double *P1,const double *E1 00510 ,const double *P2,const double *E2 00511 , double *PJ, double *EJ) 00512 { 00513 enum {nP1=3,nE1=6,nP2=kNPars,nE2=kNErrs}; 00514 double E1m[nE1],E2m[nE2]; 00515 00516 TCL::vzero(E1m ,nE1); 00517 TCL::ucopy(E2,E2m ,nE2); 00518 TCL::vadd (E1,E2m,E2m,nE1); 00519 return joinTwo(nP1,P1,E1m,nP2,P2,E2m,PJ,EJ); 00520 00521 } 00522 #endif//0 00523 #if 0 00524 //______________________________________________________________________________ 00525 double StiTrackNodeHelper::joinVtx(const double *P1,const double *E1 00526 ,const double *P2,const double *E2 00527 , double *PJ, double *EJ) 00528 { 00529 enum {nP1=3,nE1=6,nP2=kNPars,nE2=kNErrs}; 00530 00531 TArrayD ard(nE2*7); 00532 double *p = ard.GetArray(); 00533 double *sumE = (p); 00534 double *sumEI = (p+=nE2); 00535 double *sumEI_E1_sumEI = (p+=nE2); 00536 double *E2_sumEI_E2 = (p+=nE2); 00537 double *E2_sumEI_E1_sumEI_E2 = (p+=nE2); 00538 double *sumEIsubP = (p+=nE2); 00539 double *subP = (p+=nE2); 00540 00541 double chi2=3e33; 00542 00543 00544 do {//empty loop 00545 // Join errors 00546 TCL::ucopy(E2,sumE,nE1); 00547 int ign0re = sumE[0]<=0; 00548 if (ign0re) sumE[0] = 1; 00549 TCL::trsinv(sumE,sumEI,nP1); 00550 if (ign0re) {sumE[0] = 0; sumEI[0] = 0;} 00551 TCL::vsub(P1 ,P2 ,subP ,nP1); 00552 TCL::trasat(subP,sumEI,&chi2,1,nP1); 00553 if (!EJ) break; 00554 TCL::trqsq (E2 ,sumEI,E2_sumEI_E2,nP2); 00555 TCL::vsub(E2,E2_sumEI_E2,EJ,nE2); 00556 // Now account errors of vertex itself 00557 TCL::trqsq (sumEI,E1,sumEI_E1_sumEI,nP1); 00558 TCL::trqsq (E2,sumEI_E1_sumEI,E2_sumEI_E1_sumEI_E2,nP2); 00559 TCL::vadd(EJ,E2_sumEI_E1_sumEI_E2,EJ,nE2); 00560 00561 } while(0); 00562 // Join params 00563 if (PJ) { 00564 TCL::tras(subP ,sumEI,sumEIsubP,1,nP1); 00565 TCL::tras(sumEIsubP,E2 ,PJ ,1,nP2); 00566 TCL::vadd(PJ ,P2 ,PJ ,nP2); 00567 } 00568 return chi2; 00569 } 00570 #endif//1 00571 //______________________________________________________________________________ 00572 double StiTrackNodeHelper::joinVtx(const double *Y,const StiHitErrs &B 00573 ,const StiNodePars &X,const StiNodeErrs &A 00574 , StiNodePars *M, StiNodeErrs *C) 00575 { 00576 // Y : x,y,z of vertex. B: error matrix of Y 00577 // X : track parameters. A: error matrix of X 00578 // M : track parameters of X+Y. C: error matrix of M 00579 00580 00581 enum {nP1=3,nE1=6,nP2=6,nE2=21}; 00582 00583 StiNodeErrs Ai=A; //Inverted A 00584 00585 Ai._cXX=1; 00586 TCL::trsinv(Ai.A,Ai.A,nP2); 00587 Ai._cXX=0; 00588 00589 00590 double Ai11i[6],Ai10[3][3],T10[3][3],dif[6],m[6]; 00591 Ai.get11(Ai11i); 00592 Ai.get10(Ai10[0]); 00593 TCL::trsinv(Ai11i ,Ai11i,3); 00594 TCL::trsa (Ai11i,Ai10[0],T10[0],3,3); //Ai11*Ai10 00595 TCL::ucopy(Y,m,3); 00596 TCL::vsub (X.P,Y,dif,3); 00597 TCL::mxmpy(T10[0],dif,m+3,3,3,1); 00598 TCL::vadd(X.P+3,m+3,m+3,3); //m: resulting params 00599 if (M) {TCL::ucopy(m,M->P,nP2); M->ready();} //fill resulting params 00600 00601 TCL::vsub(X.P,m,dif,nP2); //dif = X - M 00602 double chi2; 00603 TCL::trasat(dif,Ai.A,&chi2,1,nP2); //calc chi2 00604 if (!C) return chi2; 00605 // Error matrix calculation 00606 C->reset(); 00607 00608 double TX[nP1][nP2];memset(TX[0],0,sizeof(TX)); 00609 for (int i=0;i<3;i++) {TCL::ucopy(T10[i],TX[i],3); TX[i][i+3]=1;} 00610 double C11[nE1]; 00611 TCL::trasat(TX[0],A.A,C11,nP1,nP2); 00612 C->set11(C11); 00613 00614 double TY[nP2][nP1];memset(TY[0],0,sizeof(TX)); 00615 for (int i=0;i<3;i++) {TCL::vcopyn(T10[i],TY[i+3],3); TY[i][i]=1;} 00616 TY[0][0] = 0; 00617 double CYY[nE2]; 00618 TCL::trasat(TY[0],B.A,CYY,nP2,nP1); 00619 TCL::vadd(CYY,C->A,C->A,nE2); 00620 return chi2; 00621 } 00622 //______________________________________________________________________________ 00623 int StiTrackNodeHelper::save() 00624 { 00625 assert(fabs(mPredPars._hz-mTargetHz)<=1e-10); 00626 assert(fabs(mFitdPars._hz-mTargetHz)<=1e-10); 00627 assert(fabs(mTargetNode->getHz()-mTargetHz)<=1e-10); 00628 00629 mTargetNode->mPP() = mPredPars; 00630 mTargetNode->mFP = mFitdPars; 00631 mTargetNode->mPE() = mPredErrs; 00632 mTargetNode->mFE = mFitdErrs; 00633 mTargetNode->_state = mState; 00634 if (mHit && ((mChi2<1000) != (mTargetNode->getChi2()<1000))) mTargetNode->mFlipFlop++; 00635 if (mTargetNode!=mVertexNode) mTargetNode->setChi2(mChi2); 00636 return 0; 00637 } 00638 //______________________________________________________________________________ 00647 //______________________________________________________________________________ 00648 int StiTrackNodeHelper::propagateMCS() 00649 { 00650 mMcs.reset(); 00651 if (!mDetector) return 0; 00652 if (fabs(mBestPars._ptin)<=1e-3) return 0; 00653 double pt = 1./fabs(mBestPars._ptin); 00654 00655 double relRadThickness; 00656 // Half path length in previous node 00657 double pL1,pL2,pL3,d1,d2,d3,dxEloss,dx; 00658 pL1=0.5*pathIn(mParentNode->getDetector(),&mBestParentPars); 00659 // Half path length in this node 00660 pL3=0.5*pathIn(mDetector,&mBestPars); 00661 // Gap path length 00662 pL2= fabs(dl); 00663 double x0p =-1; 00664 double x0Gas=-1; 00665 double x0=-1; 00666 dx = mBestPars._x - mBestParentRotPars._x; 00667 d1 = mParentNode->getDensity(); 00668 x0p = mParentNode->getX0(); 00669 d3 = mDetector->getMaterial()->getDensity(); 00670 x0 = mDetector->getMaterial()->getX0(); 00671 00672 00673 if (pL2> (pL1+pL3)) { 00674 00675 pL2=pL2-pL1-pL3; 00676 if (dx>0) { 00677 x0Gas = mDetector->getGas()->getX0(); 00678 d2 = mDetector->getGas()->getDensity(); 00679 } else { 00680 x0Gas = mParentNode->getGasX0(); 00681 d2 = mParentNode->getGasDensity(); 00682 } 00683 00684 relRadThickness = 0.; 00685 dxEloss = 0; 00686 00687 if (x0p>0.) { 00688 relRadThickness += pL1/x0p; 00689 dxEloss += d1*pL1; 00690 } 00691 if (x0Gas>0.) { 00692 relRadThickness += pL2/x0Gas; 00693 dxEloss += d2*pL2; 00694 } 00695 if (x0>0.){ 00696 relRadThickness += pL3/x0; 00697 dxEloss += d3*pL3; 00698 } 00699 } else { // pL2<=(pL1+pL3) 00700 relRadThickness = 0.; 00701 dxEloss = 0; 00702 if (x0p>0.) { 00703 relRadThickness += pL1/x0p; 00704 dxEloss += d1*pL1; 00705 } 00706 if (x0>0.){ 00707 relRadThickness += pL3/x0; 00708 dxEloss += d3*pL3; 00709 } 00710 } 00711 double tanl = mBestPars._tanl; 00712 double pti = mBestPars._ptin; 00713 double p2 = (1.+tanl*tanl)*pt*pt; 00714 double m = StiKalmanTrackFinderParameters::instance()->getMassHypothesis(); 00715 double m2 = m*m; 00716 double e2 = p2+m2; 00717 double beta2 = p2/e2; 00718 double theta2 = StiKalmanTrackNode::mcs2(relRadThickness,beta2,p2); 00719 double cos2Li = (1.+ tanl*tanl); // 1/cos(lamda)**2 00720 double f = mHitsErrFactor; 00721 mMcs._cEE = cos2Li *theta2*f; 00722 mMcs._cPP = tanl*tanl*pti*pti *theta2*f; 00723 mMcs._cTP = pti*tanl*cos2Li *theta2*f; 00724 mMcs._cTT = cos2Li*cos2Li *theta2*f; 00725 00726 double dE=0; 00727 double sign = (dx>0)? 1:-1; 00728 00729 // const static double I2Ar = (15.8*18) * (15.8*18) * 1e-18; // GeV**2 00730 StiElossCalculator * calculator = mDetector->getElossCalculator(); 00731 double eloss = calculator->calculate(1.,m, beta2); 00732 dE = sign*dxEloss*eloss; 00733 00734 mMcs._ptinCorr = ::sqrt(e2)*dE/p2; 00735 if (fabs(mMcs._ptinCorr)>0.1) mMcs._ptinCorr = (dE<0)? -0.1:0.1; 00736 return 0; 00737 } 00738 //______________________________________________________________________________ 00752 //______________________________________________________________________________ 00753 double StiTrackNodeHelper::evalChi2() 00754 { 00755 double r00, r01,r11,chi2; 00756 //If required, recalculate the errors of the detector hits. 00757 //Do not attempt this calculation for the main vertex. 00758 if (fabs(mPredPars._sinCA)>0.99 ) return 1e41; 00759 if (fabs(mPredPars._eta) >kMaxEta) return 1e41; 00760 if (fabs(mPredPars._curv) >kMaxCur) return 1e41; 00761 if (!mDetector) { //Primay vertex 00762 mHitPars[0] = mPredPars._x; 00763 // chi2 = joinVtx(mHitPars,mHrr.A,mPredPars.P,mPredErrs.A); 00764 chi2 = joinVtx(mHitPars,mHrr,mPredPars,mPredErrs); 00765 } else { //Normal hit 00766 00767 r00=mPredErrs._cYY+mHrr.hYY; 00768 r01=mPredErrs._cZY+mHrr.hZY; 00769 r11=mPredErrs._cZZ+mHrr.hZZ; 00770 mDetm = r00*r11 - r01*r01; 00771 if (mDetm<r00*r11*1.e-5) { 00772 LOG_DEBUG << Form("StiTrackNodeHelper::evalChi2 *** zero determinant %g",mDetm)<< endm; 00773 return 1e60; 00774 } 00775 double tmp=r00; r00=r11; r11=tmp; r01=-r01; 00776 00777 double dyt=(mPredPars._y-mHitPars[1]); 00778 double dzt=(mPredPars._z-mHitPars[2]); 00779 chi2 = (dyt*r00*dyt + 2*r01*dyt*dzt + dzt*r11*dzt)/mDetm; 00780 } 00781 return chi2; 00782 } 00783 //______________________________________________________________________________ 00784 /* 00785 * double StiTrackNodeHelper::joinChi2() 00786 * { 00787 * double chi2; 00788 * double mergPars[3],mHitPars[3],mergErrs[6]; 00789 * chi2 = joinTwo(3,mPredPars.P ,mPredErrs.A 00790 * ,3,mTargetNode->mPP().P,mTargetNode->mPE().A 00791 * ,mergPars ,mergErrs); 00792 * 00793 * mHitPars[0] = mHit->x(); 00794 * mHitPars[1] = mHit->y(); 00795 * mHitPars[2] = mHit->z(); 00796 * chi2 = joinTwo(3,mergPars,mergErrs,3,mHitPars,mHrr.A); 00797 * // Save untouched by current hit node's y,z & errors for alignment 00798 * mUnTouch.mPar[0] = mergPars[1]; 00799 * mUnTouch.mPar[1] = mergPars[2]; 00800 * mUnTouch.mErr[0] = mergErrs[2]; 00801 * mUnTouch.mErr[1] = mergErrs[4]; 00802 * mUnTouch.mErr[2] = mergErrs[5]; 00803 * return chi2; 00804 * } 00805 */ 00806 //______________________________________________________________________________ 00807 double StiTrackNodeHelper::recvChi2() 00808 { 00809 if (fabs(mJoinPars._sinCA)>0.99 ) return 1e41; 00810 if (fabs(mJoinPars._eta) >kMaxEta) return 1e41; 00811 if (fabs(mJoinPars._curv) >kMaxCur) return 1e41; 00812 if (!mDetector) {//Primary vertex 00813 // double chi2 =joinVtx(mHitPars,mHrr.A,mPredPars.P,mPredErrs.A); 00814 double chi2 =joinVtx(mHitPars,mHrr ,mPredPars ,mPredErrs ); 00815 return chi2; 00816 } 00817 00818 StiHitErrs myHrr = mHrr; 00819 StiNodeErrs recovErrs; 00820 StiNodePars recovPars; 00821 double f = -(1./mHitsErrFactor); 00822 myHrr*=f; 00823 double r11,r12,r22; 00824 if ((r11=myHrr.hYY+mJoinErrs._cYY) >=0) return 1e41; 00825 if ((r22=myHrr.hZZ+mJoinErrs._cZZ) >=0) return 1e41; 00826 r12 =myHrr.hZY+mJoinErrs._cZY; 00827 if (r11*r22-r12*r12<0) return 1e41; 00828 00829 00830 double chi2 = joinTwo(3,mHitPars , myHrr.A 00831 ,3,mJoinPars.P,mJoinErrs.A 00832 ,recovPars.P,recovErrs.A); 00833 if (fabs(recovPars._y-mHitPars[1])>10) StiDebug::Break(-1); 00834 if (fabs(recovPars._z-mHitPars[2])>10) StiDebug::Break(-1); 00835 00836 mUnTouch.set(recovPars,recovErrs); 00837 return -chi2; //account that result is negative 00838 } 00839 //______________________________________________________________________________ 00840 int StiTrackNodeHelper::setHitErrs() 00841 { 00842 getHitErrors(mHit,&mFitdPars,&mHrr); 00843 mHrr*=mHitsErrFactor; 00844 return 0; 00845 } 00846 //______________________________________________________________________________ 00866 //______________________________________________________________________________ 00867 int StiTrackNodeHelper::updateNode() 00868 { 00869 static int nCall=0; nCall++; 00870 mState = StiTrackNode::kTNFitBeg; 00871 assert(mPredErrs._cXX<1e-8); 00872 double r00,r01,r11; 00873 StiDebug::Break(mTargetNode->mId); 00874 if (!mDetector) { //Primary vertex 00875 mHitPars[0] = mPredPars._x; 00876 // double chi2 = joinVtx(mHitPars,mHrr.A,mPredPars.P,mPredErrs.A,mFitdPars.P,mFitdErrs.A); 00877 double chi2 = joinVtx(mHitPars,mHrr,mPredPars,mPredErrs,&mFitdPars,&mFitdErrs); 00878 mFitdPars._curv = mTargetHz*mFitdPars._ptin; 00879 assert(chi2>900 || fabs(mChi2-chi2)<1e-10); 00880 } else { //Normal Hit 00881 00882 r00=mHrr.hYY+mPredErrs._cYY; 00883 r01=mHrr.hZY+mPredErrs._cZY; 00884 r11=mHrr.hZZ+mPredErrs._cZZ; 00885 mDetm=r00*r11 - r01*r01; 00886 if (!(mDetm>(r00*r11)*1.e-5)) return 99; 00887 assert(mDetm>(r00*r11)*1.e-5); 00888 00889 // inverse matrix 00890 double tmp=r00; r00=r11/mDetm; r11=tmp/mDetm; r01=-r01/mDetm; 00891 // update error matrix 00892 double k00=mPredErrs._cYY*r00+mPredErrs._cZY*r01, k01=mPredErrs._cYY*r01+mPredErrs._cZY*r11; 00893 double k10=mPredErrs._cZY*r00+mPredErrs._cZZ*r01, k11=mPredErrs._cZY*r01+mPredErrs._cZZ*r11; 00894 double k20=mPredErrs._cEY*r00+mPredErrs._cEZ*r01, k21=mPredErrs._cEY*r01+mPredErrs._cEZ*r11; 00895 double k30=mPredErrs._cPY*r00+mPredErrs._cPZ*r01, k31=mPredErrs._cPY*r01+mPredErrs._cPZ*r11; 00896 double k40=mPredErrs._cTY*r00+mPredErrs._cTZ*r01, k41=mPredErrs._cTY*r01+mPredErrs._cTZ*r11; 00897 00898 double myY = mPredPars._y; 00899 double myZ = mPredPars._z; 00900 double dyt = mHitPars[1] - myY; 00901 double dzt = mHitPars[2] - myZ; 00902 double dPt = k30*dyt + k31*dzt; 00903 double dEt = k20*dyt + k21*dzt; 00904 double dTa = k40*dyt + k41*dzt; 00905 double eta = NICE(mPredPars._eta + dEt); 00906 double pti = mPredPars._ptin+ dPt; 00907 double tanl = mPredPars._tanl + dTa; 00908 // Check if any of the quantities required to pursue the update 00909 // are infinite. If so, it means the tracks cannot be update/propagated 00910 // any longer and should therefore be abandoned. Just return. This is 00911 // not a big but rather a feature of the fact a helicoidal tracks!!! 00912 // update Kalman state 00913 double p0 = myY + k00*dyt + k01*dzt; 00914 double p1 = myZ + k10*dyt + k11*dzt; 00915 //mPredPars._tanl += k40*dyt + k41*dzt; 00916 double sinCA = sin(eta); 00917 // The following test introduces a track propagation error but happens 00918 // only when the track should be aborted so we don't care... 00919 00920 mFitdPars._hz = mTargetHz; 00921 mFitdPars._x = mPredPars._x; 00922 mFitdPars._y = p0; 00923 mFitdPars._z = p1; 00924 mFitdPars._eta = eta; 00925 mFitdPars._ptin = pti; 00926 mFitdPars._curv = mTargetHz*pti; 00927 mFitdPars._tanl = tanl; 00928 mFitdPars._sinCA = sinCA; 00929 mFitdPars._cosCA = ::sqrt((1.-mFitdPars._sinCA)*(1.+mFitdPars._sinCA)); 00930 if (!mDetector) 00931 assert(fabs(mFitdPars._y-mHitPars[1])>1e-10 || fabs(mHitPars[0])<4); 00932 //?? cutStep(&mFitdPars,&mPredPars); 00933 //?? cutStep(&mFitdPars,&mBestPars); 00934 if (mFitdPars.check()) return -11; 00935 // update error matrix 00936 double c00=mPredErrs._cYY; 00937 double c10=mPredErrs._cZY, c11=mPredErrs._cZZ; 00938 double c20=mPredErrs._cEY, c21=mPredErrs._cEZ;//, c22=mPredErrs._cEE; 00939 double c30=mPredErrs._cPY, c31=mPredErrs._cPZ;//, c32=mPredErrs._cPE, c33=mPredErrs._cPP; 00940 double c40=mPredErrs._cTY, c41=mPredErrs._cTZ;//, c42=mPredErrs._cTE, c43=mPredErrs._cTP, c44=mPredErrs._cTT; 00941 mFitdErrs._cYY-=k00*c00+k01*c10; 00942 mFitdErrs._cZY-=k10*c00+k11*c10;mFitdErrs._cZZ-=k10*c10+k11*c11; 00943 mFitdErrs._cEY-=k20*c00+k21*c10;mFitdErrs._cEZ-=k20*c10+k21*c11;mFitdErrs._cEE-=k20*c20+k21*c21; 00944 mFitdErrs._cPY-=k30*c00+k31*c10;mFitdErrs._cPZ-=k30*c10+k31*c11;mFitdErrs._cPE-=k30*c20+k31*c21;mFitdErrs._cPP-=k30*c30+k31*c31; 00945 mFitdErrs._cTY-=k40*c00+k41*c10;mFitdErrs._cTZ-=k40*c10+k41*c11;mFitdErrs._cTE-=k40*c20+k41*c21;mFitdErrs._cTP-=k40*c30+k41*c31;mFitdErrs._cTT-=k40*c40+k41*c41; 00946 } 00947 if (mFitdErrs.check()) return -12; 00948 // mFitdErrs.recov(); 00949 00950 00951 static int ERRTEST=0; 00952 if(ERRTEST) errTest(mPredPars,mPredErrs,mHit,mHrr,mFitdPars,mFitdErrs,mChi2); 00953 00954 //prod assert(mHrr.hYY > mFitdErrs._cYY); 00955 //prod assert(mHrr.hZZ > mFitdErrs._cZZ); 00956 if (mDetector) { //Not a primary 00957 if (mHrr.hYY <= mFitdErrs._cYY) { 00958 LOG_DEBUG << Form("StiTrackNodeHelper::updateNode() WRONG hYY(%g) < nYY(%g)" 00959 ,mHrr.hYY,mFitdErrs._cYY)<< endm; 00960 return -13; 00961 } 00962 if (mHrr.hZZ <= mFitdErrs._cZZ) { 00963 LOG_DEBUG << Form("StiTrackNodeHelper::updateNode() WRONG hZZ(%g) < nZZ(%g)" 00964 ,mHrr.hZZ,mFitdErrs._cZZ)<< endm; 00965 return -14; 00966 } 00967 } //EndIf Not a primary 00968 if (mTargetNode && debug() & 8) mTargetNode->PrintpT("U"); 00969 mState = StiTrackNode::kTNFitEnd; 00970 return 0; 00971 } 00972 00973 //______________________________________________________________________________ 00974 void StiTrackNodeHelper::resetError(double fk) 00975 { 00976 if (fk) do {//fake loop 00977 if(mPredErrs._cYY>DY*DY) break; 00978 if(mPredErrs._cZZ>DZ*DZ) break; 00979 if(mPredErrs._cEE>DEta*DEta) break; 00980 if(mPredErrs._cPP>DPti*DPti) break; 00981 if(mPredErrs._cTT>DTan*DTan) break; 00982 mPredErrs*=fk; 00983 return; 00984 }while(0); 00985 00986 mPredErrs.reset(); 00987 mPredErrs._cYY=DY*DY; 00988 mPredErrs._cZZ=DZ*DZ; 00989 mPredErrs._cEE=DEta*DEta; 00990 mPredErrs._cPP=DPti*DPti; 00991 mPredErrs._cTT=DTan*DTan; 00992 } 00993 //______________________________________________________________________________ 00994 int StiTrackNodeHelper::cutStep(StiNodePars *pars,StiNodePars *base) 00995 { 00996 double fact=1,dif,fak; 00997 double cuts[kNPars]; 00998 memcpy(cuts,MAXSTEP,sizeof(cuts)); 00999 if (mBestDelta<DY) {cuts[1]=mBestDelta;cuts[2]=mBestDelta;} 01000 01001 for (int jx=0;jx<kNPars;jx++) { 01002 dif =(fabs((*pars)[jx]-(*base)[jx])); 01003 fak = (dif >cuts[jx]) ? cuts[jx]/dif:1; 01004 if (fak < fact) fact = fak; 01005 } 01006 if (fact>=1) return 0; 01007 mgCutStep++; 01008 for (int jx=0;jx<kNPars;jx++) { 01009 dif =(*pars)[jx]-(*base)[jx]; 01010 (*pars)[jx] = (*base)[jx] +dif*fact; 01011 } 01012 pars->ready(); 01013 return 1; 01014 } 01015 //______________________________________________________________________________ 01016 int StiTrackNodeHelper::nudge() 01017 { 01018 if(!mHit) return 0; 01019 StiNodePars *pars = &mBestPars; 01020 for (int i=0;i<2;i++,pars =&mPredPars) { 01021 double deltaX = mHitPars[0]-pars->_x; 01022 if (fabs(deltaX) <1e-6) continue; 01023 double deltaL = deltaX/pars->_cosCA; 01024 double deltaE = pars->_curv*deltaL; 01025 pars->_x = mHitPars[0]; 01026 pars->_y += pars->_sinCA *deltaL; 01027 pars->_z += pars->_tanl *deltaL; 01028 pars->_eta += deltaE; 01029 pars->_cosCA -= pars->_sinCA *deltaE; 01030 pars->_sinCA += pars->_cosCA *deltaE; 01031 if (pars->_cosCA>=1.) pars->ready(); 01032 if (pars->check()) return 1; 01033 } 01034 return 0; 01035 } 01036 //______________________________________________________________________________ 01037 double StiTrackNodeHelper::pathIn(const StiDetector *det,StiNodePars *pars) 01038 { 01039 if (!det) return 0.; 01040 double thickness = det->getShape()->getThickness(); 01041 double t = pars->_tanl; 01042 double c = pars->_cosCA; 01043 return (thickness*::sqrt(1.+t*t)) / c; 01044 } 01045 //______________________________________________________________________________ 01046 int StiTrackNodeHelper::getHitErrors(const StiHit *hit,const StiNodePars *pars,StiHitErrs *hrr) 01047 { 01048 hrr->reset(); 01049 const StiDetector *det = hit->detector(); 01050 const StiHitErrorCalculator *calc = (det)? det->getHitErrorCalculator():0; 01051 if (calc) {//calculate it 01052 calc->calculateError(pars,hrr->hYY,hrr->hZZ); 01053 } else {//get from hit 01054 const float *ermx = hit->errMtx(); 01055 for (int i=0;i<6;i++){hrr->A[i]=ermx[i];} 01056 } 01057 return (!det); 01058 } 01059 //______________________________________________________________________________ 01060 int errTest(StiNodePars &predP,StiNodeErrs &predE, 01061 const StiHit *hit ,StiHitErrs &hitErr, 01062 StiNodePars &fitdP,StiNodeErrs &fitdE, double chi2) 01063 { 01064 01065 StiNodePars mineP,hittP; 01066 StiNodeErrs mineE,hittE; 01067 hittP._x = hit->x(); 01068 hittP._y = hit->y(); 01069 hittP._z = hit->z(); 01070 memcpy(hittE.A,hitErr.A,sizeof(StiNodeErrs)); 01071 01072 double myChi2 = StiTrackNodeHelper::joinTwo( 01073 3,hittP.P,hittE.A, 01074 6,predP.P,predE.A, 01075 mineP.P,mineE.A); 01076 01077 01078 int ndif = 0; 01079 for (int i=0;i<kNPars;i++) { 01080 double diff = fabs(mineP.P[i]-fitdP.P[i]); 01081 if (diff < 1e-10) continue; 01082 diff/=0.5*(fabs(mineP.P[i])+fabs(fitdP.P[i])); 01083 if (diff < 1e-5 ) continue; 01084 ndif++; 01085 LOG_DEBUG << Form("errTest(P): %g(%d) - %g(%d) = %g",mineE.A[i],i,fitdE.A[i],i,diff)<< endm; 01086 } 01087 if (ndif){ mineP.print();fitdP.print();} 01088 01089 for (int i=0;i<kNErrs;i++) { 01090 double diff = fabs(mineE.A[i]-fitdE.A[i]); 01091 if (diff < 1e-10) continue; 01092 diff/=0.5*(fabs(mineE.A[i])+fabs(fitdE.A[i])); 01093 if (diff < 1e-5 ) continue; 01094 ndif+=100; 01095 LOG_DEBUG << Form("errTest(E): %g(%d) - %g(%d) = %g",mineE.A[i],i,fitdE.A[i],i,diff)<< endm; 01096 } 01097 if (ndif>=100){ mineE.print();fitdE.print();} 01098 01099 double diff = fabs((chi2-myChi2)/(chi2+myChi2)); 01100 if (diff > 1e-5 ) { 01101 ndif+=1000; 01102 LOG_DEBUG << Form("errTest(C): %g() - %g() = %g",myChi2,chi2,diff)<< endm; 01103 } 01104 01105 return ndif; 01106 } 01107 01108 //_____________________________________________________________________________ 01109 void QaFit::add(double val) 01110 { 01111 double v = val; double p = mPrev; mPrev = v; 01112 int n = (mTally)? 2:1; 01113 mTally++; 01114 for (int i=0;i<n;i++) { 01115 mAver[i] +=v; 01116 mErrr[i] +=v*v; 01117 if (mMaxi[i]<fabs(v)) mMaxi[i]=fabs(v); 01118 if (v<0) mNega[i]++; 01119 v = v*p; 01120 } 01121 } 01122 //_____________________________________________________________________________ 01123 void QaFit::finish() 01124 { 01125 if( mEnded) return; 01126 if(!mTally) return; 01127 mEnded = 1; 01128 int n = mTally; 01129 for (int i=0;i<2;i++) { 01130 mAver[i]/= n; 01131 mErrr[i]/= n; 01132 if (mErrr[i]<1e-10) mErrr[i]=1e-10; 01133 mErrr[i] = sqrt(mErrr[i]); 01134 if (!(--n)) break; 01135 } 01136 } 01137 //_____________________________________________________________________________ 01138 double QaFit::getAccu(int k) 01139 { 01140 finish(); 01141 if (mTally-k<=1) return 0; 01142 return mErrr[k]; 01143 } 01144 //_____________________________________________________________________________ 01145 double QaFit::getNStd(int k) 01146 { 01147 finish(); 01148 int n = mTally-k; 01149 if (!n) return 0; 01150 return mAver[k]*sqrt(double(n)); 01151 } 01152 //_____________________________________________________________________________ 01153 double QaFit::getNSgn(int k) 01154 { 01155 finish(); 01156 int n = mTally-k; 01157 if (!n) return 0; 01158 return (n-2*mNega[k])/sqrt(double(n)); 01159 } 01160 //_____________________________________________________________________________ 01161 void QaFit::getInfo(double *info) 01162 { 01163 for (int i=0;i<2;i++) { 01164 int l = i*4; 01165 info[l+0]=getAccu(i); 01166 info[l+1]=getNStd(i); 01167 info[l+2]=getNSgn(i); 01168 info[l+3]=getMaxi(i); 01169 } 01170 } 01171 01172 01173