#include <TCosmicFitter.h>
Inheritance diagram for TCosmicFitter:
Public Member Functions | |
| TCosmicFitter (const std::string &name) | |
| Constructor. | |
| virtual | ~TCosmicFitter () |
| Destructor. | |
| void | dump (const std::string &message=std::string(""), const std::string &prefix=std::string("")) const |
| dumps debug information. | |
| int | fit (TTrackBase &) const |
| int | fit (TTrackBase &, float t0Offset) const |
| int | fitWithCathode (TTrackBase &, float t0Offset=0., float windowSize=0.6, int SysCorr=0) |
| const std::string & | name (void) const |
| returns name. | |
Protected Member Functions | |
| void | fitDone (TTrackBase &) const |
| sets the fitted flag. (Bad implementation) | |
Private Member Functions | |
| int | dxda (const TMLink &link, const Helix &helix, double dPhi, HepVector &dxda, HepVector &dyda, HepVector &dzda, int doSagCorrection) const |
| calculates dXda. 'link' and 'dPhi' are inputs. Others are outputs. | |
Private Attributes | |
| IMagneticFieldSvc * | m_pmgnIMF |
Definition at line 45 of file TCosmicFitter.h.
| TCosmicFitter::TCosmicFitter | ( | const std::string & | name | ) |
Constructor.
Definition at line 84 of file TCosmicFitter.cxx.
References m_pmgnIMF.
00084 : TMFitter(name) { 00085 //jialk 00086 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF); 00087 if(scmgn!=StatusCode::SUCCESS) { 00088 std::cout<< "Unable to open Magnetic field service"<<std::endl; 00089 } 00090 00091 }
| TCosmicFitter::~TCosmicFitter | ( | ) | [virtual] |
| void TCosmicFitter::dump | ( | const std::string & | message = std::string(""), |
|
| const std::string & | prefix = std::string("") | |||
| ) | const |
| int TCosmicFitter::fit | ( | TTrackBase & | , | |
| float | t0Offset | |||
| ) | const |
Definition at line 113 of file TCosmicFitter.cxx.
References Convergence, DBL_MAX, check_raw_filter::dist, IMdcCalibFunSvc::driftTimeToDist(), dxda(), showlog::err, IMdcCalibFunSvc::getSigma(), IMdcCalibFunSvc::getT0(), IMdcCalibFunSvc::getTimeWalk(), genRecEmupikp::i, ganga-rec::j, TTrackBase::links(), msgSvc(), NTrailMax, TTrackBase::objectType(), t(), TFitErrorFewHits, TFitFailed, TFitUnavailable, Track, v, Bes_Common::WARNING, WireHitFittingValid, WireHitInvalidForFit, WireHitLeft, WireHitRight, and x.
00113 { 00114 00115 #ifdef TRKRECO_DEBUG_DETAIL 00116 std::cout << " TCosmicFitter::fit ..." << std::endl; 00117 #endif 00118 00119 IMdcCalibFunSvc* l_mdcCalFunSvc; 00120 Gaudi::svcLocator() -> service("MdcCalibFunSvc", l_mdcCalFunSvc); 00121 00122 //...Check # of hits... 00123 if (b.links().length() < 5) return TFitErrorFewHits; 00124 unsigned nValid = 0; 00125 for (unsigned i = 0; i < b.links().length(); i++) { 00126 unsigned state = b.links()[i]->hit()->state(); 00127 if (state & WireHitInvalidForFit) continue; 00128 if (state & WireHitFittingValid) ++nValid; 00129 } 00130 if (nValid < 5) return TFitErrorFewHits; 00131 00132 //...Type check... 00133 // if (b.type() != Track) return TFitUnavailable; 00134 if (b.objectType() != Track) return TFitUnavailable; 00135 TTrack & t = (TTrack &) b; 00136 00137 //read t0 from TDS-------------------------------------// 00138 double _t0_bes = -1.; 00139 IMessageSvc *msgSvc; 00140 Gaudi::svcLocator()->service("MessageSvc", msgSvc); 00141 MsgStream log(msgSvc, "TCosmicFitter"); 00142 00143 IDataProviderSvc* eventSvc = NULL; 00144 Gaudi::svcLocator()->service("EventDataSvc", eventSvc); 00145 00146 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc,"/Event/Recon/RecEsTimeCol"); 00147 if (aevtimeCol) { 00148 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin(); 00149 _t0_bes = (*iter_evt)->getTest(); 00150 // cout<<" tev: "<<setw(4)<<_t0_bes<<endl; 00151 }else{ 00152 log << MSG::WARNING << "Could not find RecEsTimeCol" << endreq; 00153 } 00154 // if (_t0_bes < 7.) _t0_bes = 7.; 00155 //----------------------------------------------------// 00156 00157 //...Setup... 00158 Vector a(5), da(5); 00159 a = t.helix().a(); 00160 Vector dxda(5); 00161 Vector dyda(5); 00162 Vector dzda(5); 00163 Vector dDda(5); 00164 Vector dchi2da(5); 00165 SymMatrix d2chi2d2a(5, 0); 00166 double chi2; 00167 // double chi2Old = 10e99; 00168 double chi2Old = DBL_MAX; 00169 const double convergence = Convergence; 00170 bool allAxial = true; 00171 Matrix e(3, 3); 00172 Vector f(3); 00173 int err = 0; 00174 double factor = 1.0;//jtanaka0715 00175 00176 Vector maxDouble(5); 00177 for (unsigned i = 0; i < 5; i++) maxDouble[i] = (FLT_MAX); 00178 00179 //...Fitting loop... 00180 unsigned nTrial = 0; 00181 while (nTrial < NTrailMax) { 00182 00183 //...Set up... 00184 chi2 = 0.; 00185 for (unsigned j = 0; j < 5; j++) dchi2da[j] = 0.; 00186 d2chi2d2a = SymMatrix(5, 0); 00187 00188 #if CAL_TOF_HELIX 00189 //cal the time pass through the chamber 00190 const float Rad = 81; // cm 00191 float dRho = t.helix().a()[0]; 00192 float Lproj = sqrt(Rad*Rad - dRho*dRho); 00193 float tlmd = t.helix().a()[4]; 00194 float fct = sqrt(1. + tlmd * tlmd); 00195 #endif 00196 00197 //...Loop with hits... 00198 unsigned i = 0; 00199 while (TMLink * l = t.links()[i++]) { 00200 const TMDCWireHit & h = * l->hit(); 00201 00202 //...Check state... 00203 if (h.state() & WireHitInvalidForFit) continue; 00204 if (! (h.state() & WireHitFittingValid)) continue; 00205 00206 //...Check wire... 00207 if (! nTrial) 00208 if (h.wire()->stereo()) allAxial = false; 00209 00210 //...Cal. closest points... 00211 int doSagCorrection = 0; 00212 // if(nTrial && !allAxial ) doSagCorrection = 1; //Liuqg 00213 t.approach(* l, doSagCorrection); 00214 double dPhi = l->dPhi(); 00215 const HepPoint3D & onTrack = l->positionOnTrack(); 00216 const HepPoint3D & onWire = l->positionOnWire(); 00217 00218 #ifdef TRKRECO_DEBUG_DETAIL 00219 // std::cout << " in fit " << onTrack << ", " << onWire; 00220 // h.dump(); 00221 #endif 00222 00223 //...Obtain drift distance and its error... 00224 unsigned leftRight = WireHitRight; 00225 if (onWire.cross(onTrack).z() < 0.) leftRight = WireHitLeft; 00226 double distance = h.drift(leftRight); 00227 double eDistance = h.dDrift(leftRight); 00228 //... 00229 if(nTrial && !allAxial){ 00230 int side = leftRight; 00231 00232 double Tfly = _t0_bes/220.*(110.-onWire.y()); 00233 #if CAL_TOF_HELIX 00234 float Rad_wire = sqrt(x[0]*x[0]+x[1]*x[1]); // cm 00235 float L_wire = sqrt(Rad_wire*Rad_wire - dRho*dRho); 00236 float flyLength = Lproj - L_wire; 00237 if (x[1]<0) flyLength = Lproj + L_wire; 00238 Tfly = flyLength/29.98*sqrt(1.0+(0.105/0.5)*(0.105/0.5))*fct; //approxiate... cm/ns 00239 #endif 00240 float time = l->hit()->reccdc()->tdc - Tfly; 00241 00242 int wire = h.wire()->localId(); 00243 int layerId = h.wire()->layerId(); 00244 float dist = distance; 00245 float edist = eDistance; 00246 double T0 = l_mdcCalFunSvc->getT0(layerId,wire); 00247 double rawadc = 0.; 00248 rawadc =l->hit()->reccdc()->adc; 00249 00250 double timeWalk = l_mdcCalFunSvc->getTimeWalk(layerId, rawadc); 00251 double drifttime =time -T0-timeWalk; 00252 dist = l_mdcCalFunSvc->driftTimeToDist(drifttime,layerId, wire, side, 0.0); 00253 edist = l_mdcCalFunSvc->getSigma(layerId, side, dist, 0.0); 00254 dist = dist/10.0; //mm->cm 00255 edist = edist/10.0; 00256 // if(layerId<20) edist = 9999.; 00257 00258 //zsl calcdc_driftdist_(&doPropDelay,&wire,&side,p,x,&time,&dist,&edist); 00259 00260 distance = (double) dist; 00261 eDistance = (double) edist; 00262 00263 l->drift(distance,0); //update 00264 l->drift(distance,1); 00265 l->dDrift(eDistance,0); 00266 l->dDrift(eDistance,1); 00267 l->tof(Tfly); 00268 } 00269 double eDistance2 = eDistance * eDistance; 00270 00271 //...Residual... 00272 HepVector3D v = onTrack - onWire; 00273 double vmag = v.mag(); 00274 double dDistance = vmag - distance; 00275 00276 //...dxda... 00277 this->dxda(*l, t.helix(), dPhi, dxda, dyda, dzda, doSagCorrection); 00278 00279 //...Chi2 related... 00280 // dDda = (v.x() * dxda + v.y() * dyda + v.z() * dzda) / vmag; 00281 HepVector3D vw = h.wire()->direction(); 00282 dDda = (vmag > 0.) 00283 ? ((v.x() * (1. - vw.x() * vw.x()) - 00284 v.y() * vw.x() * vw.y() - v.z() * vw.x() * vw.z()) 00285 * dxda + 00286 (v.y() * (1. - vw.y() * vw.y()) - 00287 v.z() * vw.y() * vw.z() - v.x() * vw.y() * vw.x()) 00288 * dyda + 00289 (v.z() * (1. - vw.z() * vw.z()) - 00290 v.x() * vw.z() * vw.x() - v.y() * vw.z() * vw.y()) 00291 * dzda) / vmag 00292 : Vector(5,0); 00293 if(vmag<=0.0) { 00294 std::cout << " in fit " << onTrack << ", " << onWire; 00295 h.dump(); 00296 } 00297 dchi2da += (dDistance / eDistance2) * dDda; 00298 d2chi2d2a += vT_times_v(dDda) / eDistance2; 00299 double pChi2 = dDistance * dDistance / eDistance2; 00300 chi2 += pChi2; 00301 00302 //...Store results... 00303 l->update(onTrack, onWire, leftRight, pChi2); 00304 } 00305 00306 //...Check condition... 00307 double change = chi2Old - chi2; 00308 if (fabs(change) < convergence) break; 00309 //if (change < 0.) break; 00310 //jtanaka -- from traffs -- Ozaki-san added this part to traffs. 00311 if (change < 0.){ 00312 #ifdef TRKRECO_DEBUG_DETAIL 00313 std::cout << "chi2Old, chi2=" << chi2Old <<" "<< chi2 << std::endl; 00314 #endif 00315 //change to the old value. 00316 a += factor*da; 00317 // a[2] = 0.000000001; 00318 t._helix->a(a); 00319 00320 chi2 = 0.; 00321 for (unsigned j = 0; j < 5; j++) dchi2da[j] = 0.; 00322 d2chi2d2a = SymMatrix(5, 0); 00323 00324 //...Loop with hits... 00325 unsigned i = 0; 00326 while (TMLink * l = t.links()[i++]) { 00327 const TMDCWireHit & h = * l->hit(); 00328 00329 //...Check state... 00330 if (h.state() & WireHitInvalidForFit) continue; 00331 if (! (h.state() & WireHitFittingValid)) continue; 00332 00333 //...Check wire... 00334 if (! nTrial) 00335 if (h.wire()->stereo()) allAxial = false; 00336 00337 //...Cal. closest points... 00338 int doSagCorrection = 0; 00339 // if( nTrial && !allAxial ) doSagCorrection = 1; //Liuqg 00340 t.approach(* l, doSagCorrection); 00341 double dPhi = l->dPhi(); 00342 const HepPoint3D & onTrack = l->positionOnTrack(); 00343 const HepPoint3D & onWire = l->positionOnWire(); 00344 00345 #ifdef TRKRECO_DEBUG_DETAIL 00346 // std::cout << " in fit in case of change < 0. " << onTrack << ", " << onWire; 00347 // h.dump(); 00348 #endif 00349 00350 //...Obtain drift distance and its error... 00351 unsigned leftRight = WireHitRight; 00352 if (onWire.cross(onTrack).z() < 0.) leftRight = WireHitLeft; 00353 double distance = h.drift(leftRight); 00354 double eDistance = h.dDrift(leftRight); 00355 if(nTrial && !allAxial){ 00356 int side = leftRight; 00357 00358 double Tfly = _t0_bes/220.*(110.-onWire.y()); 00359 #if CAL_TOF_HELIX 00360 float Rad_wire = sqrt(x[0]*x[0]+x[1]*x[1]); // cm 00361 float L_wire = sqrt(Rad_wire*Rad_wire - dRho*dRho); 00362 float flyLength = Lproj - L_wire; 00363 if (x[1]<0) flyLength = Lproj + L_wire; 00364 Tfly = flyLength/29.98*sqrt(1.0+(0.105/10)*(0.105/10))*fct; //approxiate... cm/ns 00365 #endif 00366 00367 float time = l->hit()->reccdc()->tdc - Tfly; 00368 00369 int wire = h.wire()->localId(); 00370 int layerId = h.wire()->layerId(); 00371 float dist= distance; 00372 float edist = eDistance; 00373 // int doPropDelay = 1; // 00374 00375 double T0 = l_mdcCalFunSvc->getT0(layerId,wire); 00376 double rawadc = 0.; 00377 rawadc =l->hit()->reccdc()->adc; 00378 double timeWalk = l_mdcCalFunSvc->getTimeWalk(layerId, rawadc); 00379 double drifttime =time -T0-timeWalk; 00380 dist = l_mdcCalFunSvc->driftTimeToDist(drifttime,layerId, wire, side, 0.0); 00381 edist = l_mdcCalFunSvc->getSigma(layerId, side, dist, 0.0); 00382 dist = dist/10.0; //mm->cm 00383 edist = edist/10.0; 00384 // if (layerId<20) edist = 9999.; 00385 00386 //zsl calcdc_driftdist_(&doPropDelay,&wire,&side,p,x,&time,&dist,&edist); 00387 00388 distance = (double) dist; 00389 eDistance = (double) edist; 00390 00391 l->drift(distance,0); //update 00392 l->drift(distance,1); 00393 l->dDrift(eDistance,0); 00394 l->dDrift(eDistance,1); 00395 l->tof(Tfly); 00396 } 00397 double eDistance2 = eDistance * eDistance; 00398 00399 //...Residual... 00400 HepVector3D v = onTrack - onWire; 00401 double vmag = v.mag(); 00402 double dDistance = vmag - distance; 00403 00404 //...dxda... 00405 this->dxda(*l, t.helix(), dPhi, dxda, dyda, dzda, doSagCorrection); 00406 00407 //...Chi2 related... 00408 //dDda = (v.x() * dxda + v.y() * dyda + v.z() * dzda) / vmag; 00409 HepVector3D vw = h.wire()->direction(); 00410 dDda = (vmag > 0.) 00411 ? ((v.x() * (1. - vw.x() * vw.x()) - 00412 v.y() * vw.x() * vw.y() - v.z() * vw.x() * vw.z()) 00413 * dxda + 00414 (v.y() * (1. - vw.y() * vw.y()) - 00415 v.z() * vw.y() * vw.z() - v.x() * vw.y() * vw.x()) 00416 * dyda + 00417 (v.z() * (1. - vw.z() * vw.z()) - 00418 v.x() * vw.z() * vw.x() - v.y() * vw.z() * vw.y()) 00419 * dzda) / vmag 00420 : Vector(5,0); 00421 if(vmag<=0.0) { 00422 std::cout << " in fit " << onTrack << ", " << onWire; 00423 h.dump(); 00424 } 00425 dchi2da += (dDistance / eDistance2) * dDda; 00426 d2chi2d2a += vT_times_v(dDda) / eDistance2; 00427 double pChi2 = dDistance * dDistance / eDistance2; 00428 chi2 += pChi2; 00429 00430 //...Store results... 00431 l->update(onTrack, onWire, leftRight, pChi2); 00432 } 00433 //break; 00434 factor *= 0.75; 00435 #ifdef TRKRECO_DEBUG_DETAIL 00436 std::cout << "factor = " << factor << std::endl; 00437 std::cout << "chi2 = " << chi2 << std::endl; 00438 #endif 00439 if(factor < 0.01)break; 00440 } 00441 00442 chi2Old = chi2; 00443 00444 //...Cal. helix parameters for next loop... 00445 if (allAxial) { 00446 f = dchi2da.sub(1, 3); 00447 e = d2chi2d2a.sub(1, 3); 00448 f = solve(e, f); 00449 da[0] = f[0]; 00450 da[1] = f[1]; 00451 da[2] = f[2]; 00452 da[3] = 0.; 00453 da[4] = 0.; 00454 } 00455 else { 00456 da = solve(d2chi2d2a, dchi2da); 00457 } 00458 00459 #ifdef TRKRECO_DEBUG_DETAIL 00460 //std::cout << " fit " << nTrial << " : " << da << std::endl; 00461 //std::cout << " d2chi " << d2chi2d2a << std::endl; 00462 //std::cout << " dchi2 " << dchi2da << std::endl; 00463 #endif 00464 00465 a -= factor*da; 00466 // a[2] = 0.000000001; 00467 t._helix->a(a); 00468 ++nTrial; 00469 } 00470 00471 //...Cal. error matrix... 00472 SymMatrix Ea(5, 0); 00473 unsigned dim; 00474 if (allAxial) { 00475 dim = 3; 00476 SymMatrix Eb(3, 0), Ec(3, 0); 00477 Eb = d2chi2d2a.sub(1, 3); 00478 Ec = Eb.inverse(err); 00479 Ea[0][0] = Ec[0][0]; 00480 Ea[0][1] = Ec[0][1]; 00481 Ea[0][2] = Ec[0][2]; 00482 Ea[1][1] = Ec[1][1]; 00483 Ea[1][2] = Ec[1][2]; 00484 Ea[2][2] = Ec[2][2]; 00485 } 00486 else { 00487 dim = 5; 00488 Ea = d2chi2d2a.inverse(err); 00489 } 00490 00491 //...Store information... 00492 if (! err) { 00493 t._helix->a(a); 00494 t._helix->Ea(Ea); 00495 } 00496 else { 00497 err = TFitFailed; 00498 } 00499 00500 t._ndf = nValid - dim; 00501 t._chi2 = chi2; 00502 // t._fitted = true; 00503 00504 return err; 00505 }
| int TCosmicFitter::fit | ( | TTrackBase & | ) | const [virtual] |
Implements TMFitter.
Definition at line 97 of file TCosmicFitter.cxx.
References TTrackBase::_fitted, showlog::err, TTrackBase::fitted(), and TFitAlreadyFitted.
Referenced by TBuilderCosmic::buildStereo(), and TBuilder0::fit().
00097 { 00098 00099 // double t0_bes; 00100 // TrkReco::t0(t0_bes); 00101 // const double t0_bes = TrkReco::t0(); 00102 00103 //...Already fitted ?... 00104 if (b.fitted()) return TFitAlreadyFitted; 00105 00106 int err = fit(b, 0.); 00107 if (! err) b._fitted = true; 00108 00109 return err; 00110 }
| void TMFitter::fitDone | ( | TTrackBase & | ) | const [protected, inherited] |
sets the fitted flag. (Bad implementation)
Definition at line 24 of file TMFitter.cxx.
References t().
Referenced by TRobustLineFitter::fit(), TLineFitter::fit(), and TCircleFitter::fit().
00024 { 00025 t._fitted = true; 00026 }
| int TCosmicFitter::fitWithCathode | ( | TTrackBase & | , | |
| float | t0Offset = 0., |
|||
| float | windowSize = 0.6, |
|||
| int | SysCorr = 0 | |||
| ) |
| const std::string & TMFitter::name | ( | void | ) | const [inline, inherited] |
returns name.
Definition at line 73 of file TMFitter.h.
References TMFitter::_name.
00073 { 00074 return _name; 00075 }
IMagneticFieldSvc* TCosmicFitter::m_pmgnIMF [private] |
1.4.7