BOSS 7.0.7
BESIII Offline Software System
Loading...
Searching...
No Matches
TCosmicFitter.cxx
Go to the documentation of this file.
1//-----------------------------------------------------------------------------
2// $Id: TCosmicFitter.cxx,v 1.12 2011/12/05 04:49:50 maoh Exp $
3//-----------------------------------------------------------------------------
4// Filename : TCosmicFitter.cc
5// Section : Tracking
6// Owner : Yoshi Iwasaki
7// Email : [email protected]
8//-----------------------------------------------------------------------------
9// Description : A class to fit a TTrackBase object to a helix.
10// See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
11//-----------------------------------------------------------------------------
12
13#define HEP_SHORT_NAMES
14#include <float.h>
15//#include "tables/cdc.h"
16#include "CLHEP/Matrix/Vector.h"
17#include "CLHEP/Matrix/SymMatrix.h"
18#include "CLHEP/Matrix/DiagMatrix.h"
19#include "CLHEP/Matrix/Matrix.h"
20//#include "panther/panther.h"
21//#include "panther/panther_manager.h"
22#include "TrkReco/TCosmicFitter.h"
23#include "TrkReco/TTrack.h"
24
25#include "TrkReco/TrkReco.h"
26
27#include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
28#include "MdcCalibFunSvc/MdcCalibFunSvc.h"
29
30#include "GaudiKernel/MsgStream.h"
31#include "GaudiKernel/AlgFactory.h"
32#include "GaudiKernel/ISvcLocator.h"
33#include "GaudiKernel/IMessageSvc.h"
34#include "GaudiKernel/IDataProviderSvc.h"
35#include "GaudiKernel/IDataManagerSvc.h"
36#include "GaudiKernel/SmartDataPtr.h"
37#include "GaudiKernel/IDataProviderSvc.h"
38#include "GaudiKernel/PropertyMgr.h"
39#include "GaudiKernel/IJobOptionsSvc.h"
40#include "GaudiKernel/IMessageSvc.h"
41#include "GaudiKernel/Bootstrap.h"
42#include "GaudiKernel/StatusCode.h"
43
44#include "GaudiKernel/ContainedObject.h"
45#include "GaudiKernel/SmartRef.h"
46#include "GaudiKernel/SmartRefVector.h"
47#include "GaudiKernel/ObjectVector.h"
48#include "EventModel/EventModel.h"
49
50#include "EvTimeEvent/RecEsTime.h"
51
52#include "CLHEP/Matrix/Matrix.h"
53#include "GaudiKernel/StatusCode.h"
54#include "GaudiKernel/IInterface.h"
55#include "GaudiKernel/Kernel.h"
56#include "GaudiKernel/Service.h"
57#include "GaudiKernel/ISvcLocator.h"
58#include "GaudiKernel/SvcFactory.h"
59#include "GaudiKernel/IDataProviderSvc.h"
60#include "GaudiKernel/Bootstrap.h"
61#include "GaudiKernel/MsgStream.h"
62#include "GaudiKernel/SmartDataPtr.h"
63#include "GaudiKernel/IMessageSvc.h"
64
65
66using CLHEP::HepVector;
67using CLHEP::HepSymMatrix;
68using CLHEP::HepDiagMatrix;
69using CLHEP::HepMatrix;
70
71typedef CLHEP::HepMatrix Matrix;
72
73#define CAL_TOF_HELIX 0
74
75#define NTrailMax 100
76#define Convergence 1.0e-5
77
78/*
79extern "C" void calcdc_sag3_(int*, float*, float[3], float*, float*, float*);
80extern "C" void calcdc_driftdist_(int*, int*, int*,
81 float[3], float[3], float*, float*, float*);
82*/
83
84TCosmicFitter::TCosmicFitter(const std::string & name) : TMFitter(name) {
85//jialk
86 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
87 if(scmgn!=StatusCode::SUCCESS) {
88 std::cout<< "Unable to open Magnetic field service"<<std::endl;
89 }
90
91}
92
93TCosmicFitter::~TCosmicFitter() {
94}
95
96int
97TCosmicFitter::fit(TTrackBase & b) const {
98
99// double t0_bes;
100// TrkReco::t0(t0_bes);
101// const double t0_bes = TrkReco::t0();
102
103 //...Already fitted ?...
104 if (b.fitted()) return TFitAlreadyFitted;
105
106 int err = fit(b, 0.);
107 if (! err) b._fitted = true;
108
109 return err;
110}
111
112int
113TCosmicFitter::fit(TTrackBase & b, float t0Offset) const {
114
115#ifdef TRKRECO_DEBUG_DETAIL
116 std::cout << " TCosmicFitter::fit ..." << std::endl;
117#endif
118
119 IMdcCalibFunSvc* l_mdcCalFunSvc;
120 Gaudi::svcLocator() -> service("MdcCalibFunSvc", l_mdcCalFunSvc);
121
122 //...Check # of hits...
123 if (b.links().length() < 5) return TFitErrorFewHits;
124 unsigned nValid = 0;
125 for (unsigned i = 0; i < b.links().length(); i++) {
126 unsigned state = b.links()[i]->hit()->state();
127 if (state & WireHitInvalidForFit) continue;
128 if (state & WireHitFittingValid) ++nValid;
129 }
130 if (nValid < 5) return TFitErrorFewHits;
131
132 //...Type check...
133 // if (b.type() != Track) return TFitUnavailable;
134 if (b.objectType() != Track) return TFitUnavailable;
135 TTrack & t = (TTrack &) b;
136
137 //read t0 from TDS-------------------------------------//
138 double _t0_bes = -1.;
139 IMessageSvc *msgSvc;
140 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
141 MsgStream log(msgSvc, "TCosmicFitter");
142
143 IDataProviderSvc* eventSvc = NULL;
144 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
145
146 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc,"/Event/Recon/RecEsTimeCol");
147 if (aevtimeCol) {
148 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
149 _t0_bes = (*iter_evt)->getTest();
150// cout<<" tev: "<<setw(4)<<_t0_bes<<endl;
151 }else{
152 log << MSG::WARNING << "Could not find RecEsTimeCol" << endreq;
153 }
154// if (_t0_bes < 7.) _t0_bes = 7.;
155 //----------------------------------------------------//
156
157 //...Setup...
158 Vector a(5), da(5);
159 a = t.helix().a();
160 Vector dxda(5);
161 Vector dyda(5);
162 Vector dzda(5);
163 Vector dDda(5);
164 Vector dchi2da(5);
165 SymMatrix d2chi2d2a(5, 0);
166 double chi2;
167 // double chi2Old = 10e99;
168 double chi2Old = DBL_MAX;
169 const double convergence = Convergence;
170 bool allAxial = true;
171 Matrix e(3, 3);
172 Vector f(3);
173 int err = 0;
174 double factor = 1.0;//jtanaka0715
175
176 Vector maxDouble(5);
177 for (unsigned i = 0; i < 5; i++) maxDouble[i] = (FLT_MAX);
178
179 //...Fitting loop...
180 unsigned nTrial = 0;
181 while (nTrial < NTrailMax) {
182
183 //...Set up...
184 chi2 = 0.;
185 for (unsigned j = 0; j < 5; j++) dchi2da[j] = 0.;
186 d2chi2d2a = SymMatrix(5, 0);
187
188#if CAL_TOF_HELIX
189 //cal the time pass through the chamber
190 const float Rad = 81; // cm
191 float dRho = t.helix().a()[0];
192 float Lproj = sqrt(Rad*Rad - dRho*dRho);
193 float tlmd = t.helix().a()[4];
194 float fct = sqrt(1. + tlmd * tlmd);
195#endif
196
197 //...Loop with hits...
198 unsigned i = 0;
199 while (TMLink * l = t.links()[i++]) {
200 const TMDCWireHit & h = * l->hit();
201
202 //...Check state...
203 if (h.state() & WireHitInvalidForFit) continue;
204 if (! (h.state() & WireHitFittingValid)) continue;
205
206 //...Check wire...
207 if (! nTrial)
208 if (h.wire()->stereo()) allAxial = false;
209
210 //...Cal. closest points...
211 int doSagCorrection = 0;
212// if(nTrial && !allAxial ) doSagCorrection = 1; //Liuqg
213 t.approach(* l, doSagCorrection);
214 double dPhi = l->dPhi();
215 const HepPoint3D & onTrack = l->positionOnTrack();
216 const HepPoint3D & onWire = l->positionOnWire();
217
218#ifdef TRKRECO_DEBUG_DETAIL
219// std::cout << " in fit " << onTrack << ", " << onWire;
220// h.dump();
221#endif
222
223 //...Obtain drift distance and its error...
224 unsigned leftRight = WireHitRight;
225 if (onWire.cross(onTrack).z() < 0.) leftRight = WireHitLeft;
226 double distance = h.drift(leftRight);
227 double eDistance = h.dDrift(leftRight);
228 //...
229 if(nTrial && !allAxial){
230 int side = leftRight;
231
232 double Tfly = _t0_bes/220.*(110.-onWire.y());
233#if CAL_TOF_HELIX
234 float Rad_wire = sqrt(x[0]*x[0]+x[1]*x[1]); // cm
235 float L_wire = sqrt(Rad_wire*Rad_wire - dRho*dRho);
236 float flyLength = Lproj - L_wire;
237 if (x[1]<0) flyLength = Lproj + L_wire;
238 Tfly = flyLength/29.98*sqrt(1.0+(0.105/0.5)*(0.105/0.5))*fct; //approxiate... cm/ns
239#endif
240 float time = l->hit()->reccdc()->tdc - Tfly;
241
242 int wire = h.wire()->localId();
243 int layerId = h.wire()->layerId();
244 float dist = distance;
245 float edist = eDistance;
246 double T0 = l_mdcCalFunSvc->getT0(layerId,wire);
247 double rawadc = 0.;
248 rawadc =l->hit()->reccdc()->adc;
249
250 double timeWalk = l_mdcCalFunSvc->getTimeWalk(layerId, rawadc);
251 double drifttime =time -T0-timeWalk;
252 dist = l_mdcCalFunSvc->driftTimeToDist(drifttime,layerId, wire, side, 0.0);
253 edist = l_mdcCalFunSvc->getSigma(layerId, side, dist, 0.0);
254 dist = dist/10.0; //mm->cm
255 edist = edist/10.0;
256// if(layerId<20) edist = 9999.;
257
258//zsl calcdc_driftdist_(&doPropDelay,&wire,&side,p,x,&time,&dist,&edist);
259
260 distance = (double) dist;
261 eDistance = (double) edist;
262
263 l->drift(distance,0); //update
264 l->drift(distance,1);
265 l->dDrift(eDistance,0);
266 l->dDrift(eDistance,1);
267 l->tof(Tfly);
268 }
269 double eDistance2 = eDistance * eDistance;
270
271 //...Residual...
272 HepVector3D v = onTrack - onWire;
273 double vmag = v.mag();
274 double dDistance = vmag - distance;
275
276 //...dxda...
277 this->dxda(*l, t.helix(), dPhi, dxda, dyda, dzda, doSagCorrection);
278
279 //...Chi2 related...
280 // dDda = (v.x() * dxda + v.y() * dyda + v.z() * dzda) / vmag;
281 HepVector3D vw = h.wire()->direction();
282 dDda = (vmag > 0.)
283 ? ((v.x() * (1. - vw.x() * vw.x()) -
284 v.y() * vw.x() * vw.y() - v.z() * vw.x() * vw.z())
285 * dxda +
286 (v.y() * (1. - vw.y() * vw.y()) -
287 v.z() * vw.y() * vw.z() - v.x() * vw.y() * vw.x())
288 * dyda +
289 (v.z() * (1. - vw.z() * vw.z()) -
290 v.x() * vw.z() * vw.x() - v.y() * vw.z() * vw.y())
291 * dzda) / vmag
292 : Vector(5,0);
293 if(vmag<=0.0) {
294 std::cout << " in fit " << onTrack << ", " << onWire;
295 h.dump();
296 }
297 dchi2da += (dDistance / eDistance2) * dDda;
298 d2chi2d2a += vT_times_v(dDda) / eDistance2;
299 double pChi2 = dDistance * dDistance / eDistance2;
300 chi2 += pChi2;
301
302 //...Store results...
303 l->update(onTrack, onWire, leftRight, pChi2);
304 }
305
306 //...Check condition...
307 double change = chi2Old - chi2;
308 if (fabs(change) < convergence) break;
309 //if (change < 0.) break;
310 //jtanaka -- from traffs -- Ozaki-san added this part to traffs.
311 if (change < 0.){
312#ifdef TRKRECO_DEBUG_DETAIL
313 std::cout << "chi2Old, chi2=" << chi2Old <<" "<< chi2 << std::endl;
314#endif
315 //change to the old value.
316 a += factor*da;
317// a[2] = 0.000000001;
318 t._helix->a(a);
319
320 chi2 = 0.;
321 for (unsigned j = 0; j < 5; j++) dchi2da[j] = 0.;
322 d2chi2d2a = SymMatrix(5, 0);
323
324 //...Loop with hits...
325 unsigned i = 0;
326 while (TMLink * l = t.links()[i++]) {
327 const TMDCWireHit & h = * l->hit();
328
329 //...Check state...
330 if (h.state() & WireHitInvalidForFit) continue;
331 if (! (h.state() & WireHitFittingValid)) continue;
332
333 //...Check wire...
334 if (! nTrial)
335 if (h.wire()->stereo()) allAxial = false;
336
337 //...Cal. closest points...
338 int doSagCorrection = 0;
339// if( nTrial && !allAxial ) doSagCorrection = 1; //Liuqg
340 t.approach(* l, doSagCorrection);
341 double dPhi = l->dPhi();
342 const HepPoint3D & onTrack = l->positionOnTrack();
343 const HepPoint3D & onWire = l->positionOnWire();
344
345#ifdef TRKRECO_DEBUG_DETAIL
346 // std::cout << " in fit in case of change < 0. " << onTrack << ", " << onWire;
347 // h.dump();
348#endif
349
350 //...Obtain drift distance and its error...
351 unsigned leftRight = WireHitRight;
352 if (onWire.cross(onTrack).z() < 0.) leftRight = WireHitLeft;
353 double distance = h.drift(leftRight);
354 double eDistance = h.dDrift(leftRight);
355 if(nTrial && !allAxial){
356 int side = leftRight;
357
358 double Tfly = _t0_bes/220.*(110.-onWire.y());
359#if CAL_TOF_HELIX
360 float Rad_wire = sqrt(x[0]*x[0]+x[1]*x[1]); // cm
361 float L_wire = sqrt(Rad_wire*Rad_wire - dRho*dRho);
362 float flyLength = Lproj - L_wire;
363 if (x[1]<0) flyLength = Lproj + L_wire;
364 Tfly = flyLength/29.98*sqrt(1.0+(0.105/10)*(0.105/10))*fct; //approxiate... cm/ns
365#endif
366
367 float time = l->hit()->reccdc()->tdc - Tfly;
368
369 int wire = h.wire()->localId();
370 int layerId = h.wire()->layerId();
371 float dist= distance;
372 float edist = eDistance;
373// int doPropDelay = 1; //
374
375 double T0 = l_mdcCalFunSvc->getT0(layerId,wire);
376 double rawadc = 0.;
377 rawadc =l->hit()->reccdc()->adc;
378 double timeWalk = l_mdcCalFunSvc->getTimeWalk(layerId, rawadc);
379 double drifttime =time -T0-timeWalk;
380 dist = l_mdcCalFunSvc->driftTimeToDist(drifttime,layerId, wire, side, 0.0);
381 edist = l_mdcCalFunSvc->getSigma(layerId, side, dist, 0.0);
382 dist = dist/10.0; //mm->cm
383 edist = edist/10.0;
384// if (layerId<20) edist = 9999.;
385
386//zsl calcdc_driftdist_(&doPropDelay,&wire,&side,p,x,&time,&dist,&edist);
387
388 distance = (double) dist;
389 eDistance = (double) edist;
390
391 l->drift(distance,0); //update
392 l->drift(distance,1);
393 l->dDrift(eDistance,0);
394 l->dDrift(eDistance,1);
395 l->tof(Tfly);
396 }
397 double eDistance2 = eDistance * eDistance;
398
399 //...Residual...
400 HepVector3D v = onTrack - onWire;
401 double vmag = v.mag();
402 double dDistance = vmag - distance;
403
404 //...dxda...
405 this->dxda(*l, t.helix(), dPhi, dxda, dyda, dzda, doSagCorrection);
406
407 //...Chi2 related...
408 //dDda = (v.x() * dxda + v.y() * dyda + v.z() * dzda) / vmag;
409 HepVector3D vw = h.wire()->direction();
410 dDda = (vmag > 0.)
411 ? ((v.x() * (1. - vw.x() * vw.x()) -
412 v.y() * vw.x() * vw.y() - v.z() * vw.x() * vw.z())
413 * dxda +
414 (v.y() * (1. - vw.y() * vw.y()) -
415 v.z() * vw.y() * vw.z() - v.x() * vw.y() * vw.x())
416 * dyda +
417 (v.z() * (1. - vw.z() * vw.z()) -
418 v.x() * vw.z() * vw.x() - v.y() * vw.z() * vw.y())
419 * dzda) / vmag
420 : Vector(5,0);
421 if(vmag<=0.0) {
422 std::cout << " in fit " << onTrack << ", " << onWire;
423 h.dump();
424 }
425 dchi2da += (dDistance / eDistance2) * dDda;
426 d2chi2d2a += vT_times_v(dDda) / eDistance2;
427 double pChi2 = dDistance * dDistance / eDistance2;
428 chi2 += pChi2;
429
430 //...Store results...
431 l->update(onTrack, onWire, leftRight, pChi2);
432 }
433 //break;
434 factor *= 0.75;
435#ifdef TRKRECO_DEBUG_DETAIL
436 std::cout << "factor = " << factor << std::endl;
437 std::cout << "chi2 = " << chi2 << std::endl;
438#endif
439 if(factor < 0.01)break;
440 }
441
442 chi2Old = chi2;
443
444 //...Cal. helix parameters for next loop...
445 if (allAxial) {
446 f = dchi2da.sub(1, 3);
447 e = d2chi2d2a.sub(1, 3);
448 f = solve(e, f);
449 da[0] = f[0];
450 da[1] = f[1];
451 da[2] = f[2];
452 da[3] = 0.;
453 da[4] = 0.;
454 }
455 else {
456 da = solve(d2chi2d2a, dchi2da);
457 }
458
459#ifdef TRKRECO_DEBUG_DETAIL
460 //std::cout << " fit " << nTrial << " : " << da << std::endl;
461 //std::cout << " d2chi " << d2chi2d2a << std::endl;
462 //std::cout << " dchi2 " << dchi2da << std::endl;
463#endif
464
465 a -= factor*da;
466// a[2] = 0.000000001;
467 t._helix->a(a);
468 ++nTrial;
469 }
470
471 //...Cal. error matrix...
472 SymMatrix Ea(5, 0);
473 unsigned dim;
474 if (allAxial) {
475 dim = 3;
476 SymMatrix Eb(3, 0), Ec(3, 0);
477 Eb = d2chi2d2a.sub(1, 3);
478 Ec = Eb.inverse(err);
479 Ea[0][0] = Ec[0][0];
480 Ea[0][1] = Ec[0][1];
481 Ea[0][2] = Ec[0][2];
482 Ea[1][1] = Ec[1][1];
483 Ea[1][2] = Ec[1][2];
484 Ea[2][2] = Ec[2][2];
485 }
486 else {
487 dim = 5;
488 Ea = d2chi2d2a.inverse(err);
489 }
490
491 //...Store information...
492 if (! err) {
493 t._helix->a(a);
494 t._helix->Ea(Ea);
495 }
496 else {
497 err = TFitFailed;
498 }
499
500 t._ndf = nValid - dim;
501 t._chi2 = chi2;
502 // t._fitted = true;
503
504 return err;
505}
506
507/*
508
509// addition by matsu ( 1999/07/05 )
510int
511TCosmicFitter::fitWithCathode( TTrackBase &b, float t0Offset,
512 float windowSize , int SysCorr ) {
513
514#ifdef TRKRECO_DEBUG_DETAIL
515 std::cout << " TCosmicFitter::fitCathode ..." << std::endl;
516#endif
517
518 //...Already fitted ? ...
519 if (b.fittedWithCathode()) return 0;
520
521 //...Check # of his...
522 if( b.nCores() < 5 ) return -1;
523
524 //...Type check...
525 if (b.objectType() != Track) return TFitUnavailable;
526 TTrack & t = (TTrack &) b;
527
528 //...for cathode ndf...
529 int NusedCathode(0);
530
531 //...Setup...
532 unsigned nTrial = 0;
533 Vector a(5), da(5);
534 a = t.helix().a();
535 Vector dxda(5);
536 Vector dyda(5);
537 Vector dzda(5);
538 Vector dDda(5);
539 Vector dDzda(5); // for cathode z
540 Vector dchi2da(5);
541 SymMatrix d2chi2d2a(5, 0);
542 double chi2;
543 double chi2Old = DBL_MAX;
544 const double convergence = Convergence;
545 bool allAxial = true;
546 int LayerStat(0); // layer status axial:0 stereo:1 cathode:2
547 Matrix e(3, 3);
548 Vector f(3);
549 int err = 0;
550 double factor = 1.0;
551
552 const AList<TMDCCatHit> & chits = t.catHits();
553
554 Vector maxDouble(5);
555 for (unsigned i = 0; i < 5; i++) maxDouble[i] = (FLT_MAX);
556
557 //...Fitting loop...
558 while(nTrial < NTrailMax){
559
560 //...Set up...
561 chi2 = 0.;
562 NusedCathode = 0;
563 for (unsigned j = 0; j < 5; j++) dchi2da[j] = 0.;
564 d2chi2d2a = SymMatrix(5, 0);
565
566 //...Loop with hits...
567 unsigned i = 0;
568 AList<TMLink> cores = t.cores();
569 while (TMLink * l = cores[i++]) {
570
571 const TMDCWireHit & h = * l->hit();
572
573 // Check layer status ( cathode added )
574 LayerStat = 0;
575 if ( h.wire()->stereo() ) LayerStat = 1;
576 unsigned nlayer = h.wire()->layerId();
577 if (nlayer == 0 || nlayer == 1 || nlayer == 2 ) LayerStat = 2;
578
579 //...Check wire...
580 if (! nTrial)
581 if (h.wire()->stereo() || LayerStat == 2 ) allAxial = false;
582
583 //...Cal. closest points...
584 int doSagCorrection = 0;
585 if(nTrial && !allAxial ) doSagCorrection = 1;
586 t.approach(* l, doSagCorrection );
587 double dPhi = l->dPhi();
588 const HepPoint3D & onTrack = l->positionOnTrack();
589 const HepPoint3D & onWire = l->positionOnWire();
590
591#ifdef TRKRECO_DEBUG_DETAIL
592 // std::cout << " in fitCathode " << onTrack << ", " << onWire;
593 // h.dump();
594#endif
595
596 //...Obtain drift distance and its error...
597 unsigned leftRight = WireHitRight;
598 if (onWire.cross(onTrack).z() < 0.) leftRight = WireHitLeft;
599 double distance = h.drift(leftRight);
600 double eDistance = h.dDrift(leftRight);
601 //...
602 if(nTrial && !allAxial){
603 int wire = h.wire()->id();
604 int side = leftRight;
605 if(side==0) side = -1;
606 float x[3]={ onWire.x(), onWire.y(), onWire.z()};
607 float p[3]={t.p().x(), t.p().y(), t.p().z()};
608 float time = l->hit()->reccdc()->tdc + t0Offset;
609 float dist = distance;
610 float edist = eDistance;
611 int doPropDelay = 1;
612 calcdc_driftdist_(
613 &doPropDelay,&wire,&side,p,x,&time,&dist,&edist);
614 distance = (double) dist;
615 eDistance = (double) edist;
616 }
617 double eDistance2 = eDistance * eDistance;
618
619 //...Residual...
620 HepVector3D v = onTrack - onWire;
621 double vmag = v.mag();
622 double dDistance = vmag - distance;
623
624 //...dxda...
625 this->dxda( *l, t.helix(), dPhi, dxda, dyda, dzda , doSagCorrection );
626
627 // ... Chi2 related ...
628 double pChi2(0.);
629
630 if( LayerStat == 0 || LayerStat == 1 ){
631 // dDda = (v.x() * dxda + v.y() * dyda + v.z() * dzda) / vmag;
632 Vector3 vw = h.wire()->direction();
633 dDda = (vmag > 0.)
634 ? ((v.x() * (1. - vw.x() * vw.x()) -
635 v.y() * vw.x() * vw.y() - v.z() * vw.x() * vw.z())
636 * dxda +
637 (v.y() * (1. - vw.y() * vw.y()) -
638 v.z() * vw.y() * vw.z() - v.x() * vw.y() * vw.x())
639 * dyda +
640 (v.z() * (1. - vw.z() * vw.z()) -
641 v.x() * vw.z() * vw.x() - v.y() * vw.z() * vw.y())
642 * dzda) / vmag
643 : Vector(5,0);
644 if(vmag<=0.0) {
645 std::cout << " in fit " << onTrack << ", " << onWire;
646 h.dump();
647 }
648 dchi2da += (dDistance / eDistance2) * dDda;
649 d2chi2d2a += vT_times_v(dDda) / eDistance2;
650 double pChi2 = dDistance * dDistance / eDistance2;
651 chi2 += pChi2;
652
653 } else {
654
655
656 dDda = ( v.x() * dxda + v.y() * dyda )/v.perp();
657 dchi2da += (dDistance/eDistance2) * dDda;
658 d2chi2d2a += vT_times_v(dDda)/eDistance2;
659
660 pChi2 = 0;
661 pChi2 += (dDistance/eDistance) * (dDistance/eDistance) ;
662
663 if ( l->usecathode() >= 3 ) {
664
665 TMDCClust * mclust = l->getmclust();
666
667 double dDistanceZ(t.helix().x(dPhi).z());
668
669 if( SysCorr ){
670 if( !nTrial ) {
671 mclust->zcalc( atan( t.helix().tanl()) );
672 mclust->esterz( atan( t.helix().tanl()) );
673 }
674
675 dDistanceZ -= mclust->zclustWithSysCorr();
676 } else {
677 dDistanceZ -= mclust->zclust();
678 }
679
680 double eDistanceZ(mclust->erz());
681 if( !eDistanceZ ) eDistanceZ = 99999.;
682
683
684 double eDistance2Z = eDistanceZ * eDistanceZ;
685 double pChi2z = 0;
686 pChi2z = (dDistanceZ/eDistanceZ);
687 pChi2z *= pChi2z;
688
689#ifdef TRKRECO_DEBUG_DETAIL
690 std::cout << "dDistanceZ = " << dDistanceZ << std::endl;
691#endif
692
693 //.... requirement for use of cluster
694 if( nTrial == 0 &&
695 fabs(dDistanceZ)< windowSize &&
696 mclust->chits().length() == 1
697 ) l->setusecathode(4 );
698
699 if( l->usecathode() == 4 ){
700 NusedCathode++;
701 dDzda = dzda ;
702 dchi2da += (dDistanceZ/eDistance2Z) * dDzda;
703 d2chi2d2a += vT_times_v(dDzda)/eDistance2Z;
704 pChi2 += pChi2z ;
705
706 }
707 }
708
709 } // end Chi2 related
710
711 chi2 += pChi2;
712
713
714 //...Store results...
715 l->update(onTrack, onWire, leftRight, pChi2);
716
717
718 } // TMLink loop end
719
720 //...Check condition...
721 double change = chi2Old - chi2;
722 //if(TRKRECO_DEBUG_DETAIL>0) std::cout << "chi2 change = " <<change << std::endl;
723
724 if (fabs(change) < convergence) break;
725 if (change < 0.) {
726
727#ifdef TRKRECO_DEBUG_DETAIL
728 std::cout << "chi2Old, chi2=" << chi2Old <<" "<< chi2 << std::endl;
729#endif
730
731 NusedCathode = 0;
732 //change to the old value.
733 a += factor*da;
734 t._helix->a(a);
735
736
737 chi2 = 0.;
738 for (unsigned j = 0; j < 5; j++) dchi2da[j] = 0.;
739 d2chi2d2a = SymMatrix(5, 0);
740
741
742 //...Loop with hits...
743 unsigned i = 0;
744 while (TMLink * l = cores[i++]) {
745
746 const TMDCWireHit & h = * l->hit();
747
748 // Check layer status ( cathode added )
749 LayerStat = 0;
750 if ( h.wire()->stereo() ) LayerStat = 1;
751 unsigned nlayer = h.wire()->layerId();
752 if (nlayer == 0 || nlayer == 1 || nlayer == 2 ) LayerStat = 2;
753
754 //...Cal. closest points...
755 int doSagCorrection = 0;
756 if( nTrial && !allAxial ) doSagCorrection = 1;
757 t.approach(* l , doSagCorrection );
758 double dPhi = l->dPhi();
759 const HepPoint3D & onTrack = l->positionOnTrack();
760 const HepPoint3D & onWire = l->positionOnWire();
761
762#ifdef TRKRECO_DEBUG_DETAIL
763 // std::cout << " in fitCathode " << onTrack << ", " << onWire;
764 // h.dump();
765#endif
766
767 //...Obtain drift distance and its error...
768 unsigned leftRight = WireHitRight;
769 if (onWire.cross(onTrack).z() < 0.) leftRight = WireHitLeft;
770 double distance = h.drift(leftRight);
771 double eDistance = h.dDrift(leftRight);
772 if(nTrial && !allAxial){
773 int wire = h.wire()->id();
774 int side = leftRight;
775 if(side==0) side = -1;
776 float x[3]={ onWire.x(), onWire.y(), onWire.z()};
777 float p[3]={t.p().x(), t.p().y(), t.p().z()};
778 float time = l->hit()->reccdc()->tdc + t0Offset;
779 float dist= distance;
780 float edist = eDistance;
781 int doPropDelay = 1; //
782 calcdc_driftdist_(
783 &doPropDelay,&wire,&side,p,x,&time,&dist,&edist);
784 distance = (double) dist;
785 eDistance = (double) edist;
786 }
787 double eDistance2 = eDistance * eDistance;
788
789 //...Residual...
790 HepVector3D v = onTrack - onWire;
791 double vmag = v.mag();
792 double dDistance = vmag - distance;
793
794 //...dxda...
795 this->dxda(* l, t.helix(), dPhi, dxda, dyda, dzda , doSagCorrection );
796
797 // ... Chi2 related ...
798 double pChi2(0.);
799
800
801 if( LayerStat == 0 || LayerStat == 1 ){
802 // dDda = (v.x() * dxda + v.y() * dyda + v.z() * dzda) / vmag;
803 Vector3 vw = h.wire()->direction();
804 dDda = (vmag > 0.)
805 ? ((v.x() * (1. - vw.x() * vw.x()) -
806 v.y() * vw.x() * vw.y() - v.z() * vw.x() * vw.z())
807 * dxda +
808 (v.y() * (1. - vw.y() * vw.y()) -
809 v.z() * vw.y() * vw.z() - v.x() * vw.y() * vw.x())
810 * dyda +
811 (v.z() * (1. - vw.z() * vw.z()) -
812 v.x() * vw.z() * vw.x() - v.y() * vw.z() * vw.y())
813 * dzda) / vmag
814 : Vector(5,0);
815 if(vmag<=0.0) {
816 std::cout << " in fit " << onTrack << ", " << onWire;
817 h.dump();
818 }
819 dchi2da += (dDistance / eDistance2) * dDda;
820 d2chi2d2a += vT_times_v(dDda) / eDistance2;
821 double pChi2 = dDistance * dDistance / eDistance2;
822 chi2 += pChi2;
823
824 } else {
825
826 dDda = ( v.x() * dxda + v.y() * dyda )/v.perp();
827 dchi2da += (dDistance/eDistance2) * dDda;
828 d2chi2d2a += vT_times_v(dDda)/eDistance2;
829
830 pChi2 = 0;
831 pChi2 += (dDistance/eDistance) * (dDistance/eDistance) ;
832
833 if( l->usecathode() == 4 ){
834
835 TMDCClust * mclust = l->getmclust();
836
837 if( mclust ){
838 NusedCathode++;
839
840 double dDistanceZ(t.helix().x(dPhi).z());
841 if( SysCorr ) dDistanceZ -= mclust->zclustWithSysCorr();
842 else dDistanceZ -= mclust->zclust();
843
844 double eDistanceZ(99999.);
845 if( mclust->erz() != 0. ) eDistanceZ = mclust->erz();
846
847 double eDistance2Z = eDistanceZ * eDistanceZ;
848 double pChi2z = 0;
849 pChi2z = (dDistanceZ/eDistanceZ);
850 pChi2z *= pChi2z;
851
852 dDzda = dzda ;
853 dchi2da += (dDistanceZ/eDistance2Z) * dDzda;
854 d2chi2d2a += vT_times_v(dDzda)/eDistance2Z;
855 pChi2 += pChi2z ;
856 }
857
858 }
859
860 } // end Chi2 related
861
862 chi2 += pChi2;
863
864
865 //...Store results...
866 l->update(onTrack, onWire, leftRight, pChi2);
867
868 }
869
870 //break;
871
872 factor *= 0.75;
873#ifdef TRKRECO_DEBUG_DETAIL
874 std::cout << "factor = " << factor << std::endl;
875 std::cout << "chi2 = " << chi2 << std::endl;
876#endif
877 if(factor < 0.01)break;
878
879 }
880
881 chi2Old = chi2;
882
883 //...Cal. helix parameters for next loop...
884 if (allAxial) {
885 f = dchi2da.sub(1, 3);
886 e = d2chi2d2a.sub(1, 3);
887 f = solve(e, f);
888 da[0] = f[0];
889 da[1] = f[1];
890 da[2] = f[2];
891 da[3] = 0.;
892 da[4] = 0.;
893 }
894 else {
895 da = solve(d2chi2d2a, dchi2da);
896 }
897
898#ifdef TRKRECO_DEBUG_DETAIL
899 //std::cout << " fit " << nTrial << " : " << da << std::endl;
900 //std::cout << " d2chi " << d2chi2d2a << std::endl;
901 //std::cout << " dchi2 " << dchi2da << std::endl;
902#endif
903
904 a -= da;
905 t._helix->a(a);
906 ++nTrial;
907 }
908
909
910 //...Cal. error matrix...
911 SymMatrix Ea(5, 0);
912 unsigned dim;
913 if (allAxial) {
914 dim = 3;
915 SymMatrix Eb(3, 0), Ec(3, 0);
916 Eb = d2chi2d2a.sub(1, 3);
917 Ec = Eb.inverse(err);
918 Ea[0][0] = Ec[0][0];
919 Ea[0][1] = Ec[0][1];
920 Ea[0][2] = Ec[0][2];
921 Ea[1][1] = Ec[1][1];
922 Ea[1][2] = Ec[1][2];
923 Ea[2][2] = Ec[2][2];
924 }
925 else {
926 dim = 5;
927 Ea = d2chi2d2a.inverse(err);
928 }
929
930 //...Store information...
931 if (! err) {
932 t._helix->a(a);
933 t._helix->Ea(Ea);
934 } else {
935 err = TFitFailed;
936 }
937
938 t._ndf = t.nCores() + NusedCathode - dim;
939
940 t._chi2 = chi2;
941
942 t._fittedWithCathode = true;
943
944 return err;
945}
946
947// end of addition
948
949*/
950
951
952int
953TCosmicFitter::dxda(const TMLink & link,
954 const Helix & h,
955 double dPhi,
956 Vector & dxda,
957 Vector & dyda,
958 Vector & dzda,
959 int doSagCorrection) const {
960
961 //...Setup...
962 const TMDCWire & w = * link.wire();
963 Vector a = h.a();
964 double dRho = a[0];
965 double phi0 = a[1];
966 double kappa = a[2];
967 // double rho = Helix::ConstantAlpha / kappa;
968 // double rho = 333.564095 / kappa;
969 double rho = 333.564095/(-1000*(m_pmgnIMF->getReferField())) / kappa;
970cout<<"TCosmicFitter::rho: "<<333.564095/(-1000*(m_pmgnIMF->getReferField()))<<endl;
971 double tanLambda = a[4];
972
973 double sinPhi0 = sin(phi0);
974 double cosPhi0 = cos(phi0);
975 double sinPhi0dPhi = sin(phi0 + dPhi);
976 double cosPhi0dPhi = cos(phi0 + dPhi);
977 Vector dphida(5);
978 //... sag correction
979 HepPoint3D xw = w.xyPosition();
980 HepPoint3D wireBackwardPosition = w.backwardPosition();
981 HepVector3D v = w.direction();
982
983 if(doSagCorrection ){
984 cout<<"doSagCorrection in TCosmicFitter!"<<endl;
985 int wireID = w.id();
986 float wirePosition[3];
987 float wireZ = link.positionOnTrack().z();
988 float dydz =0;
989 float ybSag = 0;
990 float yfSag = 0;
991 if(wireZ >w.backwardPosition().z() &&
992 wireZ < w.forwardPosition().z() ){
993
994//zsl calcdc_sag3_(&wireID, &wireZ, wirePosition, &dydz, &ybSag, &yfSag);
995
996 //... wire Position
997 xw.setX( (double) wirePosition[0]);
998 xw.setY( (double) wirePosition[1]);
999 xw.setZ( (double) wirePosition[2]);
1000 //...
1001 wireBackwardPosition.setY((double) ybSag);
1002 //...
1003 //... v.setY((double) dydz);
1004 HepVector3D v_aux(w.forwardPosition().x() - w.backwardPosition().x(),
1005 yfSag - ybSag,
1006 w.forwardPosition().z() - w.backwardPosition().z());
1007 v = v_aux.unit();
1008 }
1009 }
1010
1011 //...Axial case...
1012 if (w.axial()) {
1013 // HepPoint3D d = h.center() - w.xyPosition();
1014 HepPoint3D d = h.center() - xw;
1015 double dmag2 = d.mag2();
1016
1017 dphida[0] = (sinPhi0 * d.x() - cosPhi0 * d.y()) / dmag2;
1018 dphida[1] = (dRho + rho) * (cosPhi0 * d.x() + sinPhi0 * d.y())
1019 / dmag2 - 1.;
1020 dphida[2] = (- rho / kappa) * (sinPhi0 * d.x() - cosPhi0 * d.y())
1021 / dmag2;
1022 dphida[3] = 0.;
1023 dphida[4] = 0.;
1024 }
1025
1026 //...Stereo case...
1027 else {
1028 HepPoint3D onTrack = h.x(dPhi);
1029 Vector c(3);
1030 //c = HepPoint3D(w.backwardPosition() - (v * w.backwardPosition()) * v);
1031 c = HepPoint3D(wireBackwardPosition - (v * wireBackwardPosition) * v);
1032
1033 Vector x(3);
1034 x = onTrack;
1035
1036 Vector dxdphi(3);
1037 dxdphi[0] = rho * sinPhi0dPhi;
1038 dxdphi[1] = - rho * cosPhi0dPhi;
1039 dxdphi[2] = - rho * tanLambda;
1040
1041 Vector d2xdphi2(3);
1042 d2xdphi2[0] = rho * cosPhi0dPhi;
1043 d2xdphi2[1] = rho * sinPhi0dPhi;
1044 d2xdphi2[2] = 0.;
1045
1046 double dxdphi_dot_v = dxdphi[0]*v.x() + dxdphi[1]*v.y() + dxdphi[2]*v.z();
1047 double x_dot_v = x[0]*v.x() + x[1]*v.y() + x[2]*v.z();
1048
1049 double dfdphi = - (dxdphi[0] - dxdphi_dot_v*v.x()) * dxdphi[0]
1050 - (dxdphi[1] - dxdphi_dot_v*v.y()) * dxdphi[1]
1051 - (dxdphi[2] - dxdphi_dot_v*v.z()) * dxdphi[2]
1052 - (x[0] - c[0] - x_dot_v*v.x()) * d2xdphi2[0]
1053 - (x[1] - c[1] - x_dot_v*v.y()) * d2xdphi2[1];
1054 /* - (x[2] - c[2] - x_dot_v*v.z()) * d2xdphi2[2]; = 0. */
1055
1056
1057 //dxda_phi, dyda_phi, dzda_phi : phi is fixed
1058 Vector dxda_phi(5);
1059 dxda_phi[0] = cosPhi0;
1060 dxda_phi[1] = -(dRho + rho)*sinPhi0 + rho*sinPhi0dPhi;
1061 dxda_phi[2] = -(rho/kappa)*( cosPhi0 - cosPhi0dPhi );
1062 dxda_phi[3] = 0.;
1063 dxda_phi[4] = 0.;
1064
1065 Vector dyda_phi(5);
1066 dyda_phi[0] = sinPhi0;
1067 dyda_phi[1] = (dRho + rho)*cosPhi0 - rho*cosPhi0dPhi;
1068 dyda_phi[2] = -(rho/kappa)*( sinPhi0 - sinPhi0dPhi );
1069 dyda_phi[3] = 0.;
1070 dyda_phi[4] = 0.;
1071
1072 Vector dzda_phi(5);
1073 dzda_phi[0] = 0.;
1074 dzda_phi[1] = 0.;
1075 dzda_phi[2] = (rho/kappa)*tanLambda*dPhi;
1076 dzda_phi[3] = 1.;
1077 dzda_phi[4] = -rho*dPhi;
1078
1079 Vector d2xdphida(5);
1080 d2xdphida[0] = 0.;
1081 d2xdphida[1] = rho*cosPhi0dPhi;
1082 d2xdphida[2] = -(rho/kappa)*sinPhi0dPhi;
1083 d2xdphida[3] = 0.;
1084 d2xdphida[4] = 0.;
1085
1086 Vector d2ydphida(5);
1087 d2ydphida[0] = 0.;
1088 d2ydphida[1] = rho*sinPhi0dPhi;
1089 d2ydphida[2] = (rho/kappa)*cosPhi0dPhi;
1090 d2ydphida[3] = 0.;
1091 d2ydphida[4] = 0.;
1092
1093 Vector d2zdphida(5);
1094 d2zdphida[0] = 0.;
1095 d2zdphida[1] = 0.;
1096 d2zdphida[2] = (rho/kappa)*tanLambda;
1097 d2zdphida[3] = 0.;
1098 d2zdphida[4] = -rho;
1099
1100 Vector dfda(5);
1101 for( int i = 0; i < 5; i++ ){
1102 double d_dot_v = v.x()*dxda_phi[i]
1103 + v.y()*dyda_phi[i]
1104 + v.z()*dzda_phi[i];
1105 dfda[i] = - (dxda_phi[i] - d_dot_v*v.x()) * dxdphi[0]
1106 - (dyda_phi[i] - d_dot_v*v.y()) * dxdphi[1]
1107 - (dzda_phi[i] - d_dot_v*v.z()) * dxdphi[2]
1108 - (x[0] - c[0] - x_dot_v*v.x()) * d2xdphida[i]
1109 - (x[1] - c[1] - x_dot_v*v.y()) * d2ydphida[i]
1110 - (x[2] - c[2] - x_dot_v*v.z()) * d2zdphida[i];
1111 dphida[i] = - dfda[i] /dfdphi;
1112 }
1113 }
1114
1115 dxda[0] = cosPhi0 + rho * sinPhi0dPhi * dphida[0];
1116 dxda[1] = - (dRho + rho) * sinPhi0 + rho * sinPhi0dPhi * (1. + dphida[1]);
1117 dxda[2] = - rho / kappa * (cosPhi0 - cosPhi0dPhi)
1118 + rho * sinPhi0dPhi * dphida[2];
1119 dxda[3] = rho * sinPhi0dPhi * dphida[3];
1120 dxda[4] = rho * sinPhi0dPhi * dphida[4];
1121
1122 dyda[0] = sinPhi0 - rho * cosPhi0dPhi * dphida[0];
1123 dyda[1] = (dRho + rho) * cosPhi0 - rho * cosPhi0dPhi * (1. + dphida[1]);
1124 dyda[2] = - rho / kappa * (sinPhi0 - sinPhi0dPhi)
1125 - rho * cosPhi0dPhi * dphida[2];
1126 dyda[3] = - rho * cosPhi0dPhi * dphida[3];
1127 dyda[4] = - rho * cosPhi0dPhi * dphida[4];
1128
1129 dzda[0] = - rho * tanLambda * dphida[0];
1130 dzda[1] = - rho * tanLambda * dphida[1];
1131 dzda[2] = rho / kappa * tanLambda * dPhi - rho * tanLambda * dphida[2];
1132 dzda[3] = 1. - rho * tanLambda * dphida[3];
1133 dzda[4] = - rho * dPhi - rho * tanLambda * dphida[4];
1134
1135 return 0;
1136}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
w
double w
Definition: EvtPhokharaDef.hh:26
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: Gam4pikp.cxx:37
DBL_MAX
#define DBL_MAX
Definition: KalFitAlg.h:13
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition: KarLud.h:35
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
Convergence
#define Convergence
Definition: TCosmicFitter.cxx:76
Matrix
CLHEP::HepMatrix Matrix
Definition: TCosmicFitter.cxx:71
NTrailMax
#define NTrailMax
Definition: TCosmicFitter.cxx:75
WireHitFittingValid
#define WireHitFittingValid
Definition: TMDCWireHit.h:29
WireHitLeft
#define WireHitLeft
Definition: TMDCWireHit.h:21
WireHitRight
#define WireHitRight
Definition: TMDCWireHit.h:22
WireHitInvalidForFit
#define WireHitInvalidForFit
Definition: TMDCWireHit.h:55
TFitAlreadyFitted
#define TFitAlreadyFitted
Definition: TMFitter.h:28
TFitFailed
#define TFitFailed
Definition: TMFitter.h:30
TFitUnavailable
#define TFitUnavailable
Definition: TMFitter.h:31
TFitErrorFewHits
#define TFitErrorFewHits
Definition: TMFitter.h:29
Track
#define Track
Definition: TTrackBase.h:30
t
TTree * t
Definition: binning.cxx:23
Helix::center
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:203
Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/src/Helix.cxx:129
Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:263
IMagneticFieldSvc::getReferField
virtual double getReferField()=0
IMdcCalibFunSvc::driftTimeToDist
virtual double driftTimeToDist(double drifttime, int layid, int cellid, int lr, double entrance=0.0) const =0
IMdcCalibFunSvc::getSigma
virtual double getSigma(int layid, int lr, double dist, double entrance=0.0, double tanlam=0.0, double z=0.0, double Q=1000.0) const =0
IMdcCalibFunSvc::getT0
virtual double getT0(int layid, int cellid) const =0
IMdcCalibFunSvc::getTimeWalk
virtual double getTimeWalk(int layid, double Q) const =0
TCosmicFitter::~TCosmicFitter
virtual ~TCosmicFitter()
Destructor.
Definition: TCosmicFitter.cxx:93
TCosmicFitter::TCosmicFitter
TCosmicFitter(const std::string &name)
Constructor.
Definition: TCosmicFitter.cxx:84
TCosmicFitter::fit
int fit(TTrackBase &) const
Definition: TCosmicFitter.cxx:97
TMDCWireHit::dump
void dump(const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
dumps debug information.
Definition: TMDCWireHit.cxx:64
TMDCWireHit::state
unsigned state(void) const
returns state.
Definition: TMDCWireHit.h:230
TMDCWireHit::dDrift
float dDrift(unsigned) const
returns drift distance error.
Definition: TMDCWireHit.h:243
TMDCWireHit::drift
float drift(unsigned) const
returns drift distance.
Definition: TMDCWireHit.h:236
TMDCWireHit::wire
const TMDCWire *const wire(void) const
returns a pointer to a TMDCWire.
Definition: TMDCWireHit.h:218
TMDCWire
A class to represent a wire in MDC.
Definition: TMDCWire.h:55
TMDCWire::localId
unsigned localId(void) const
returns local id in a wire layer.
Definition: TMDCWire.h:213
TMDCWire::layerId
unsigned layerId(void) const
returns layer id.
Definition: TMDCWire.h:219
TMDCWire::direction
const HepVector3D & direction(void) const
returns direction vector of the wire.
Definition: TMDCWire.h:342
TMDCWire::stereo
bool stereo(void) const
returns true if this wire is in a stereo layer.
Definition: TMDCWire.h:354
TMFitter
A class to fit a TTrackBase object.
Definition: TMFitter.h:34
TTrackBase
A virtual class for a track class in tracking.
Definition: TTrackBase.h:46
TTrack
A class to represent a track in tracking.
Definition: TTrack.h:129
x
double x[1000]
Definition: draw_charge_space_coarse_II.cxx:8
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
b
const double b
Definition: slope.cxx:9