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