BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
TTrack.cxx
Go to the documentation of this file.
1#define OPTNK
2//-----------------------------------------------------------------------------
3// $Id: TTrack.cxx,v 1.28 2022/01/28 22:14:13 maqm Exp $
4//-----------------------------------------------------------------------------
5// Filename : TTrack.h
6// Section : Tracking
7// Owner : Yoshi Iwasaki
8// Email : yoshihito.iwasaki@kek.jp
9//-----------------------------------------------------------------------------
10// Description : A class to represent a track in tracking.
11// See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
12//-----------------------------------------------------------------------------
13
14/* for isnan */
15#if defined( __sparc )
16# if defined( __EXTENSIONS__ )
17# include <cmath>
18# else
19# define __EXTENSIONS__
20# include <cmath>
21# undef __EXTENSIONS__
22# endif
23#elif defined( __GNUC__ )
24# if defined( _XOPEN_SOURCE )
25# include <cmath>
26# else
27# define _XOPEN_SOURCE
28# include <cmath>
29# undef _XOPEN_SOURCE
30# endif
31#endif
32#include <cfloat>
33// #define HEP_SHORT_NAMES
34// #include "CLHEP/String/Strings.h"
35#include "CLHEP/Matrix/DiagMatrix.h"
36#include "CLHEP/Matrix/Matrix.h"
37#include "CLHEP/Matrix/SymMatrix.h"
38#include "CLHEP/Matrix/Vector.h"
39#include "TrkReco/TCircle.h"
40#include "TrkReco/TMDCUtil.h"
41#include "TrkReco/TMDCWire.h"
42#include "TrkReco/TMDCWireHit.h"
43#include "TrkReco/TMDCWireHitMC.h"
44#include "TrkReco/TMLine.h"
45#include "TrkReco/TMLink.h"
46// #include "TrkReco/TMDCStrip.h"
47#include "TrkReco/TSegment.h"
48#include "TrkReco/TTrack.h"
49// #include "tables/cdc.h"
50// #include "tables/trk.h"
51// #include "tables/mdst.h"
52// #include "tables/hepevt.h"
53#include "MdcTables/HepevtTables.h"
54#include "MdcTables/MdcTables.h"
55#include "MdcTables/MdstTables.h"
56#include "MdcTables/TrkTables.h"
57
58#include "CLHEP/Matrix/Matrix.h"
59#include "GaudiKernel/Bootstrap.h"
60#include "GaudiKernel/IDataProviderSvc.h"
61#include "GaudiKernel/IInterface.h"
62#include "GaudiKernel/IMessageSvc.h"
63#include "GaudiKernel/ISvcLocator.h"
64#include "GaudiKernel/Kernel.h"
65#include "GaudiKernel/MsgStream.h"
66#include "GaudiKernel/Service.h"
67#include "GaudiKernel/SmartDataPtr.h"
68#include "GaudiKernel/StatusCode.h"
69
70using CLHEP::HepDiagMatrix;
71using CLHEP::HepMatrix;
72using CLHEP::HepSymMatrix;
73using CLHEP::HepVector;
74
75const THelixFitter TTrack::_fitter = THelixFitter( "TTrack Default Helix Fitter" );
76
77TTrack::TTrack( const TCircle& c )
78 : TTrackBase( c.links() )
79 , _helix( new Helix( ORIGIN, Vector( 5, 0 ), SymMatrix( 5, 0 ) ) )
80 , _charge( c.charge() )
81 , _ndf( 0 )
82 , _chi2( 0. )
83 , _name( "none" )
84 , _type( 0 )
85 , _state( 0 )
86 , _mother( 0 )
87 , _daughter( 0 )
88 , m_pmgnIMF( nullptr ) {
89
90 //_fitter.setMagneticFieldPointer(m_pmgnIMF);//yzhang add 2012-05-04
91 // cout<<"TTrack: "<<getPmgnIMF()->getReferField()<<endl;
92 //...Set a defualt fitter...
93 fitter( &TTrack::_fitter );
94
95 //...Calculate helix parameters...
96 Vector a( 5 );
97 a[1] = fmod( atan2( _charge * ( c.center().y() - ORIGIN.y() ),
98 _charge * ( c.center().x() - ORIGIN.x() ) ) +
99 4. * M_PI,
100 2. * M_PI );
101 // a[2] = Helix::ConstantAlpha / c.radius();
102 // a[2] = 333.564095 / c.radius();
103 const double Bz = -1000 * getPmgnIMF()->getReferField();
104 a[2] = 333.564095 / ( c.radius() * Bz );
105 a[0] = ( c.center().x() - ORIGIN.x() ) / cos( a[1] ) - c.radius();
106 a[3] = 0.;
107 a[4] = 0.;
108 _helix->a( a );
109
110 //...Update links...
111 unsigned n = _links.length();
112 for ( unsigned i = 0; i < n; i++ ) _links[i]->track( this );
113
114 _fitted = false;
115 _fittedWithCathode = false;
116 /*
117 //jialk
118 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
119 if(scmgn!=StatusCode::SUCCESS) {
120 std::cout<< __FILE__<<" "<<__LINE__<<" Unable to open Magnetic field
121 service"<<std::endl;
122 }
123 */
124}
125
126TTrack::TTrack( const TTrack& a )
127 : TTrackBase( (TTrackBase&)a )
128 , _charge( a._charge )
129 , _segments( a._segments )
130 , _helix( new Helix( *a._helix ) )
131 , _ndf( a._ndf )
132 , _chi2( a._chi2 )
133 , _name( "copy of " + a._name )
134 , _type( a._type )
135 , _state( a._state )
136 , _mother( a._mother )
137 , _daughter( a._daughter )
138 , m_pmgnIMF( nullptr ) {}
139
140/*TTrack::TTrack(const T3DLine & a)
141: TTrackBase((TTrackBase &) a),
142 _charge(0),
143 _helix(new Helix(* a.helix())),
144 _ndf(a.ndf()),
145 _chi2(a.chi2()),
146 _name(0),
147 _type(a.type()),
148 _state(1),
149 _mother(0),
150 _daughter(0) {
151}*/
152
153TTrack::TTrack( const TRunge& a )
154 : TTrackBase( (TTrackBase&)a )
155 , _helix( new Helix( a.helix() ) )
156 , _ndf( a.ndf() )
157 , _chi2( a.chi2() )
158 , _name( "no" )
159 , _type( 0 )
160 , _state( 0 )
161 , _mother( 0 )
162 , _daughter( 0 )
163 , m_pmgnIMF( nullptr ) {
164 if ( _helix->kappa() > 0. ) _charge = 1.;
165 else _charge = -1.;
166}
167
168TTrack::TTrack( const Helix& h )
169 : TTrackBase()
170 , _helix( new Helix( h ) )
171 , _ndf( 0 )
172 , _chi2( 0. )
173 , _name( "none" )
174 , _type( 0 )
175 , _state( 0 )
176 , _mother( 0 )
177 , _daughter( 0 )
178 , m_pmgnIMF( nullptr ) {
179
180 //...Set a defualt fitter...
181 fitter( &TTrack::_fitter );
182
183 if ( _helix->kappa() > 0. ) _charge = 1.;
184 else _charge = -1.;
185
186 _fitted = false;
187 _fittedWithCathode = false;
188}
189
190TTrack::TTrack()
191 : TTrackBase()
192 , _charge( 1. )
193 , _helix( new Helix( ORIGIN, Vector( 5, 0 ), SymMatrix( 5, 0 ) ) )
194 , _ndf( 0 )
195 , _chi2( 0. )
196 , _name( "empty track" )
197 , _state( 0 )
198 , _mother( 0 )
199 , _daughter( 0 )
200 , m_pmgnIMF( nullptr ) {}
201
202TTrack::~TTrack() { delete _helix; }
203
204void TTrack::dump( const std::string& msg, const std::string& pre0 ) const {
205 bool def = false;
206 if ( msg == "" ) def = true;
207 std::string pre = pre0;
208 std::string tab;
209 for ( unsigned i = 0; i < pre.length(); i++ ) tab += " ";
210 if ( def || msg.find( "track" ) != std::string::npos ||
211 msg.find( "detail" ) != std::string::npos )
212 {
213 std::cout << pre << p() << "(pt=" << pt() << ")" << std::endl
214 << tab << "links=" << _links.length() << "(cores=" << nCores()
215 << "),chrg=" << _charge << ",ndf=" << _ndf << ",chi2=" << _chi2 << std::endl;
216 pre = tab;
217 }
218 if ( msg.find( "helix" ) != std::string::npos || msg.find( "detail" ) != std::string::npos )
219 {
220 std::cout << pre << "pivot=" << _helix->pivot() << ",center=" << _helix->center()
221 << std::endl
222 << pre << "helix=(" << _helix->a()[0] << "," << _helix->a()[1] << ","
223 << _helix->a()[2] << "," << _helix->a()[3] << "," << _helix->a()[4] << ")"
224 << std::endl;
225 pre = tab;
226 }
227 if ( !def ) TTrackBase::dump( msg, pre );
228}
229
230void TTrack::movePivot( void ) {
231 unsigned n = _links.length();
232 if ( !n )
233 {
234 std::cout << "TTrack::movePivot !!! can't move a pivot"
235 << " because of no link";
236 std::cout << std::endl;
237 return;
238 }
239
240 //...Check cores...
241 const AList<TMLink>& cores = TTrackBase::cores();
242 const AList<TMLink>* links = &cores;
243 if ( cores.length() == 0 ) links = &_links;
244
245 //...Hit loop...
246 unsigned innerMost = 0;
247 unsigned innerMostLayer = ( *links )[0]->wire()->layerId();
248 n = links->length();
249 for ( unsigned i = 1; i < n; i++ )
250 {
251 TMLink* l = ( *links )[i];
252 if ( l->wire()->layerId() < innerMostLayer )
253 {
254 innerMost = i;
255 innerMostLayer = l->wire()->layerId();
256 }
257 }
258
259 //...Move pivot...
260 HepPoint3D newPivot = ( *links )[innerMost]->positionOnWire();
261 if ( quality() & TrackQuality2D ) newPivot.setZ( 0. );
262 _helix->pivot( newPivot );
263}
264
265int TTrack::approach( TMLink& l ) const { return approach( l, false ); }
266
267void TTrack::refine2D( AList<TMLink>& list, float maxSigma ) {
268 unsigned n = _links.length();
269 AList<TMLink> bad;
270 for ( unsigned i = 0; i < n; ++i )
271 {
272 if ( _links[i]->pull() > maxSigma ) bad.append( _links[i] );
273 }
274 _links.remove( bad );
275 if ( bad.length() )
276 {
277 _fitted = false;
278 fit2D();
279 }
280 list.append( bad );
281}
282
283int TTrack::HelCyl( double rhole, double rCyl, double zb, double zf, double epsl, double& phi,
284 HepPoint3D& xp ) const {
285
286 int status( 0 ); // return value
287 //---------------------------------------------------------------------
288 // value | ext | status
289 //---------------------------------------------------------------------
290 // 1. | OK |
291 // -1. | NO | charge = 0
292 // 0. | NO | | tanl | < 0.1 ( neglect | lamda | < 5.7 deg. )
293 // | | or | dPhi | > 2 pi ( neglect curly track )
294 // | | or cannot reach to r=rhole at z = zb or zf.
295 // 2. | OK | backward , ext point set on z = zb
296 // 3. | OK | forward , ext point set on z = zf
297 //---------------------------------------------------------------------
298 // * when value = 0,2,3 , ext(z) <= zb or ext(z) >= zf
299
300 //--- debug
301 // std::cout << " " << std::endl;
302 // std::cout << "HelCyl called .. rhole=" << rhole << " rCyl=" << rCyl ;
303 // std::cout << " zb=" << zb << " zf=" << zf << " epsl=" << epsl << std::endl;
304 //--- debug end
305
306 // Check of Charge
307
308 if ( int( _charge ) == 0 )
309 {
310 std::cout << "HelCyl gets a straight line !!" << std::endl;
311 return -1;
312 }
313
314 // parameters
315
316 HepPoint3D CenterCyl( 0., 0., 0. );
317 HepPoint3D CenterTrk = _helix->center();
318 double rTrk = fabs( _helix->radius() );
319
320 double rPivot = fabs( _helix->pivot().perp() );
321
322 double phi0 = _helix->a()[1];
323 double tanl = _helix->a()[4];
324 // double zdz = _helix->pivot().z() + _helix->a()[3];
325 double dPhi;
326 double zee;
327
328 // Calculate intersections between cylinder and track
329 // if return value = 2 track hitting barrel part
330
331 HepPoint3D Cross1, Cross2;
332
333 if ( intersection( CenterTrk, _charge * rTrk, CenterCyl, rCyl, epsl, Cross1, Cross2 ) == 2 )
334 {
335
336 double phiCyl = atan2( _charge * ( CenterTrk.y() - Cross1.y() ),
337 _charge * ( CenterTrk.x() - Cross1.x() ) );
338 phiCyl = ( phiCyl > 0. ) ? phiCyl : phiCyl + 2. * M_PI;
339
340 dPhi = phiCyl - phi0;
341
342 // dPhi region ( at cylinder )
343 // -pi <= dPhi < pi
344
345 double PhiYobun = 1. / fabs( _helix->radius() );
346 double zero = 0.00001;
347
348 if ( _charge >= 0. )
349 {
350 if ( dPhi > PhiYobun ) dPhi -= 2. * M_PI;
351 if ( -2. * M_PI - zero <= dPhi <= ( -2. * M_PI + PhiYobun ) ) dPhi += 2. * M_PI;
352 }
353
354 if ( _charge < 0. )
355 {
356 if ( dPhi < -PhiYobun ) dPhi += 2. * M_PI;
357 if ( 2. * M_PI + zero >= dPhi >= ( 2. * M_PI - PhiYobun ) ) dPhi -= 2. * M_PI;
358 }
359
360 if ( dPhi < -M_PI ) dPhi += 2. * M_PI;
361 if ( dPhi >= M_PI ) dPhi -= 2. * M_PI;
362
363 //--debug
364 // std::cout << "dPhi = " << dPhi << std::endl;
365 //--debug end
366
367 xp.setX( Cross1.x() );
368 xp.setY( Cross1.y() );
369 xp.setZ( _helix->x( dPhi ).z() );
370 // xp.setZ( zdz - _charge * rTrk * tanl * dPhi );
371
372 if ( xp.z() > zb && xp.z() < zf )
373 {
374 phi = dPhi;
375 //--- debug ---
376 // std::cout << "return1 ( ext success )" << std::endl;
377 // std::cout << " xp:= " << xp.x() << ", " << xp.y() << ", " << xp.z() << std::endl;
378 //--- debug ----
379 return 1;
380 }
381 }
382
383 // tracks hitting endcaps
384
385 if ( fabs( tanl ) < 0.1 )
386 {
387 //--- debug ---
388 // std::cout << "return0 ( ext failed , |tanl| < 0.1 )" << std::endl;
389 // std::cout << " xp:= " << xp.x() << ", " << xp.y() << ", " << xp.z() << std::endl;
390 //--- debug ---
391 return 0;
392 }
393
394 if ( tanl > 0. )
395 {
396 zee = zf;
397 status = 3;
398 }
399 else
400 {
401 zee = zb;
402 status = 2;
403 }
404
405 dPhi = _charge * ( _helix->x( 0. ).z() - zee ) / rTrk / tanl;
406 // dPhi = _charge * ( zdz - zee )/rTrk/tanl;
407
408 // Requre dPhi < 2*pi
409
410 if ( fabs( dPhi ) > 2. * M_PI )
411 {
412 //--- debug ---
413 // std::cout << " return0 ( ext failed , dPhi > 2pi )" << std::endl;
414 // std::cout << " xp:= " << xp.x() << ", " << xp.y() << ", " << xp.z() << std::endl;
415 //--- debug ---
416 return 0;
417 }
418
419 xp.setX( _helix->x( dPhi ).x() );
420 xp.setY( _helix->x( dPhi ).y() );
421 xp.setZ( zee );
422
423 double test = xp.perp2() - rhole * rhole;
424 if ( test < 0. )
425 {
426 //--- debug ---
427 // std::cout << "return0 ( cannot reach to rhole at z=edge )" << std::endl;
428 // std::cout << " xp:= " << xp.x() << ", " << xp.y() << ", " << xp.z() << std::endl;
429 //--- debug ---
430 return 0;
431 }
432
433 phi = dPhi;
434 //--- debug ---
435 // std::cout << "return" << status << std::endl;
436 // std::cout << " xp:= " << xp.x() << ", " << xp.y() << ", " << xp.z() << std::endl;
437 //--- debug ---
438
439 return status;
440}
441
442/*
443void
444TTrack::findCatHit(unsigned trackid) {
445
446 unsigned i = 0;
447 // while ( TMDC01CatHit * chit = _cathits.last() ) {
448 // _cathits.remove(chit);
449 // delete chit;
450 // }
451
452 //--- debug ---
453 //std::cout << std::endl << "FindCatHit called !! " << std::endl << std::endl;
454 //--- debug end ----
455
456 // set cylinder geometry
457 double rcyl[3];
458 rcyl[0] = 8.80 ;
459 rcyl[1] = 9.80 ;
460 rcyl[2] = 10.85 ;
461 double rhole = 8.00 ;
462 double zb = -26.17 ;
463 double zf = 45.87 ;
464 double epsl = 0.01 ;
465
466 //
467 HepPoint3D xp ;
468 double phi ;
469 TMDCCatHit * chit;
470
471 // loop over layers and find cathits
472
473 for ( unsigned layer = 0 ; layer < 3 ; layer++ ) {
474 if ( HelCyl (rhole, rcyl[layer], zb,zf,epsl, phi, xp ) == 1 ) {
475 chit = new TMDCCatHit( this, trackid, layer, xp.x(), xp.y(), xp.z()) ;
476 _catHits.append(chit);
477 }
478 }
479
480 //--- debug ---
481 //if(_cathits.length()) {
482 // std::cout << std::endl;
483 // for ( int j = 0 ; j < _cathits.length() ; j++ ) {
484 // _cathits[j]->dump(" ", " ") ;
485 // }
486 //}
487 // chit->dump(" ", " ") ;
488 //--- debug end ----
489
490}
491//--- Nagoya mods. 19981225 end ---------------------------
492
493
494// int
495// TTrack::fitx(void) {
496
497// #ifdef TRKRECO_DEBUG_DETAIL
498// std::cout << "TTrack::fit !!! This is obsolete !!!" << std::endl;
499// #endif
500// if (_fitted) return 0;
501
502// //...Check # of hits...
503// if (_links.length() < 5) return -1;
504
505// //...Setup...
506// unsigned nTrial = 0;
507// Vector a(5), da(5);
508// a = _helix->a();
509// Vector dxda(5);
510// Vector dyda(5);
511// Vector dzda(5);
512// Vector dDda(5);
513// Vector dchi2da(5);
514// SymMatrix d2chi2d2a(5, 0);
515// double chi2;
516// double chi2Old = 10e99;
517// const double convergence = 1.0e-5;
518// bool allAxial = true;
519// Matrix e(3, 3);
520// Vector f(3);
521// int err = 0;
522// double factor = 1.0;//jtanaka0715
523
524// Vector maxDouble(5);
525// for (unsigned i = 0; i < 5; i++) maxDouble[i] = (FLT_MAX);
526
527// //...Fitting loop...
528// while (nTrial < 100) {
529
530// //...Set up...
531// chi2 = 0.;
532// for (unsigned j = 0; j < 5; j++) dchi2da[j] = 0.;
533// d2chi2d2a = SymMatrix(5, 0);
534
535// //...Loop with hits...
536// unsigned i = 0;
537// while (TMLink * l = _links[i++]) {
538// const TMDCWireHit & h = * l->hit();
539
540// //...Check state...
541// if (h.state() & WireHitInvalidForFit) continue;
542
543// //...Check wire...
544// if (! nTrial)
545// if (h.wire()->stereo()) allAxial = false;
546
547// //...Cal. closest points...
548// approach(* l);
549// double dPhi = l->dPhi();
550// const HepPoint3D & onTrack = l->positionOnTrack();
551// const HepPoint3D & onWire = l->positionOnWire();
552
553// #ifdef TRKRECO_DEBUG_DETAIL
554// // std::cout << " in fit " << onTrack << ", " << onWire;
555// // h.dump();
556// #endif
557
558// //...Obtain drift distance and its error...
559// unsigned leftRight = WireHitRight;
560// if (onWire.cross(onTrack).z() < 0.) leftRight = WireHitLeft;
561// double distance = h.distance(leftRight);
562// double eDistance = h.dDistance(leftRight);
563// double eDistance2 = eDistance * eDistance;
564
565// //...Residual...
566// HepVector3D v = onTrack - onWire;
567// double vmag = v.mag();
568// double dDistance = vmag - distance;
569
570// //...dxda...
571// this->dxda(* l, dPhi, dxda, dyda, dzda);
572
573// //...Chi2 related...
574// // dDda = (v.x() * dxda + v.y() * dyda + v.z() * dzda) / vmag;
575// Vector3 vw = h.wire()->direction();
576// dDda = (vmag > 0.)
577// ? ((v.x() * (1. - vw.x() * vw.x()) -
578// v.y() * vw.x() * vw.y() - v.z() * vw.x() * vw.z())
579// * dxda +
580// (v.y() * (1. - vw.y() * vw.y()) -
581// v.z() * vw.y() * vw.z() - v.x() * vw.y() * vw.x())
582// * dyda +
583// (v.z() * (1. - vw.z() * vw.z()) -
584// v.x() * vw.z() * vw.x() - v.y() * vw.z() * vw.y())
585// * dzda) / vmag
586// : Vector(5,0);
587// if(vmag<=0.0) {
588// std::cout << " in fit " << onTrack << ", " << onWire;
589// h.dump();
590// }
591// dchi2da += (dDistance / eDistance2) * dDda;
592// d2chi2d2a += vT_times_v(dDda) / eDistance2;
593// double pChi2 = dDistance * dDistance / eDistance2;
594// chi2 += pChi2;
595
596// //...Store results...
597// l->update(onTrack, onWire, leftRight, pChi2);
598// }
599
600// //...Check condition...
601// double change = chi2Old - chi2;
602// if (fabs(change) < convergence) break;
603// //if (change < 0.) break;
604// //jtanaka -- from traffs -- Ozaki-san added this part to traffs.
605// if (change < 0.){
606// #ifdef TRKRECO_DEBUG_DETAIL
607// std::cout << "chi2Old, chi2=" << chi2Old <<" "<< chi2 << std::endl;
608// #endif
609// //change to the old value.
610// a += factor*da;
611// _helix->a(a);
612
613// chi2 = 0.;
614// for (unsigned j = 0; j < 5; j++) dchi2da[j] = 0.;
615// d2chi2d2a = SymMatrix(5, 0);
616
617// //...Loop with hits...
618// unsigned i = 0;
619// while (TMLink * l = _links[i++]) {
620// const TMDCWireHit & h = * l->hit();
621
622// //...Check wire...
623// if (! nTrial)
624// if (h.wire()->stereo()) allAxial = false;
625
626// //...Cal. closest points...
627// approach(* l);
628// double dPhi = l->dPhi();
629// const HepPoint3D & onTrack = l->positionOnTrack();
630// const HepPoint3D & onWire = l->positionOnWire();
631
632// #ifdef TRKRECO_DEBUG_DETAIL
633// // std::cout << " in fit in case of change < 0. " << onTrack << ", " <<
634onWire;
635// // h.dump();
636// #endif
637
638// //...Obtain drift distance and its error...
639// unsigned leftRight = WireHitRight;
640// if (onWire.cross(onTrack).z() < 0.) leftRight = WireHitLeft;
641// double distance = h.distance(leftRight);
642// double eDistance = h.dDistance(leftRight);
643// double eDistance2 = eDistance * eDistance;
644
645// //...Residual...
646// HepVector3D v = onTrack - onWire;
647// double vmag = v.mag();
648// double dDistance = vmag - distance;
649
650// //...dxda...
651// this->dxda(* l, dPhi, dxda, dyda, dzda);
652
653// //...Chi2 related...
654// //dDda = (v.x() * dxda + v.y() * dyda + v.z() * dzda) / vmag;
655// Vector3 vw = h.wire()->direction();
656// dDda = (vmag > 0.)
657// ? ((v.x() * (1. - vw.x() * vw.x()) -
658// v.y() * vw.x() * vw.y() - v.z() * vw.x() * vw.z())
659// * dxda +
660// (v.y() * (1. - vw.y() * vw.y()) -
661// v.z() * vw.y() * vw.z() - v.x() * vw.y() * vw.x())
662// * dyda +
663// (v.z() * (1. - vw.z() * vw.z()) -
664// v.x() * vw.z() * vw.x() - v.y() * vw.z() * vw.y())
665// * dzda) / vmag
666// : Vector(5,0);
667// if(vmag<=0.0) {
668// std::cout << " in fit " << onTrack << ", " << onWire;
669// h.dump();
670// }
671// dchi2da += (dDistance / eDistance2) * dDda;
672// d2chi2d2a += vT_times_v(dDda) / eDistance2;
673// double pChi2 = dDistance * dDistance / eDistance2;
674// chi2 += pChi2;
675
676// //...Store results...
677// l->update(onTrack, onWire, leftRight, pChi2);
678// }
679// //break;
680// factor *= 0.75;
681// #ifdef TRKRECO_DEBUG_DETAIL
682// std::cout << "factor = " << factor << std::endl;
683// std::cout << "chi2 = " << chi2 << std::endl;
684// #endif
685// if(factor < 0.01)break;
686// }
687
688// chi2Old = chi2;
689
690// //...Cal. helix parameters for next loop...
691// if (allAxial) {
692// f = dchi2da.sub(1, 3);
693// e = d2chi2d2a.sub(1, 3);
694// f = solve(e, f);
695// da[0] = f[0];
696// da[1] = f[1];
697// da[2] = f[2];
698// da[3] = 0.;
699// da[4] = 0.;
700// }
701// else {
702// da = solve(d2chi2d2a, dchi2da);
703// }
704
705// #ifdef TRKRECO_DEBUG_DETAIL
706// //std::cout << " fit " << nTrial << " : " << da << std::endl;
707// //std::cout << " d2chi " << d2chi2d2a << std::endl;
708// //std::cout << " dchi2 " << dchi2da << std::endl;
709// #endif
710
711// a -= factor*da;
712// _helix->a(a);
713// ++nTrial;
714// }
715
716// //...Cal. error matrix...
717// SymMatrix Ea(5, 0);
718// unsigned dim;
719// if (allAxial) {
720// dim = 3;
721// SymMatrix Eb(3, 0), Ec(3, 0);
722// Eb = d2chi2d2a.sub(1, 3);
723// Ec = Eb.inverse(err);
724// Ea[0][0] = Ec[0][0];
725// Ea[0][1] = Ec[0][1];
726// Ea[0][2] = Ec[0][2];
727// Ea[1][1] = Ec[1][1];
728// Ea[1][2] = Ec[1][2];
729// Ea[2][2] = Ec[2][2];
730// }
731// else {
732// dim = 5;
733// Ea = d2chi2d2a.inverse(err);
734// }
735
736// //...Store information...
737// if (! err) {
738// _helix->a(a);
739// _helix->Ea(Ea);
740// }
741
742// _ndf = _links.length() - dim;
743// _chi2 = chi2;
744
745// _fitted = true;
746// return err;
747// }
748
749*/
750
751int TTrack::dxda( const TMLink& link, double dPhi, Vector& dxda, Vector& dyda,
752 Vector& dzda ) const {
753
754 //...Setup...
755 const TMDCWire& w = *link.wire();
756 Vector a = _helix->a();
757 double dRho = a[0];
758 double phi0 = a[1];
759 double kappa = a[2];
760 // double rho = Helix::ConstantAlpha / kappa;
761 // double rho = 333.564095 / kappa;
762 const double Bz = -1000 * getPmgnIMF()->getReferField();
763 double rho = 333.564095 / ( kappa * Bz );
764 double tanLambda = a[4];
765
766 double sinPhi0 = sin( phi0 );
767 double cosPhi0 = cos( phi0 );
768 double sinPhi0dPhi = sin( phi0 + dPhi );
769 double cosPhi0dPhi = cos( phi0 + dPhi );
770 Vector dphida( 5 );
771
772 // yzhang 2012-08-30
773 //...Axial case...
774 // if (w.axial()) {
775 // HepPoint3D d = _helix->center() - w.xyPosition();
776 // double dmag2 = d.mag2();
777 //
778 // dphida[0] = (sinPhi0 * d.x() - cosPhi0 * d.y()) / dmag2;
779 // dphida[1] = (dRho + rho) * (cosPhi0 * d.x() + sinPhi0 * d.y())
780 // / dmag2 - 1.;
781 // dphida[2] = (- rho / kappa) * (sinPhi0 * d.x() - cosPhi0 * d.y())
782 // / dmag2;
783 // dphida[3] = 0.;
784 // dphida[4] = 0.;
785 // }
786 //
787 // //...Stereo case...
788 // else {
789 HepPoint3D onTrack = _helix->x( dPhi );
790 HepVector3D v = w.direction();
791 Vector c( 3 );
792 c = HepPoint3D( w.backwardPosition() - ( v * w.backwardPosition() ) * v );
793
794 Vector x( 3 );
795 x = onTrack;
796
797 Vector dxdphi( 3 );
798 dxdphi[0] = rho * sinPhi0dPhi;
799 dxdphi[1] = -rho * cosPhi0dPhi;
800 dxdphi[2] = -rho * tanLambda;
801
802 Vector d2xdphi2( 3 );
803 d2xdphi2[0] = rho * cosPhi0dPhi;
804 d2xdphi2[1] = rho * sinPhi0dPhi;
805 d2xdphi2[2] = 0.;
806
807 double dxdphi_dot_v = dxdphi[0] * v.x() + dxdphi[1] * v.y() + dxdphi[2] * v.z();
808 double x_dot_v = x[0] * v.x() + x[1] * v.y() + x[2] * v.z();
809
810 double dfdphi = -( dxdphi[0] - dxdphi_dot_v * v.x() ) * dxdphi[0] -
811 ( dxdphi[1] - dxdphi_dot_v * v.y() ) * dxdphi[1] -
812 ( dxdphi[2] - dxdphi_dot_v * v.z() ) * dxdphi[2] -
813 ( x[0] - c[0] - x_dot_v * v.x() ) * d2xdphi2[0] -
814 ( x[1] - c[1] - x_dot_v * v.y() ) * d2xdphi2[1];
815 /* - (x[2] - c[2] - x_dot_v*v.z()) * d2xdphi2[2]; = 0. */
816
817 // dxda_phi, dyda_phi, dzda_phi : phi is fixed
818 Vector dxda_phi( 5 );
819 dxda_phi[0] = cosPhi0;
820 dxda_phi[1] = -( dRho + rho ) * sinPhi0 + rho * sinPhi0dPhi;
821 dxda_phi[2] = -( rho / kappa ) * ( cosPhi0 - cosPhi0dPhi );
822 dxda_phi[3] = 0.;
823 dxda_phi[4] = 0.;
824
825 Vector dyda_phi( 5 );
826 dyda_phi[0] = sinPhi0;
827 dyda_phi[1] = ( dRho + rho ) * cosPhi0 - rho * cosPhi0dPhi;
828 dyda_phi[2] = -( rho / kappa ) * ( sinPhi0 - sinPhi0dPhi );
829 dyda_phi[3] = 0.;
830 dyda_phi[4] = 0.;
831
832 Vector dzda_phi( 5 );
833 dzda_phi[0] = 0.;
834 dzda_phi[1] = 0.;
835 dzda_phi[2] = ( rho / kappa ) * tanLambda * dPhi;
836 dzda_phi[3] = 1.;
837 dzda_phi[4] = -rho * dPhi;
838
839 Vector d2xdphida( 5 );
840 d2xdphida[0] = 0.;
841 d2xdphida[1] = rho * cosPhi0dPhi;
842 d2xdphida[2] = -( rho / kappa ) * sinPhi0dPhi;
843 d2xdphida[3] = 0.;
844 d2xdphida[4] = 0.;
845
846 Vector d2ydphida( 5 );
847 d2ydphida[0] = 0.;
848 d2ydphida[1] = rho * sinPhi0dPhi;
849 d2ydphida[2] = ( rho / kappa ) * cosPhi0dPhi;
850 d2ydphida[3] = 0.;
851 d2ydphida[4] = 0.;
852
853 Vector d2zdphida( 5 );
854 d2zdphida[0] = 0.;
855 d2zdphida[1] = 0.;
856 d2zdphida[2] = ( rho / kappa ) * tanLambda;
857 d2zdphida[3] = 0.;
858 d2zdphida[4] = -rho;
859
860 Vector dfda( 5 );
861 for ( int i = 0; i < 5; i++ )
862 {
863 double d_dot_v = v.x() * dxda_phi[i] + v.y() * dyda_phi[i] + v.z() * dzda_phi[i];
864 dfda[i] = -( dxda_phi[i] - d_dot_v * v.x() ) * dxdphi[0] -
865 ( dyda_phi[i] - d_dot_v * v.y() ) * dxdphi[1] -
866 ( dzda_phi[i] - d_dot_v * v.z() ) * dxdphi[2] -
867 ( x[0] - c[0] - x_dot_v * v.x() ) * d2xdphida[i] -
868 ( x[1] - c[1] - x_dot_v * v.y() ) * d2ydphida[i] -
869 ( x[2] - c[2] - x_dot_v * v.z() ) * d2zdphida[i];
870 dphida[i] = -dfda[i] / dfdphi;
871 }
872 //}
873
874 dxda[0] = cosPhi0 + rho * sinPhi0dPhi * dphida[0];
875 dxda[1] = -( dRho + rho ) * sinPhi0 + rho * sinPhi0dPhi * ( 1. + dphida[1] );
876 dxda[2] = -rho / kappa * ( cosPhi0 - cosPhi0dPhi ) + rho * sinPhi0dPhi * dphida[2];
877 dxda[3] = rho * sinPhi0dPhi * dphida[3];
878 dxda[4] = rho * sinPhi0dPhi * dphida[4];
879
880 dyda[0] = sinPhi0 - rho * cosPhi0dPhi * dphida[0];
881 dyda[1] = ( dRho + rho ) * cosPhi0 - rho * cosPhi0dPhi * ( 1. + dphida[1] );
882 dyda[2] = -rho / kappa * ( sinPhi0 - sinPhi0dPhi ) - rho * cosPhi0dPhi * dphida[2];
883 dyda[3] = -rho * cosPhi0dPhi * dphida[3];
884 dyda[4] = -rho * cosPhi0dPhi * dphida[4];
885
886 dzda[0] = -rho * tanLambda * dphida[0];
887 dzda[1] = -rho * tanLambda * dphida[1];
888 dzda[2] = rho / kappa * tanLambda * dPhi - rho * tanLambda * dphida[2];
889 dzda[3] = 1. - rho * tanLambda * dphida[3];
890 dzda[4] = -rho * dPhi - rho * tanLambda * dphida[4];
891
892 return 0;
893}
894
895#define NEW_FIT2D 1
896#if !( NEW_FIT2D )
897int TTrack::fit2D( void ) {
898# ifdef TRKRECO_DEBUG_DETAIL
899 std::cout << " TTrack::fit2D(r-phi) ..." << std::endl;
900# endif
901 if ( _fitted ) return 0;
902
903 //...Check # of hits...
904 if ( _links.length() < 4 ) return -1;
905 // std::cout << "# = " << _links.length() << std::endl;
906 //...Setup...
907 unsigned nTrial = 0;
908 Vector a = _helix->a();
909 Vector da( 5 ), dxda( 5 ), dyda( 5 ), dzda( 5 );
910 Vector dDda( 5 ), dchi2da( 5 );
911 SymMatrix d2chi2d2a( 5, 0 );
912 double chi2;
913 double chi2Old = 1.0e+99;
914 const double convergence = 1.0e-5;
915 Matrix e( 3, 3 );
916 Vector f( 3 );
917 int err = 0;
918 double factor = 1.0;
919
920 //...Fitting loop...
921 while ( nTrial < 100 )
922 {
923# ifdef TRKRECO_DEBUG_DETAIL
924 if ( a[3] != 0. || a[4] != 0. )
925 cerr << "Error in 2D functions of TTrack : a = " << a << std::endl;
926# endif
927 //...Set up...
928 chi2 = 0.;
929 for ( unsigned j = 0; j < 5; ++j ) dchi2da[j] = 0.;
930 d2chi2d2a = SymMatrix( 5, 0 );
931
932 //...Loop with hits...
933 unsigned i = 0;
934 while ( TMLink* l = _links[i++] )
935 {
936 const TMDCWireHit& h = *l->hit();
937
938 //...Cal. closest points...
939 approach2D( *l );
940 double dPhi = l->dPhi();
941 const HepPoint3D& onTrack = l->positionOnTrack();
942 const HepPoint3D& onWire = l->positionOnWire();
943# ifdef TRKRECO_DEBUG_DETAIL
944 if ( onTrack.z() != 0. )
945 cerr << "Error in 2D functions of TTrack : onTrack = " << onTrack << std::endl;
946 if ( onWire.z() != 0. )
947 cerr << "Error in 2D functions of TTrack : onWire = " << onWire << std::endl;
948# endif
949 HepPoint3D onTrack2( onTrack.x(), onTrack.y(), 0. );
950 HepPoint3D onWire2( onWire.x(), onWire.y(), 0. );
951
952 //...Obtain drift distance and its error...
953 unsigned leftRight = WireHitRight;
954 if ( onWire2.cross( onTrack2 ).z() < 0. ) leftRight = WireHitLeft;
955 double distance = h.distance( leftRight );
956 double eDistance = h.dDistance( leftRight );
957 double eDistance2 = eDistance * eDistance;
958
959 //...Residual...
960 HepVector3D v = onTrack2 - onWire2;
961 double vmag = v.mag();
962 double dDistance = vmag - distance;
963
964 //...dxda...
965 this->dxda2D( *l, dPhi, dxda, dyda, dzda );
966
967 //...Chi2 related...
968 // Vector3 vw(0.,0.,1.);
969 dDda = ( vmag > 0. ) ? ( v.x() * dxda + v.y() * dyda ) / vmag : Vector( 5, 0 );
970 // dDda = (vmag > 0.)
971 //? ((v.x() * (1. - vw.x() * vw.x()) -
972 // v.y() * vw.x() * vw.y() - v.z() * vw.x() * vw.z())
973 // * dxda +
974 // (v.y() * (1. - vw.y() * vw.y()) -
975 // v.z() * vw.y() * vw.z() - v.x() * vw.y() * vw.x())
976 // * dyda +
977 // (v.z() * (1. - vw.z() * vw.z()) -
978 // v.x() * vw.z() * vw.x() - v.y() * vw.z() * vw.y())
979 // * dzda) / vmag
980 //: Vector(5,0);
981 if ( vmag <= 0.0 )
982 {
983 std::cout << " in fit2D " << onTrack << ", " << onWire;
984 h.dump();
985 }
986 dchi2da += ( dDistance / eDistance2 ) * dDda;
987 d2chi2d2a += vT_times_v( dDda ) / eDistance2;
988 double pChi2 = dDistance * dDistance / eDistance2;
989 chi2 += pChi2;
990
991 //...Store results...
992 l->update( onTrack2, onWire2, leftRight, pChi2 );
993 }
994
995 //...Check condition...
996 double change = chi2Old - chi2;
997 if ( fabs( change ) < convergence ) break;
998 if ( change < 0. )
999 {
1000# ifdef TRKRECO_DEBUG_DETAIL
1001 std::cout << "chi2Old, chi2=" << chi2Old << " " << chi2 << std::endl;
1002# endif
1003 // change to the old value.
1004 a += factor * da;
1005 _helix->a( a );
1006
1007 chi2 = 0.;
1008 for ( unsigned j = 0; j < 5; ++j ) dchi2da[j] = 0.;
1009 d2chi2d2a = SymMatrix( 5, 0 );
1010
1011 //...Loop with hits...
1012 unsigned i = 0;
1013 while ( TMLink* l = _links[i++] )
1014 {
1015 const TMDCWireHit& h = *l->hit();
1016
1017 //...Cal. closest points...
1018 approach2D( *l );
1019 double dPhi = l->dPhi();
1020 const HepPoint3D& onTrack = l->positionOnTrack();
1021 const HepPoint3D& onWire = l->positionOnWire();
1022 HepPoint3D onTrack2( onTrack.x(), onTrack.y(), 0. );
1023 HepPoint3D onWire2( onWire.x(), onWire.y(), 0. );
1024
1025 //...Obtain drift distance and its error...
1026 unsigned leftRight = WireHitRight;
1027 if ( onWire2.cross( onTrack2 ).z() < 0. ) leftRight = WireHitLeft;
1028 double distance = h.distance( leftRight );
1029 double eDistance = h.dDistance( leftRight );
1030 double eDistance2 = eDistance * eDistance;
1031
1032 //...Residual...
1033 HepVector3D v = onTrack2 - onWire2;
1034 double vmag = v.mag();
1035 double dDistance = vmag - distance;
1036
1037 //...dxda...
1038 this->dxda2D( *l, dPhi, dxda, dyda, dzda );
1039
1040 //...Chi2 related...
1041 dDda = ( vmag > 0. ) ? ( v.x() * dxda + v.y() * dyda ) / vmag : Vector( 5, 0 );
1042 if ( vmag <= 0.0 )
1043 {
1044 std::cout << " in fit2D " << onTrack << ", " << onWire;
1045 h.dump();
1046 }
1047 dchi2da += ( dDistance / eDistance2 ) * dDda;
1048 d2chi2d2a += vT_times_v( dDda ) / eDistance2;
1049 double pChi2 = dDistance * dDistance / eDistance2;
1050 chi2 += pChi2;
1051
1052 //...Store results...
1053 l->update( onTrack2, onWire2, leftRight, pChi2 );
1054 }
1055 // break;
1056 factor *= 0.75;
1057# ifdef TRKRECO_DEBUG_DETAIL
1058 std::cout << "factor = " << factor << std::endl;
1059 std::cout << "chi2 = " << chi2 << std::endl;
1060# endif
1061 if ( factor < 0.01 ) break;
1062 }
1063
1064 chi2Old = chi2;
1065
1066 //...Cal. helix parameters for next loop...
1067 f = dchi2da.sub( 1, 3 );
1068 e = d2chi2d2a.sub( 1, 3 );
1069 f = solve( e, f );
1070 da[0] = f[0];
1071 da[1] = f[1];
1072 da[2] = f[2];
1073 da[3] = 0.;
1074 da[4] = 0.;
1075
1076 a -= factor * da;
1077 _helix->a( a );
1078 // std::cout << nTrial << ": chi2 = " << chi2
1079 // << ", helix = (" << a[0] << ", " << a[1]
1080 // << ", " << a[2] << ", " << a[3]
1081 // << ", " << a[4] << ")" << std::endl;
1082 ++nTrial;
1083 }
1084
1085 //...Cal. error matrix...
1086 SymMatrix Ea( 5, 0 );
1087 unsigned dim = 3;
1088 SymMatrix Eb = d2chi2d2a.sub( 1, 3 );
1089 SymMatrix Ec = Eb.inverse( err );
1090 Ea[0][0] = Ec[0][0];
1091 Ea[0][1] = Ec[0][1];
1092 Ea[0][2] = Ec[0][2];
1093 Ea[1][1] = Ec[1][1];
1094 Ea[1][2] = Ec[1][2];
1095 Ea[2][2] = Ec[2][2];
1096
1097 //...Store information...
1098 if ( !err )
1099 {
1100 _helix->a( a );
1101 _helix->Ea( Ea );
1102 }
1103
1104 _ndf = _links.length() - dim;
1105 _chi2 = chi2;
1106
1107 _fitted = true;
1108 return err;
1109}
1110#endif
1111
1112/*
1113int TTrack::fitWithCathode(float window, int SysCorr ) {
1114
1115#ifdef TRKRECO_DEBUG_DETAIL
1116 std::cout << " TTrack::fitCathode ..." << std::endl;
1117#endif
1118
1119 if (_fittedWithCathode) return 0;
1120
1121 //...Check # of hits...
1122 if (nCores() < 5) return -1;
1123
1124 //...for cathode ndf...
1125 int NusedCathode(0);
1126
1127 //...Setup...
1128 unsigned nTrial = 0;
1129 Vector a(5), da(5);
1130 a = _helix->a();
1131 Vector dxda(5);
1132 Vector dyda(5);
1133 Vector dzda(5);
1134 Vector dDda(5);
1135 Vector dDzda(5); // for cathode z
1136 Vector dchi2da(5);
1137 SymMatrix d2chi2d2a(5, 0);
1138 double chi2;
1139 double chi2Old = 10e99;
1140 const double convergence = 1.0e-5;
1141 bool allAxial = true;
1142 int LayerStat(0); // layer status axial:0 stereo:1 cathode:2
1143 Matrix e(3, 3);
1144 Vector f(3);
1145 int err = 0;
1146 double factor = 1.0;
1147
1148 const AList<TMDCCatHit> & chits = _catHits;
1149
1150 Vector maxDouble(5);
1151 for (unsigned i = 0; i < 5; i++) maxDouble[i] = (FLT_MAX);
1152
1153 //...Fitting loop...
1154 while (nTrial < 100) {
1155
1156 //...Set up...
1157 chi2 = 0.;
1158 NusedCathode = 0;
1159 for (unsigned j = 0; j < 5; j++) dchi2da[j] = 0.;
1160 d2chi2d2a = SymMatrix(5, 0);
1161
1162 //...Loop with hits...
1163 unsigned i = 0;
1164 while (TMLink * l = cores()[i++]) {
1165
1166 const TMDCWireHit & h = * l->hit();
1167
1168 // Check layer status ( cathode added )
1169 LayerStat = 0;
1170 if ( h.wire()->stereo() ) LayerStat = 1;
1171 unsigned nlayer = h.wire()->layerId();
1172 if (nlayer == 0 || nlayer == 1 || nlayer == 2 ) LayerStat = 2;
1173
1174 //...Check wire...
1175 if (! nTrial)
1176 if (h.wire()->stereo() || LayerStat == 2 ) allAxial = false;
1177
1178 //...Cal. closest points...
1179 approach(* l);
1180 double dPhi = l->dPhi();
1181 const HepPoint3D & onTrack = l->positionOnTrack();
1182 const HepPoint3D & onWire = l->positionOnWire();
1183
1184#ifdef TRKRECO_DEBUG_DETAIL
1185 // std::cout << " in fitCathode " << onTrack << ", " << onWire;
1186 // h.dump();
1187#endif
1188
1189 //...Obtain drift distance and its error...
1190 unsigned leftRight = WireHitRight;
1191 if (onWire.cross(onTrack).z() < 0.) leftRight = WireHitLeft;
1192 double distance = h.drift(leftRight);
1193 double eDistance = h.dDrift(leftRight);
1194 double eDistance2 = eDistance * eDistance;
1195
1196 //...Residual...
1197 HepVector3D v = onTrack - onWire;
1198 double vmag = v.mag();
1199 double dDistance = vmag - distance;
1200
1201 //...dxda...
1202 this->dxda(* l, dPhi, dxda, dyda, dzda);
1203
1204 // ... Chi2 related ...
1205 double pChi2(0.);
1206
1207 if( LayerStat == 0 || LayerStat == 1 ){
1208 // dDda = (v.x() * dxda + v.y() * dyda + v.z() * dzda) / vmag;
1209 Vector3 vw = h.wire()->direction();
1210 dDda = (vmag > 0.)
1211 ? ((v.x() * (1. - vw.x() * vw.x()) -
1212 v.y() * vw.x() * vw.y() - v.z() * vw.x() * vw.z())
1213 * dxda +
1214 (v.y() * (1. - vw.y() * vw.y()) -
1215 v.z() * vw.y() * vw.z() - v.x() * vw.y() * vw.x())
1216 * dyda +
1217 (v.z() * (1. - vw.z() * vw.z()) -
1218 v.x() * vw.z() * vw.x() - v.y() * vw.z() * vw.y())
1219 * dzda) / vmag
1220 : Vector(5,0);
1221 if(vmag<=0.0) {
1222 std::cout << " in fit " << onTrack << ", " << onWire;
1223 h.dump();
1224 }
1225 dchi2da += (dDistance / eDistance2) * dDda;
1226 d2chi2d2a += vT_times_v(dDda) / eDistance2;
1227 double pChi2 = dDistance * dDistance / eDistance2;
1228 chi2 += pChi2;
1229
1230 } else {
1231
1232
1233 dDda = ( v.x() * dxda + v.y() * dyda )/v.perp();
1234 dchi2da += (dDistance/eDistance2) * dDda;
1235 d2chi2d2a += vT_times_v(dDda)/eDistance2;
1236
1237 pChi2 = 0;
1238 pChi2 += (dDistance/eDistance) * (dDistance/eDistance) ;
1239
1240 if ( l->usecathode() >= 3 ) {
1241
1242 TMDCClust * mclust = l->getmclust();
1243
1244 double dDistanceZ(this->helix().x(dPhi).z());
1245
1246 if( SysCorr ){
1247 if( !nTrial ) {
1248 mclust->zcalc( atan( this->helix().tanl()) );
1249 mclust->esterz( atan( this->helix().tanl()) );
1250 }
1251
1252 dDistanceZ -= mclust->zclustWithSysCorr();
1253 } else {
1254 dDistanceZ -= mclust->zclust();
1255 }
1256
1257 double eDistanceZ(mclust->erz());
1258 if( !eDistanceZ ) eDistanceZ = 99999.;
1259
1260 double eDistance2Z = eDistanceZ * eDistanceZ;
1261 double pChi2z = 0;
1262 pChi2z = (dDistanceZ/eDistanceZ);
1263 pChi2z *= pChi2z;
1264
1265#ifdef TRKRECO_DEBUG_DETAIL
1266 std::cout << "dDistanceZ = " << dDistanceZ << std::endl;
1267#endif
1268
1269 //.... requirement for use of cluster
1270
1271 if( nTrial == 0 &&
1272 fabs(dDistanceZ)< window &&
1273 mclust->chits().length() == 1
1274 ) l->setusecathode(4);
1275
1276 if( l->usecathode() == 4 ){
1277 NusedCathode++;
1278 dDzda = dzda ;
1279 dchi2da += (dDistanceZ/eDistance2Z) * dDzda;
1280 d2chi2d2a += vT_times_v(dDzda)/eDistance2Z;
1281 pChi2 += pChi2z ;
1282
1283 }
1284 }
1285
1286 } // end Chi2 related
1287
1288 chi2 += pChi2;
1289
1290
1291 //...Store results...
1292 l->update(onTrack, onWire, leftRight, pChi2);
1293
1294
1295 } // Tlink loop end
1296
1297 //...Check condition...
1298 double change = chi2Old - chi2;
1299 //if(TRKRECO_DEBUG_DETAIL>0) std::cout << "chi2 change = " <<change << std::endl;
1300
1301 if (fabs(change) < convergence) break;
1302 if (change < 0.) {
1303
1304#ifdef TRKRECO_DEBUG_DETAIL
1305 std::cout << "chi2Old, chi2=" << chi2Old <<" "<< chi2 << std::endl;
1306#endif
1307
1308 NusedCathode = 0;
1309 //change to the old value.
1310 a += factor*da;
1311 _helix->a(a);
1312
1313 chi2 = 0.;
1314 for (unsigned j = 0; j < 5; j++) dchi2da[j] = 0.;
1315 d2chi2d2a = SymMatrix(5, 0);
1316
1317
1318 //...Loop with hits...
1319 unsigned i = 0;
1320 while (TMLink * l = _links[i++]) {
1321 const TMDCWireHit & h = * l->hit();
1322
1323 // Check layer status ( cathode added )
1324 LayerStat = 0;
1325 if ( h.wire()->stereo() ) LayerStat = 1;
1326 unsigned nlayer = h.wire()->layerId();
1327 if (nlayer == 0 || nlayer == 1 || nlayer == 2 ) LayerStat = 2;
1328
1329 //...Cal. closest points...
1330 approach(* l);
1331 double dPhi = l->dPhi();
1332 const HepPoint3D & onTrack = l->positionOnTrack();
1333 const HepPoint3D & onWire = l->positionOnWire();
1334
1335#ifdef TRKRECO_DEBUG_DETAIL
1336 // std::cout << " in fitCathode " << onTrack << ", " << onWire;
1337 // h.dump();
1338#endif
1339
1340 //...Obtain drift distance and its error...
1341 unsigned leftRight = WireHitRight;
1342 if (onWire.cross(onTrack).z() < 0.) leftRight = WireHitLeft;
1343 double distance = h.drift(leftRight);
1344 double eDistance = h.dDrift(leftRight);
1345 double eDistance2 = eDistance * eDistance;
1346
1347 //...Residual...
1348 HepVector3D v = onTrack - onWire;
1349 double vmag = v.mag();
1350 double dDistance = vmag - distance;
1351
1352 //...dxda...
1353 this->dxda(* l, dPhi, dxda, dyda, dzda);
1354
1355 // ... Chi2 related ...
1356 double pChi2(0.);
1357
1358 if( LayerStat == 0 || LayerStat == 1 ){
1359 // dDda = (v.x() * dxda + v.y() * dyda + v.z() * dzda) / vmag;
1360 Vector3 vw = h.wire()->direction();
1361 dDda = (vmag > 0.)
1362 ? ((v.x() * (1. - vw.x() * vw.x()) -
1363 v.y() * vw.x() * vw.y() - v.z() * vw.x() * vw.z())
1364 * dxda +
1365 (v.y() * (1. - vw.y() * vw.y()) -
1366 v.z() * vw.y() * vw.z() - v.x() * vw.y() * vw.x())
1367 * dyda +
1368 (v.z() * (1. - vw.z() * vw.z()) -
1369 v.x() * vw.z() * vw.x() - v.y() * vw.z() * vw.y())
1370 * dzda) / vmag
1371 : Vector(5,0);
1372 if(vmag<=0.0) {
1373 std::cout << " in fit " << onTrack << ", " << onWire;
1374 h.dump();
1375 }
1376 dchi2da += (dDistance / eDistance2) * dDda;
1377 d2chi2d2a += vT_times_v(dDda) / eDistance2;
1378 double pChi2 = dDistance * dDistance / eDistance2;
1379 chi2 += pChi2;
1380
1381 } else {
1382
1383 dDda = ( v.x() * dxda + v.y() * dyda )/v.perp();
1384 dchi2da += (dDistance/eDistance2) * dDda;
1385 d2chi2d2a += vT_times_v(dDda)/eDistance2;
1386
1387 pChi2 = 0;
1388 pChi2 += (dDistance/eDistance) * (dDistance/eDistance) ;
1389
1390 if( l->usecathode() == 4 ){
1391
1392 TMDCClust * mclust = l->getmclust();
1393
1394 if( mclust ){
1395 NusedCathode++;
1396
1397 double dDistanceZ(this->helix().x(dPhi).z());
1398 if( SysCorr ) dDistanceZ -= mclust->zclustWithSysCorr();
1399 else dDistanceZ -= mclust->zclust();
1400
1401 double eDistanceZ(99999.);
1402 if( mclust->erz() != 0. ) eDistanceZ = mclust->erz();
1403
1404 double eDistance2Z = eDistanceZ * eDistanceZ;
1405 double pChi2z = 0;
1406 pChi2z = (dDistanceZ/eDistanceZ);
1407 pChi2z *= pChi2z;
1408
1409 dDzda = dzda ;
1410 dchi2da += (dDistanceZ/eDistance2Z) * dDzda;
1411 d2chi2d2a += vT_times_v(dDzda)/eDistance2Z;
1412 pChi2 += pChi2z ;
1413 }
1414
1415 }
1416
1417 } // end Chi2 related
1418
1419 chi2 += pChi2;
1420
1421
1422 //...Store results...
1423 l->update(onTrack, onWire, leftRight, pChi2);
1424
1425 }
1426
1427 //break;
1428
1429 factor *= 0.75;
1430#ifdef TRKRECO_DEBUG_DETAIL
1431 std::cout << "factor = " << factor << std::endl;
1432 std::cout << "chi2 = " << chi2 << std::endl;
1433#endif
1434 if(factor < 0.01)break;
1435
1436 }
1437
1438 chi2Old = chi2;
1439
1440 //...Cal. helix parameters for next loop...
1441 if (allAxial) {
1442 f = dchi2da.sub(1, 3);
1443 e = d2chi2d2a.sub(1, 3);
1444 f = solve(e, f);
1445 da[0] = f[0];
1446 da[1] = f[1];
1447 da[2] = f[2];
1448 da[3] = 0.;
1449 da[4] = 0.;
1450 }
1451 else {
1452 da = solve(d2chi2d2a, dchi2da);
1453 }
1454
1455#ifdef TRKRECO_DEBUG_DETAIL
1456 //std::cout << " fit " << nTrial << " : " << da << std::endl;
1457 //std::cout << " d2chi " << d2chi2d2a << std::endl;
1458 //std::cout << " dchi2 " << dchi2da << std::endl;
1459#endif
1460
1461 a -= da;
1462 _helix->a(a);
1463 ++nTrial;
1464 }
1465
1466 //...Cal. error matrix...
1467 SymMatrix Ea(5, 0);
1468 unsigned dim;
1469 if (allAxial) {
1470 dim = 3;
1471 SymMatrix Eb(3, 0), Ec(3, 0);
1472 Eb = d2chi2d2a.sub(1, 3);
1473 Ec = Eb.inverse(err);
1474 Ea[0][0] = Ec[0][0];
1475 Ea[0][1] = Ec[0][1];
1476 Ea[0][2] = Ec[0][2];
1477 Ea[1][1] = Ec[1][1];
1478 Ea[1][2] = Ec[1][2];
1479 Ea[2][2] = Ec[2][2];
1480 }
1481 else {
1482 dim = 5;
1483 Ea = d2chi2d2a.inverse(err);
1484 }
1485
1486 //...Store information...
1487 if (! err) {
1488 _helix->a(a);
1489 _helix->Ea(Ea);
1490 }
1491
1492 _ndf = nCores() + NusedCathode - dim;
1493 _chi2 = chi2;
1494
1495 _fittedWithCathode = true;
1496
1497 return err;
1498}
1499*/
1500
1501#if OLD_STEREO
1502int TTrack::stereoHitForCurl( TMLink& link, AList<HepPoint3D>& arcZList ) const {
1503 /*
1504 ATTENTION!!!!!
1505 Elements of arcZList are made by "new". We must delete them out of this function!!!!!
1506
1507 This function is used in "stereoHitForCurl(TMLink &link1, TMLink &link2)" and
1508 "stereoHitForCurl(TMLink &link1, TMLink &link2, TMLink &link3)" only.
1509 (I(jtanaka) checked it. --> no memory leak!!!)
1510
1511 */
1512
1513 // std::cout << "\n\nWire ID = " << link.hit()->wire()->id() << std::endl;
1514 const TMDCWireHit& h = *link.hit();
1515 HepVector3D X = 0.5 * ( h.wire()->forwardPosition() + h.wire()->backwardPosition() );
1516
1517 HepPoint3D center = _helix->center();
1518 HepPoint3D tmp( -999., -999., 0. );
1519 HepVector3D x = HepVector3D( X.x(), X.y(), 0. );
1520 HepVector3D w = x - center;
1521 HepVector3D V = h.wire()->direction();
1522 HepVector3D v = HepVector3D( V.x(), V.y(), 0. );
1523 double vmag2 = v.mag2();
1524 double vmag = sqrt( vmag2 );
1525 //...temporary
1526 link.position( tmp );
1527 arcZList.removeAll();
1528
1529 //...stereo?
1530 if ( vmag == 0. ) { return -1; }
1531
1532 double r = _helix->curv();
1533 double R[2];
1534 double drift = h.drift( WireHitLeft );
1535 R[0] = r + drift;
1536 R[1] = r - drift;
1537 double wv = w.dot( v );
1538 double d2[2];
1539 d2[0] = vmag2 * R[0] * R[0] +
1540 ( wv * wv - vmag2 * w.mag2() ); //...= v^2*(r^2 - w^2*sin()^2)...outer
1541 d2[1] = vmag2 * R[1] * R[1] +
1542 ( wv * wv - vmag2 * w.mag2() ); //...= v^2*(r^2 - w^2*sin()^2)...inner
1543
1544 //...No crossing in R/Phi plane...
1545 if ( d2[0] < 0. && d2[1] < 0. ) { return -1; }
1546
1547 bool ok_inner, ok_outer;
1548 ok_inner = true;
1549 ok_outer = true;
1550 double d[2];
1551 d[0] = -1.;
1552 d[1] = -1.;
1553 //...outer
1554 if ( d2[0] >= 0. ) { d[0] = sqrt( d2[0] ); }
1555 else { ok_outer = false; }
1556 if ( d2[1] >= 0. ) { d[1] = sqrt( d2[1] ); }
1557 else { ok_inner = false; }
1558
1559 //...Cal. length and z to crossing points...
1560 double l[2][2];
1561 double z[2][2];
1562 for (int i=0; i<2; i++){
1563 for (int j=0; j<2; j++){
1564 z[i][j] = 0.; // add initialize 25-05-15
1565 }
1566 }
1567 //...outer
1568 if ( ok_outer )
1569 {
1570 l[0][0] =
1571 ( -wv + d[0] ) / vmag2; //...= (-wvcos()+d)/v/v = (-wcos() + (r^2 - w^2*sin()^2)^0.5)/v
1572 l[1][0] =
1573 ( -wv - d[0] ) / vmag2; //...= (-wvcos()+d)/v/v = (-wcos() - (r^2 - w^2*sin()^2)^0.5)/v
1574 z[0][0] = X.z() + l[0][0] * V.z();
1575 z[1][0] = X.z() + l[1][0] * V.z();
1576 }
1577 //...inner
1578 if ( ok_inner )
1579 {
1580 l[0][1] =
1581 ( -wv + d[1] ) / vmag2; //...= (-wvcos()+d)/v/v = (-wcos() + (r^2 - w^2*sin()^2)^0.5)/v
1582 l[1][1] =
1583 ( -wv - d[1] ) / vmag2; //...= (-wvcos()+d)/v/v = (-wcos() - (r^2 - w^2*sin()^2)^0.5)/v
1584 z[0][1] = X.z() + l[0][1] * V.z();
1585 z[1][1] = X.z() + l[1][1] * V.z();
1586 }
1587
1588 //...Cal. xy position of crossing points...
1589 HepVector3D p[2][2];
1590 if ( ok_outer )
1591 {
1592 p[0][0] = x + l[0][0] * v;
1593 p[1][0] = x + l[1][0] * v;
1594 HepVector3D tmp_pc = p[0][0] - center;
1595 HepVector3D pc0 = HepVector3D( tmp_pc.x(), tmp_pc.y(), 0. );
1596 p[0][0] -= drift / pc0.mag() * pc0;
1597 tmp_pc = p[1][0] - center;
1598 HepVector3D pc1 = HepVector3D( tmp_pc.x(), tmp_pc.y(), 0. );
1599 p[1][0] -= drift / pc1.mag() * pc1;
1600 }
1601 if ( ok_inner )
1602 {
1603 p[0][1] = x + l[0][1] * v;
1604 p[1][1] = x + l[1][1] * v;
1605 HepVector3D tmp_pc = p[0][1] - center;
1606 HepVector3D pc0 = HepVector3D( tmp_pc.x(), tmp_pc.y(), 0. );
1607 p[0][1] += drift / pc0.mag() * pc0;
1608 tmp_pc = p[1][1] - center;
1609 HepVector3D pc1 = HepVector3D( tmp_pc.x(), tmp_pc.y(), 0. );
1610 p[1][1] += drift / pc1.mag() * pc1;
1611 }
1612
1613 //...boolean
1614 bool ok_xy[2][2];
1615 if ( ok_outer )
1616 {
1617 ok_xy[0][0] = true;
1618 ok_xy[1][0] = true;
1619 }
1620 else
1621 {
1622 ok_xy[0][0] = false;
1623 ok_xy[1][0] = false;
1624 }
1625 if ( ok_inner )
1626 {
1627 ok_xy[0][1] = true;
1628 ok_xy[1][1] = true;
1629 }
1630 else
1631 {
1632 ok_xy[0][1] = false;
1633 ok_xy[1][1] = false;
1634 }
1635 if ( ok_outer )
1636 {
1637 if ( _charge * ( center.x() * p[0][0].y() - center.y() * p[0][0].x() ) < 0. )
1638 ok_xy[0][0] = false;
1639 if ( _charge * ( center.x() * p[1][0].y() - center.y() * p[1][0].x() ) < 0. )
1640 ok_xy[1][0] = false;
1641 }
1642 if ( ok_inner )
1643 {
1644 if ( _charge * ( center.x() * p[0][1].y() - center.y() * p[0][1].x() ) < 0. )
1645 ok_xy[0][1] = false;
1646 if ( _charge * ( center.x() * p[1][1].y() - center.y() * p[1][1].x() ) < 0. )
1647 ok_xy[1][1] = false;
1648 }
1649 if ( !ok_inner && ok_outer && ( !ok_xy[0][0] ) && ( !ok_xy[1][0] ) ) { return -1; }
1650 if ( ok_inner && !ok_outer && ( !ok_xy[0][1] ) && ( !ok_xy[1][1] ) ) { return -1; }
1651
1652# if 0
1653 std::cout << "Drift Dist. = " << h.distance(WireHitLeft) << std::endl;
1654 std::cout << "outer ok? = " << ok_outer << std::endl;
1655 std::cout << "Z cand = " << z[0][0] << ", " << z[1][0] << std::endl;
1656 std::cout << "inner ok? = " << ok_inner << std::endl;
1657 std::cout << "Z cand = " << z[0][1] << ", " << z[1][1] << std::endl;
1658# endif
1659
1660 //...Check z position...
1661 if ( ok_xy[0][0] )
1662 {
1663 if ( z[0][0] < h.wire()->backwardPosition().z() ||
1664 z[0][0] > h.wire()->forwardPosition().z() )
1665 ok_xy[0][0] = false;
1666 }
1667 if ( ok_xy[1][0] )
1668 {
1669 if ( z[1][0] < h.wire()->backwardPosition().z() ||
1670 z[1][0] > h.wire()->forwardPosition().z() )
1671 ok_xy[1][0] = false;
1672 }
1673 if ( ok_xy[0][1] )
1674 {
1675 if ( z[0][1] < h.wire()->backwardPosition().z() ||
1676 z[0][1] > h.wire()->forwardPosition().z() )
1677 ok_xy[0][1] = false;
1678 }
1679 if ( ok_xy[1][1] )
1680 {
1681 if ( z[1][1] < h.wire()->backwardPosition().z() ||
1682 z[1][1] > h.wire()->forwardPosition().z() )
1683 ok_xy[1][1] = false;
1684 }
1685 if ( ( !ok_xy[0][0] ) && ( !ok_xy[1][0] ) && ( !ok_xy[0][1] ) && ( !ok_xy[1][1] ) )
1686 { return -1; }
1687
1688 double cosdPhi, dPhi;
1689
1690 if ( ok_xy[0][0] )
1691 {
1692 //...cal. arc length...
1693 cosdPhi = -center.dot( ( p[0][0] - center ).unit() ) / center.mag();
1694 if ( fabs( cosdPhi ) < 1.0 ) { dPhi = acos( cosdPhi ); }
1695 else if ( cosdPhi >= 1.0 ) { dPhi = 0.; }
1696 else { dPhi = M_PI; }
1697
1698 HepPoint3D* arcZ = new HepPoint3D;
1699 arcZ->setX( r * dPhi );
1700 arcZ->setY( z[0][0] );
1701 arcZList.append( arcZ );
1702 }
1703 if ( ok_xy[1][0] )
1704 {
1705 //...cal. arc length...
1706 cosdPhi = -center.dot( ( p[1][0] - center ).unit() ) / center.mag();
1707 if ( fabs( cosdPhi ) < 1.0 ) { dPhi = acos( cosdPhi ); }
1708 else if ( cosdPhi >= 1.0 ) { dPhi = 0.; }
1709 else { dPhi = M_PI; }
1710
1711 HepPoint3D* arcZ = new HepPoint3D;
1712 arcZ->setX( r * dPhi );
1713 arcZ->setY( z[1][0] );
1714 arcZList.append( arcZ );
1715 }
1716 if ( ok_xy[0][1] )
1717 {
1718 //...cal. arc length...
1719 cosdPhi = -center.dot( ( p[0][1] - center ).unit() ) / center.mag();
1720 if ( cosdPhi > 1.0 ) cosdPhi = 1.0;
1721 if ( cosdPhi < -1.0 ) cosdPhi = -1.0;
1722
1723 if ( fabs( cosdPhi ) < 1.0 ) { dPhi = acos( cosdPhi ); }
1724 else if ( cosdPhi >= 1.0 ) { dPhi = 0.; }
1725 else { dPhi = M_PI; }
1726
1727 HepPoint3D* arcZ = new HepPoint3D;
1728 arcZ->setX( r * dPhi );
1729 arcZ->setY( z[0][1] );
1730 arcZList.append( arcZ );
1731 }
1732 if ( ok_xy[1][1] )
1733 {
1734 //...cal. arc length...
1735 cosdPhi = -center.dot( ( p[1][1] - center ).unit() ) / center.mag();
1736 if ( cosdPhi > 1.0 ) cosdPhi = 1.0;
1737 if ( cosdPhi < -1.0 ) cosdPhi = -1.0;
1738 if ( fabs( cosdPhi ) < 1.0 ) { dPhi = acos( cosdPhi ); }
1739 else if ( cosdPhi >= 1.0 ) { dPhi = 0.; }
1740 else { dPhi = M_PI; }
1741
1742 HepPoint3D* arcZ = new HepPoint3D;
1743 arcZ->setX( r * dPhi );
1744 arcZ->setY( z[1][1] );
1745 arcZList.append( arcZ );
1746 }
1747
1748 if ( arcZList.length() == 1 || arcZList.length() == 2 ) { return 0; }
1749 else { return -1; }
1750}
1751
1752int TTrack::stereoHitForCurl( TMLink& link1, TMLink& link2 ) const {
1753 int flag1 = 0, flag2 = 0; // add initialize 25-05-15
1754 double minZ = 9999999.;
1755
1756 AList<HepPoint3D> arcZList1;
1757 AList<HepPoint3D> arcZList2;
1758 if ( stereoHitForCurl( link1, arcZList1 ) == 0 )
1759 {
1760 if ( stereoHitForCurl( link2, arcZList2 ) == 0 )
1761 {
1762 //...finds
1763 int n1 = arcZList1.length();
1764 int n2 = arcZList2.length();
1765 for ( int i = 0; i < n1; i++ )
1766 {
1767 for ( int j = 0; j < n2; j++ )
1768 {
1769 if ( fabs( arcZList1[i]->y() - arcZList2[j]->y() ) < minZ )
1770 {
1771 minZ = fabs( arcZList1[i]->y() - arcZList2[j]->y() );
1772 flag1 = i;
1773 flag2 = j;
1774 }
1775 }
1776 }
1777 }
1778 }
1779
1780 if ( minZ == 9999999. )
1781 {
1782 //...deletes
1783 deleteListForCurl( arcZList1, arcZList2 );
1784 return -1;
1785 }
1786
1787 //...sets the best values
1788 HepPoint3D tmp1 = *arcZList1[flag1];
1789 HepPoint3D tmp2 = *arcZList2[flag2];
1790 link1.position( tmp1 );
1791 link2.position( tmp2 );
1792 //...deletes
1793 deleteListForCurl( arcZList1, arcZList2 );
1794 return 0;
1795}
1796
1797# if defined( __GNUG__ )
1798int TArcOrder( const HepPoint3D** a, const HepPoint3D** b ) {
1799 if ( ( *a )->x() < ( *b )->x() ) { return 1; }
1800 else if ( ( *a )->x() == ( *b )->x() ) { return 0; }
1801 else { return -1; }
1802}
1803# else
1804extern "C" int TArcOrder( const void* av, const void* bv ) {
1805 const HepPoint3D** a( (const HepPoint3D**)av );
1806 const HepPoint3D** b( (const HepPoint3D**)bv );
1807 if ( ( *a )->x() < ( *b )->x() ) { return 1; }
1808 else if ( ( *a )->x() == ( *b )->x() ) { return 0; }
1809 else { return -1; }
1810}
1811# endif
1812
1813int TTrack::stereoHitForCurl( TMLink& link1, TMLink& link2, TMLink& link3 ) const {
1814 //...definition of the return values
1815 // -1 = error
1816 // 0 = normal = can use 3 links
1817 // 12 = can use 1 and 2 only
1818 // 23 = can use 2 and 3 only
1819
1820 int flag1 = 0, flag2 = 0, flag3 = 0; // add initialize 25-05-15
1821 double minZ1 = 9999999.;
1822 double minZ2 = 9999999.;
1823 double minZ01 = 9999999.;
1824 double minZ12 = 9999999.;
1825
1826 AList<HepPoint3D> arcZList1;
1827 AList<HepPoint3D> arcZList2;
1828 AList<HepPoint3D> arcZList3;
1829 int ok1 = stereoHitForCurl( link1, arcZList1 );
1830 int ok2 = stereoHitForCurl( link2, arcZList2 );
1831 int ok3 = stereoHitForCurl( link3, arcZList3 );
1832
1833 if ( ( ok1 != 0 && ok2 != 0 && ok3 != 0 ) || ( ok1 == 0 && ok2 != 0 && ok3 != 0 ) ||
1834 ( ok1 != 0 && ok2 == 0 && ok3 != 0 ) || ( ok1 != 0 && ok2 != 0 && ok3 == 0 ) ||
1835 ( ok1 == 0 && ok2 != 0 && ok3 == 0 ) )
1836 {
1837 //...deletes
1838 deleteListForCurl( arcZList1, arcZList2, arcZList3 );
1839 return -1;
1840 }
1841
1842 if ( ok1 == 0 && ok2 == 0 && ok3 == 0 )
1843 {
1844 //...finds
1845 double checkArc;
1846 int n1 = arcZList1.length();
1847 int n2 = arcZList2.length();
1848 int n3 = arcZList3.length();
1849 for ( int i = 0; i < n1; i++ )
1850 {
1851 for ( int j = 0; j < n2; j++ )
1852 {
1853 for ( int k = 0; k < n3; k++ )
1854 {
1855 AList<HepPoint3D> arcZ;
1856 arcZ.append( *arcZList1[i] );
1857 arcZ.append( *arcZList2[j] );
1858 arcZ.append( *arcZList3[k] );
1859 arcZ.sort( TArcOrder );
1860 double z01 = fabs( arcZ[0]->y() - arcZ[1]->y() );
1861 double z12 = fabs( arcZ[1]->y() - arcZ[2]->y() );
1862 if ( z01 <= minZ1 && z12 <= minZ2 )
1863 {
1864 minZ1 = z01;
1865 minZ2 = z12;
1866 flag1 = i;
1867 flag2 = j;
1868 flag3 = k;
1869 }
1870 }
1871 }
1872 }
1873 if ( minZ1 == 9999999. || minZ2 == 9999999. )
1874 {
1875 deleteListForCurl( arcZList1, arcZList2, arcZList3 );
1876 return -1;
1877 }
1878 //...sets the best values
1879 HepPoint3D tmp1 = *arcZList1[flag1];
1880 HepPoint3D tmp2 = *arcZList2[flag2];
1881 HepPoint3D tmp3 = *arcZList3[flag3];
1882 link1.position( tmp1 );
1883 link2.position( tmp2 );
1884 link3.position( tmp3 );
1885 deleteListForCurl( arcZList1, arcZList2, arcZList3 );
1886 return 0;
1887 }
1888 else if ( ok1 == 0 && ok2 == 0 && ok3 != 0 )
1889 {
1890 int n1 = arcZList1.length();
1891 int n2 = arcZList2.length();
1892 for ( int i = 0; i < n1; i++ )
1893 {
1894 for ( int j = 0; j < n2; j++ )
1895 {
1896 if ( fabs( arcZList1[i]->y() - arcZList2[j]->y() ) < minZ01 )
1897 {
1898 minZ01 = fabs( arcZList1[i]->y() - arcZList2[j]->y() );
1899 flag1 = i;
1900 flag2 = j;
1901 }
1902 }
1903 }
1904 if ( minZ01 == 9999999. )
1905 {
1906 deleteListForCurl( arcZList1, arcZList2, arcZList3 );
1907 return -1;
1908 }
1909 //...sets the best values
1910 HepPoint3D tmp1 = *arcZList1[flag1];
1911 HepPoint3D tmp2 = *arcZList2[flag2];
1912 link1.position( tmp1 );
1913 link2.position( tmp2 );
1914 deleteListForCurl( arcZList1, arcZList2, arcZList3 );
1915 return 12;
1916 }
1917 else if ( ok1 != 0 && ok2 == 0 && ok3 == 0 )
1918 {
1919 int n2 = arcZList2.length();
1920 int n3 = arcZList3.length();
1921 for ( int i = 0; i < n2; i++ )
1922 {
1923 for ( int j = 0; j < n3; j++ )
1924 {
1925 if ( fabs( arcZList2[i]->y() - arcZList3[j]->y() ) < minZ12 )
1926 {
1927 minZ12 = fabs( arcZList2[i]->y() - arcZList3[j]->y() );
1928 flag2 = i;
1929 flag3 = j;
1930 }
1931 }
1932 }
1933 if ( minZ12 == 9999999. )
1934 {
1935 deleteListForCurl( arcZList1, arcZList2, arcZList3 );
1936 return -1;
1937 }
1938 //...sets the best values
1939 HepPoint3D tmp2 = *arcZList2[flag2];
1940 HepPoint3D tmp3 = *arcZList3[flag3];
1941 link1.position( tmp2 );
1942 link2.position( tmp3 );
1943 deleteListForCurl( arcZList1, arcZList2, arcZList3 );
1944 return 23;
1945 }
1946 else
1947 {
1948 deleteListForCurl( arcZList1, arcZList2, arcZList3 );
1949 return -1;
1950 }
1951}
1952
1953void TTrack::deleteListForCurl( AList<HepPoint3D>& l1, AList<HepPoint3D>& l2 ) const {
1954 HepAListDeleteAll( l1 );
1955 HepAListDeleteAll( l2 );
1956}
1957
1958void TTrack::deleteListForCurl( AList<HepPoint3D>& l1, AList<HepPoint3D>& l2,
1959 AList<HepPoint3D>& l3 ) const {
1960 HepAListDeleteAll( l1 );
1961 HepAListDeleteAll( l2 );
1962 HepAListDeleteAll( l3 );
1963}
1964#endif // OLD_STEREO
1965
1966int TTrack::stereoHitForCurl( AList<TMLink>& list ) const {
1967 if ( list.length() == 0 ) return -1;
1968
1969 HepPoint3D center = _helix->center();
1970 HepPoint3D tmp( -999., -999., 0. );
1971 double r = fabs( _helix->curv() );
1972 for ( unsigned i = 0, size = list.length(); i < size; ++i )
1973 {
1974 TMDCWireHit& h = *const_cast<TMDCWireHit*>( list[i]->hit() );
1975 HepVector3D X = 0.5 * ( h.wire()->forwardPosition() + h.wire()->backwardPosition() );
1976 HepVector3D x = HepVector3D( X.x(), X.y(), 0. );
1977 HepVector3D w = x - center;
1978 HepVector3D V = h.wire()->direction();
1979 HepVector3D v = HepVector3D( V.x(), V.y(), 0. );
1980 double vmag2 = v.mag2();
1981 double vmag = sqrt( vmag2 );
1982 //...temporary
1983 for ( unsigned j = 0; j < 4; ++j ) list[i]->arcZ( tmp, j );
1984
1985 //...stereo?
1986 if ( vmag == 0. ) continue;
1987
1988 double drift = h.drift( WireHitLeft );
1989 double R[2] = { r + drift, r - drift };
1990 double wv = w.dot( v );
1991 double d2[2];
1992 d2[0] = vmag2 * R[0] * R[0] +
1993 ( wv * wv - vmag2 * w.mag2() ); //...= v^2*(r^2 - w^2*sin()^2)...outer
1994 d2[1] = vmag2 * R[1] * R[1] +
1995 ( wv * wv - vmag2 * w.mag2() ); //...= v^2*(r^2 - w^2*sin()^2)...inner
1996
1997 //...No crossing in R/Phi plane...
1998 if ( d2[0] < 0. && d2[1] < 0. ) continue;
1999
2000 bool ok_inner( true );
2001 bool ok_outer( true );
2002 double d[2] = { -1., -1. };
2003 //...outer
2004 if ( d2[0] >= 0. ) { d[0] = sqrt( d2[0] ); }
2005 else { ok_outer = false; }
2006 if ( d2[1] >= 0. ) { d[1] = sqrt( d2[1] ); }
2007 else { ok_inner = false; }
2008
2009 //...Cal. length and z to crossing points...
2010 double l[2][2];
2011 double z[2][2];
2012 for (int i=0; i<2; i++){
2013 for (int j=0; j<2; j++){
2014 z[i][j] = 0.; // add initialize 25-05-15
2015 }
2016 }
2017 //...outer
2018 if ( ok_outer )
2019 {
2020 l[0][0] = ( -wv + d[0] ) /
2021 vmag2; //...= (-wvcos()+d)/v/v = (-wcos() + (r^2 - w^2*sin()^2)^0.5)/v
2022 l[1][0] = ( -wv - d[0] ) /
2023 vmag2; //...= (-wvcos()+d)/v/v = (-wcos() - (r^2 - w^2*sin()^2)^0.5)/v
2024 z[0][0] = X.z() + l[0][0] * V.z();
2025 z[1][0] = X.z() + l[1][0] * V.z();
2026 }
2027 //...inner
2028 if ( ok_inner )
2029 {
2030 l[0][1] = ( -wv + d[1] ) /
2031 vmag2; //...= (-wvcos()+d)/v/v = (-wcos() + (r^2 - w^2*sin()^2)^0.5)/v
2032 l[1][1] = ( -wv - d[1] ) /
2033 vmag2; //...= (-wvcos()+d)/v/v = (-wcos() - (r^2 - w^2*sin()^2)^0.5)/v
2034 z[0][1] = X.z() + l[0][1] * V.z();
2035 z[1][1] = X.z() + l[1][1] * V.z();
2036 }
2037
2038 //...Cal. xy position of crossing points...
2039 HepVector3D p[2][2];
2040#if 1
2041 HepVector3D tp[2][2];
2042#endif
2043 if ( ok_outer )
2044 {
2045 p[0][0] = x + l[0][0] * v;
2046 p[1][0] = x + l[1][0] * v;
2047#if 1
2048 tp[0][0] = p[0][0];
2049 tp[1][0] = p[1][0];
2050#endif
2051 HepVector3D tmp_pc = p[0][0] - center;
2052 HepVector3D pc0 = HepVector3D( tmp_pc.x(), tmp_pc.y(), 0. );
2053 p[0][0] -= drift / pc0.mag() * pc0;
2054 tmp_pc = p[1][0] - center;
2055 HepVector3D pc1 = HepVector3D( tmp_pc.x(), tmp_pc.y(), 0. );
2056 p[1][0] -= drift / pc1.mag() * pc1;
2057 }
2058#define JJTEST 0
2059 if ( ok_inner )
2060 {
2061 p[0][1] = x + l[0][1] * v;
2062 p[1][1] = x + l[1][1] * v;
2063#if JJTEST
2064 tp[0][1] = p[0][1];
2065 tp[1][1] = p[1][1];
2066#endif
2067 HepVector3D tmp_pc = p[0][1] - center;
2068 HepVector3D pc0 = HepVector3D( tmp_pc.x(), tmp_pc.y(), 0. );
2069 p[0][1] += drift / pc0.mag() * pc0;
2070 tmp_pc = p[1][1] - center;
2071 HepVector3D pc1 = HepVector3D( tmp_pc.x(), tmp_pc.y(), 0. );
2072 p[1][1] += drift / pc1.mag() * pc1;
2073 }
2074
2075 //...boolean
2076 bool ok_xy[2][2];
2077 if ( ok_outer )
2078 {
2079 ok_xy[0][0] = true;
2080 ok_xy[1][0] = true;
2081 }
2082 else
2083 {
2084 ok_xy[0][0] = false;
2085 ok_xy[1][0] = false;
2086 }
2087 if ( ok_inner )
2088 {
2089 ok_xy[0][1] = true;
2090 ok_xy[1][1] = true;
2091 }
2092 else
2093 {
2094 ok_xy[0][1] = false;
2095 ok_xy[1][1] = false;
2096 }
2097 if ( ok_outer )
2098 {
2099 if ( _charge * ( center.x() * p[0][0].y() - center.y() * p[0][0].x() ) < 0. )
2100 ok_xy[0][0] = false;
2101 if ( _charge * ( center.x() * p[1][0].y() - center.y() * p[1][0].x() ) < 0. )
2102 ok_xy[1][0] = false;
2103 }
2104 if ( ok_inner )
2105 {
2106 if ( _charge * ( center.x() * p[0][1].y() - center.y() * p[0][1].x() ) < 0. )
2107 ok_xy[0][1] = false;
2108 if ( _charge * ( center.x() * p[1][1].y() - center.y() * p[1][1].x() ) < 0. )
2109 ok_xy[1][1] = false;
2110 }
2111 if ( !ok_inner && ok_outer && ( !ok_xy[0][0] ) && ( !ok_xy[1][0] ) ) { continue; }
2112 if ( ok_inner && !ok_outer && ( !ok_xy[0][1] ) && ( !ok_xy[1][1] ) ) { continue; }
2113
2114 //...Check z position...
2115 if ( ok_xy[0][0] )
2116 {
2117 if ( z[0][0] < h.wire()->backwardPosition().z() ||
2118 z[0][0] > h.wire()->forwardPosition().z() )
2119 ok_xy[0][0] = false;
2120 }
2121 if ( ok_xy[1][0] )
2122 {
2123 if ( z[1][0] < h.wire()->backwardPosition().z() ||
2124 z[1][0] > h.wire()->forwardPosition().z() )
2125 ok_xy[1][0] = false;
2126 }
2127 if ( ok_xy[0][1] )
2128 {
2129 if ( z[0][1] < h.wire()->backwardPosition().z() ||
2130 z[0][1] > h.wire()->forwardPosition().z() )
2131 ok_xy[0][1] = false;
2132 }
2133 if ( ok_xy[1][1] )
2134 {
2135 if ( z[1][1] < h.wire()->backwardPosition().z() ||
2136 z[1][1] > h.wire()->forwardPosition().z() )
2137 ok_xy[1][1] = false;
2138 }
2139 if ( ( !ok_xy[0][0] ) && ( !ok_xy[1][0] ) && ( !ok_xy[0][1] ) && ( !ok_xy[1][1] ) )
2140 { continue; }
2141 double cosdPhi, dPhi;
2142 unsigned index;
2143 index = 0;
2144 if ( ok_xy[0][0] )
2145 {
2146 //...cal. arc length...
2147 cosdPhi = -center.dot( ( p[0][0] - center ).unit() ) / center.mag();
2148 if ( fabs( cosdPhi ) < 1.0 ) { dPhi = acos( cosdPhi ); }
2149 else if ( cosdPhi >= 1.0 ) { dPhi = 0.; }
2150 else { dPhi = M_PI; }
2151 list[i]->arcZ( HepPoint3D( r * dPhi, z[0][0], 0. ), index );
2152 // std::cout << r*dPhi << ", " << z[0][0] << std::endl;
2153 ++index;
2154 }
2155 if ( ok_xy[1][0] )
2156 {
2157 //...cal. arc length...
2158 cosdPhi = -center.dot( ( p[1][0] - center ).unit() ) / center.mag();
2159 if ( fabs( cosdPhi ) < 1.0 ) { dPhi = acos( cosdPhi ); }
2160 else if ( cosdPhi >= 1.0 ) { dPhi = 0.; }
2161 else { dPhi = M_PI; }
2162 list[i]->arcZ( HepPoint3D( r * dPhi, z[1][0], 0. ), index );
2163 // std::cout << r*dPhi << ", " << z[1][0] << std::endl;
2164 ++index;
2165 }
2166 if ( ok_xy[0][1] )
2167 {
2168 //...cal. arc length...
2169 cosdPhi = -center.dot( ( p[0][1] - center ).unit() ) / center.mag();
2170 if ( fabs( cosdPhi ) < 1.0 ) { dPhi = acos( cosdPhi ); }
2171 else if ( cosdPhi >= 1.0 ) { dPhi = 0.; }
2172 else { dPhi = M_PI; }
2173 list[i]->arcZ( HepPoint3D( r * dPhi, z[0][1], 0. ), index );
2174 // std::cout << r*dPhi << ", " << z[0][1] << std::endl;
2175 ++index;
2176 }
2177 if ( ok_xy[1][1] )
2178 {
2179 //...cal. arc length...
2180 cosdPhi = -center.dot( ( p[1][1] - center ).unit() ) / center.mag();
2181 if ( fabs( cosdPhi ) < 1.0 ) { dPhi = acos( cosdPhi ); }
2182 else if ( cosdPhi >= 1.0 ) { dPhi = 0.; }
2183 else { dPhi = M_PI; }
2184 list[i]->arcZ( HepPoint3D( r * dPhi, z[1][1], 0. ), index );
2185 // std::cout << r*dPhi << ", " << z[1][1] << std::endl;
2186 ++index;
2187 }
2188
2189#if JJTEST
2190 double gmaxX = -9999., gminX = 9999.;
2191 double gmaxY = -9999., gminY = 9999.;
2192 FILE * gnuplot, *data;
2193 double step = 100.;
2194 double dStep = 2. * M_PI / step;
2195 if ( ( data = fopen( "dat1", "w" ) ) != NULL )
2196 {
2197 if ( ok_xy[0][0] )
2198 {
2199 for ( int ii = 0; ii < step; ++ii )
2200 {
2201 double X = tp[0][0].x() + drift * cos( dStep * static_cast<double>( ii ) );
2202 double Y = tp[0][0].y() + drift * sin( dStep * static_cast<double>( ii ) );
2203 fprintf( data, "%lf, %lf\n", X, Y );
2204 if ( gmaxX < X ) gmaxX = X;
2205 if ( gminX > X ) gminX = X;
2206 if ( gmaxY < Y ) gmaxY = Y;
2207 if ( gminY > Y ) gminY = Y;
2208 }
2209 }
2210 fclose( data );
2211 }
2212 if ( ( data = fopen( "dat2", "w" ) ) != NULL )
2213 {
2214 if ( ok_xy[1][0] )
2215 {
2216 for ( int ii = 0; ii < step; ++ii )
2217 {
2218 double X = tp[1][0].x() + drift * cos( dStep * static_cast<double>( ii ) );
2219 double Y = tp[1][0].y() + drift * sin( dStep * static_cast<double>( ii ) );
2220 fprintf( data, "%lf, %lf\n", X, Y );
2221 if ( gmaxX < X ) gmaxX = X;
2222 if ( gminX > X ) gminX = X;
2223 if ( gmaxY < Y ) gmaxY = Y;
2224 if ( gminY > Y ) gminY = Y;
2225 }
2226 }
2227 fclose( data );
2228 }
2229 if ( ( data = fopen( "dat3", "w" ) ) != NULL )
2230 {
2231 if ( ok_xy[0][1] )
2232 {
2233 for ( int ii = 0; ii < step; ++ii )
2234 {
2235 double X = tp[0][1].x() + drift * cos( dStep * static_cast<double>( ii ) );
2236 double Y = tp[0][1].y() + drift * sin( dStep * static_cast<double>( ii ) );
2237 fprintf( data, "%lf, %lf\n", X, Y );
2238 if ( gmaxX < X ) gmaxX = X;
2239 if ( gminX > X ) gminX = X;
2240 if ( gmaxY < Y ) gmaxY = Y;
2241 if ( gminY > Y ) gminY = Y;
2242 }
2243 }
2244 fclose( data );
2245 }
2246 if ( ( data = fopen( "dat4", "w" ) ) != NULL )
2247 {
2248 if ( ok_xy[1][1] )
2249 {
2250 for ( int ii = 0; ii < step; ++ii )
2251 {
2252 double X = tp[1][1].x() + drift * cos( dStep * static_cast<double>( ii ) );
2253 double Y = tp[1][1].y() + drift * sin( dStep * static_cast<double>( ii ) );
2254 fprintf( data, "%lf, %lf\n", X, Y );
2255 if ( gmaxX < X ) gmaxX = X;
2256 if ( gminX > X ) gminX = X;
2257 if ( gmaxY < Y ) gmaxY = Y;
2258 if ( gminY > Y ) gminY = Y;
2259 }
2260 }
2261 fclose( data );
2262 }
2263 HepVector3D tX = h.wire()->forwardPosition() - h.wire()->backwardPosition();
2264 HepVector3D tDist( tX.x(), tX.y(), 0. );
2265 double tD = tDist.mag();
2266 double vvvM = 1. / v.mag();
2267 HepVector3D tDire = vvvM * v;
2268 step = 2.;
2269 dStep = tD / step;
2270 if ( ( data = fopen( "dat5", "w" ) ) != NULL )
2271 {
2272 for ( int ii = 0; ii < step + 1; ++ii )
2273 {
2274 double X =
2275 h.wire()->backwardPosition().x() + dStep * static_cast<double>( ii ) * tDire.x();
2276 double Y =
2277 h.wire()->backwardPosition().y() + dStep * static_cast<double>( ii ) * tDire.y();
2278 fprintf( data, "%lf, %lf\n", X, Y );
2279 if ( gmaxX < X ) gmaxX = X;
2280 if ( gminX > X ) gminX = X;
2281 if ( gmaxY < Y ) gmaxY = Y;
2282 if ( gminY > Y ) gminY = Y;
2283 }
2284 fclose( data );
2285 }
2286 if ( ( data = fopen( "dat6", "w" ) ) != NULL )
2287 {
2288 double X = h.wire()->backwardPosition().x();
2289 double Y = h.wire()->backwardPosition().y();
2290 fprintf( data, "%lf, %lf\n", X, Y );
2291 if ( gmaxX < X ) gmaxX = X;
2292 if ( gminX > X ) gminX = X;
2293 if ( gmaxY < Y ) gmaxY = Y;
2294 if ( gminY > Y ) gminY = Y;
2295 fclose( data );
2296 }
2297 if ( ( data = fopen( "dat7", "w" ) ) != NULL )
2298 {
2299 double X = h.wire()->forwardPosition().x();
2300 double Y = h.wire()->forwardPosition().y();
2301 fprintf( data, "%lf, %lf\n", X, Y );
2302 if ( gmaxX < X ) gmaxX = X;
2303 if ( gminX > X ) gminX = X;
2304 if ( gmaxY < Y ) gmaxY = Y;
2305 if ( gminY > Y ) gminY = Y;
2306 fclose( data );
2307 }
2308 step = 300.;
2309 dStep = 2. * M_PI / step;
2310 if ( ( data = fopen( "dat8", "w" ) ) != NULL )
2311 {
2312 for ( int ii = 0; ii < step; ++ii )
2313 {
2314 double X = center.x() + r * cos( dStep * static_cast<double>( ii ) );
2315 double Y = center.y() + r * sin( dStep * static_cast<double>( ii ) );
2316 fprintf( data, "%lf, %lf\n", X, Y );
2317 }
2318 fclose( data );
2319 }
2320 if ( ( data = fopen( "dat9", "w" ) ) != NULL )
2321 {
2322 if ( ok_xy[0][0] )
2323 {
2324 double X = p[0][0].x();
2325 double Y = p[0][0].y();
2326 fprintf( data, "%lf, %lf\n", X, Y );
2327 }
2328 fclose( data );
2329 }
2330 if ( ( data = fopen( "dat10", "w" ) ) != NULL )
2331 {
2332 if ( ok_xy[1][0] )
2333 {
2334 double X = p[1][0].x();
2335 double Y = p[1][0].y();
2336 fprintf( data, "%lf, %lf\n", X, Y );
2337 }
2338 fclose( data );
2339 }
2340 if ( ( data = fopen( "dat11", "w" ) ) != NULL )
2341 {
2342 if ( ok_xy[0][1] )
2343 {
2344 double X = p[0][1].x();
2345 double Y = p[0][1].y();
2346 fprintf( data, "%lf, %lf\n", X, Y );
2347 }
2348 fclose( data );
2349 }
2350 if ( ( data = fopen( "dat12", "w" ) ) != NULL )
2351 {
2352 if ( ok_xy[1][1] )
2353 {
2354 double X = p[1][1].x();
2355 double Y = p[1][1].y();
2356 fprintf( data, "%lf, %lf\n", X, Y );
2357 }
2358 fclose( data );
2359 }
2360 std::cout << "Drift Distance = " << drift << ", xy1cm -> z" << V.z() / v.mag()
2361 << "cm, xyDist(cm) = " << tD << std::endl;
2362 if ( tX.z() < 0. ) std::cout << "ERROR : F < B" << std::endl;
2363 if ( ( gnuplot = popen( "gnuplot", "w" ) ) != NULL )
2364 {
2365 fprintf( gnuplot, "set size 0.721,1.0 \n" );
2366 fprintf( gnuplot, "set xrange [%f:%f] \n", gminX, gmaxX );
2367 fprintf( gnuplot, "set yrange [%f:%f] \n", gminY, gmaxY );
2368 fprintf( gnuplot, "plot \"dat1\" with line, \"dat2\" with line, \"dat3\" with line, "
2369 "\"dat4\" with line, \"dat5\" with "
2370 "line, \"dat6\", \"dat7\", \"dat8\" with line, \"dat9\", \"dat10\", "
2371 "\"dat11\", \"dat12\" \n" );
2372 fflush( gnuplot );
2373 char tmp[8];
2374 gets( tmp );
2375 pclose( gnuplot );
2376 }
2377#endif
2378 }
2379 return 0;
2380}
2381
2382#if !( NEW_FIT2D )
2383int TTrack::approach2D( TMLink& l ) const {
2384
2385 //...Setup...
2386 const TMDCWire& w = *l.wire();
2387 Vector a = _helix->a();
2388 double kappa = a[2];
2389 double phi0 = a[1];
2390 HepPoint3D xc = _helix->center();
2391 HepPoint3D xw = w.xyPosition();
2392 HepPoint3D xt = _helix->x();
2393 HepVector3D v0, v1;
2394 v0 = xt - xc;
2395 v1 = xw - xc;
2396 // if (_charge > 0.) {
2397 // v0 *= -1;
2398 // v1 *= -1;
2399 // }
2400 double vCrs = v0.x() * v1.y() - v0.y() * v1.x();
2401 double vDot = v0.x() * v1.x() + v0.y() * v1.y();
2402 double dPhi = atan2( vCrs, vDot );
2403
2404 //...Cal. phi to rotate...
2405 // double phiWire = atan2(_charge * (xc.y() - xw.y()),
2406 // _charge * (xc.x() - xw.x()));
2407 // phiWire = (phiWire > 0.) ? phiWire : phiWire + 2. * M_PI;
2408 // double dPhi = phiWire - phi0;
2409
2410 //...Why this is needed?...
2411 // if ((_charge >= 0.) && (dPhi > 0.)) dPhi -= 2. * M_PI;
2412 // else if ((_charge < 0.) && (dPhi < 0.)) dPhi += 2. * M_PI;
2413
2414 // std::cout << _helix->x(dPhi) << " , " << _helix->x(dPhi+0.2) << " , " <<
2415 // _helix->x(dPhi-0.1) << std::endl;
2416
2417 l.positionOnTrack( _helix->x( dPhi ) );
2418 HepPoint3D x = xw;
2419 x.setZ( 0. );
2420 l.positionOnWire( x );
2421 l.dPhi( dPhi );
2422 return 0;
2423}
2424
2425int TTrack::dxda2D( const TMLink& link, double dPhi, Vector& dxda, Vector& dyda,
2426 Vector& dzda ) const {
2427
2428 //...Setup...
2429 const TMDCWire& w = *link.wire();
2430 Vector a = _helix->a();
2431 double dRho = a[0];
2432 double phi0 = a[1];
2433 double kappa = a[2];
2434 // double rho = Helix::ConstantAlpha / kappa;
2435 // double rho = 333.564095 / kappa;
2436 const double Bz = -1000 * getPmgnIMF()->getReferField();
2437 double rho = 333.564095 / ( kappa * Bz );
2438
2439 // double tanLambda = a[4];
2440
2441 double sinPhi0 = sin( phi0 );
2442 double cosPhi0 = cos( phi0 );
2443 double sinPhi0dPhi = sin( phi0 + dPhi );
2444 double cosPhi0dPhi = cos( phi0 + dPhi );
2445 Vector dphida( 5 );
2446
2447 HepPoint3D d = _helix->center() - w.xyPosition();
2448 double dmag2 = d.mag2();
2449
2450 dphida[0] = ( sinPhi0 * d.x() - cosPhi0 * d.y() ) / dmag2;
2451 dphida[1] = ( dRho + rho ) * ( cosPhi0 * d.x() + sinPhi0 * d.y() ) / dmag2 - 1.;
2452 dphida[2] = ( -rho / kappa ) * ( sinPhi0 * d.x() - cosPhi0 * d.y() ) / dmag2;
2453 dphida[3] = 0.;
2454 dphida[4] = 0.;
2455
2456 dxda[0] = cosPhi0 + rho * sinPhi0dPhi * dphida[0];
2457 dxda[1] = -( dRho + rho ) * sinPhi0 + rho * sinPhi0dPhi * ( 1. + dphida[1] );
2458 dxda[2] = -rho / kappa * ( cosPhi0 - cosPhi0dPhi ) + rho * sinPhi0dPhi * dphida[2];
2459 dxda[3] = rho * sinPhi0dPhi * dphida[3];
2460 dxda[4] = rho * sinPhi0dPhi * dphida[4];
2461
2462 dyda[0] = sinPhi0 - rho * cosPhi0dPhi * dphida[0];
2463 dyda[1] = ( dRho + rho ) * cosPhi0 - rho * cosPhi0dPhi * ( 1. + dphida[1] );
2464 dyda[2] = -rho / kappa * ( sinPhi0 - sinPhi0dPhi ) - rho * cosPhi0dPhi * dphida[2];
2465 dyda[3] = -rho * cosPhi0dPhi * dphida[3];
2466 dyda[4] = -rho * cosPhi0dPhi * dphida[4];
2467
2468 dzda[0] = 0.;
2469 dzda[1] = 0.;
2470 dzda[2] = 0.;
2471 dzda[3] = 1.;
2472 dzda[4] = -rho * dPhi;
2473 // dzda[0] = - rho * tanLambda * dphida[0];
2474 // dzda[1] = - rho * tanLambda * dphida[1];
2475 // dzda[2] = rho / kappa * tanLambda * dPhi - rho * tanLambda * dphida[2];
2476 // dzda[3] = 1. - rho * tanLambda * dphida[3];
2477 // dzda[4] = - rho * dPhi - rho * tanLambda * dphida[4];
2478
2479 return 0;
2480}
2481#endif
2482#if 0
2483int
2484TTrack::svdHitForCurl(AList<TSvdHit>& svdList) const
2485{
2486 HepPoint3D center = _helix->center();
2487 double r = fabs(_helix->curv());
2488
2489 for(unsigned i = 0, size = svdList.length(); i < size; ++i){
2490 HepPoint3D p(svdList[i]->x(), svdList[i]->y(), 0.);
2491 double cosdPhi = - center.dot((p - center).unit()) / center.mag();
2492 double dPhi;
2493 if(fabs(cosdPhi) < 1.0){
2494 dPhi = acos(cosdPhi);
2495 }else if(cosdPhi >= 1.0){
2496 dPhi = 0.;
2497 }else{
2498 dPhi = M_PI;
2499 }
2500 HepPoint3D tmp(r*dPhi, svdList[i]->z(), 0.);
2501 svdList[i]->arcZ(tmp);
2502 }
2503 return 0;
2504}
2505#endif
2506
2507#if defined( __GNUG__ )
2508int SortByPt( const TTrack** a, const TTrack** b ) {
2509 if ( ( *a )->pt() < ( *b )->pt() ) return 1;
2510 else if ( ( *a )->pt() == ( *b )->pt() ) return 0;
2511 else return -1;
2512}
2513#else
2514extern "C" int SortByPt( const void* av, const void* bv ) {
2515 const TTrack** a( (const TTrack**)av );
2516 const TTrack** b( (const TTrack**)bv );
2517 if ( ( *a )->pt() < ( *b )->pt() ) return 1;
2518 else if ( ( *a )->pt() == ( *b )->pt() ) return 0;
2519 else return -1;
2520}
2521#endif
2522
2523#if NEW_FIT2D
2524int TTrack::fit2D( unsigned ipFlag, double ipDistance, double ipError ) {
2525# ifdef TRKRECO_DEBUG_DETAIL
2526 std::cout << " TTrack::fit2D(r-phi) ..." << std::endl;
2527# endif
2528 // if(_fitted)return 0;
2529
2530 //...Check # of hits...
2531
2532 // std::cout << "# = " << _links.length() << std::endl;
2533 //...Setup...
2534 unsigned nTrial( 0 );
2535 Vector a( _helix->a() );
2536 double chi2;
2537 double chi2Old = 1.0e+99;
2538 Vector dchi2da( 3 );
2539 SymMatrix d2chi2d2a( 3, 0 );
2540 Vector da( 5 ), dxda( 3 ), dyda( 3 );
2541 Vector dDda( 3 );
2542 const double convergence = 1.0e-5;
2543 Vector f( 3 );
2544 int err = 0;
2545 double factor = 1.0;
2546 unsigned usedWireNumber = 0;
2547
2548 //...Fitting loop...
2549 while ( nTrial < 100 )
2550 {
2551 //...Set up...
2552 chi2 = 0.;
2553 for ( unsigned j = 0; j < 3; ++j ) dchi2da[j] = 0.;
2554 d2chi2d2a = SymMatrix( 3, 0 );
2555 usedWireNumber = 0;
2556
2557 //...Loop with hits...
2558 unsigned i = 0;
2559 while ( TMLink* l = _links[i++] )
2560 {
2561 if ( l->fit2D() == 0 ) continue;
2562 const TMDCWireHit& h = *l->hit();
2563
2564 //...Cal. closest points...
2565 if ( approach2D( *l ) != 0 ) continue;
2566 double dPhi = l->dPhi();
2567 const HepPoint3D& onTrack = l->positionOnTrack();
2568 const HepPoint3D& onWire = l->positionOnWire();
2569 HepPoint3D onTrack2( onTrack.x(), onTrack.y(), 0. );
2570 HepPoint3D onWire2( onWire.x(), onWire.y(), 0. );
2571
2572 //...Obtain drift distance and its error...
2573 unsigned leftRight = WireHitRight;
2574 if ( onWire2.cross( onTrack2 ).z() < 0. ) leftRight = WireHitLeft;
2575 double distance = h.drift( leftRight );
2576 double eDistance = h.dDrift( leftRight );
2577 double eDistance2 = eDistance * eDistance;
2578 if ( eDistance == 0. )
2579 {
2580 std::cout << "Error(?) : Drift Distance Error = 0 in fit2D of TrkReco." << std::endl;
2581 continue;
2582 }
2583
2584 //...Residual...
2585 HepVector3D v = onTrack2 - onWire2;
2586 double vmag = v.mag();
2587 double dDistance = vmag - distance;
2588
2589 //...dxda...
2590 if ( this->dxda2D( *l, dPhi, dxda, dyda ) != 0 ) continue;
2591
2592 //...Chi2 related...
2593 // Vector3 vw(0.,0.,1.);
2594 dDda = ( vmag > 0. ) ? ( v.x() * dxda + v.y() * dyda ) / vmag : Vector( 3, 0 );
2595 if ( vmag <= 0.0 )
2596 {
2597 std::cout << " in fit2D " << onTrack << ", " << onWire;
2598 h.dump();
2599 continue;
2600 }
2601 dchi2da += ( dDistance / eDistance2 ) * dDda;
2602 d2chi2d2a += vT_times_v( dDda ) / eDistance2;
2603 double pChi2 = dDistance * dDistance / eDistance2;
2604 chi2 += pChi2;
2605
2606 //...Store results...
2607 l->update( onTrack2, onWire2, leftRight, pChi2 );
2608 ++usedWireNumber;
2609 }
2610 if ( ipFlag != 0 )
2611 {
2612 double kappa = _helix->a()[2];
2613 double phi0 = _helix->a()[1];
2614 HepPoint3D xc( _helix->center() );
2615 HepPoint3D onWire( 0., 0., 0. );
2616 HepPoint3D xt( _helix->x() );
2617 HepVector3D v0( xt - xc );
2618 HepVector3D v1( onWire - xc );
2619 double vCrs = v0.x() * v1.y() - v0.y() * v1.x();
2620 double vDot = v0.x() * v1.x() + v0.y() * v1.y();
2621 double dPhi = atan2( vCrs, vDot );
2622 HepPoint3D onTrack( _helix->x( dPhi ).x(), _helix->x( dPhi ).y(), 0. );
2623 double distance = ipDistance;
2624 double eDistance = ipError;
2625 double eDistance2 = eDistance * eDistance;
2626
2627 HepVector3D v = onTrack - onWire;
2628 double vmag = v.mag();
2629 double dDistance = vmag - distance;
2630
2631 if ( this->dxda2D( dPhi, dxda, dyda ) != 0 ) goto ipOff;
2632
2633 dDda = ( vmag > 0. ) ? ( v.x() * dxda + v.y() * dyda ) / vmag : Vector( 3, 0 );
2634 if ( vmag <= 0.0 ) { goto ipOff; }
2635 dchi2da += ( dDistance / eDistance2 ) * dDda;
2636 d2chi2d2a += vT_times_v( dDda ) / eDistance2;
2637 double pChi2 = dDistance * dDistance / eDistance2;
2638 chi2 += pChi2;
2639
2640 ++usedWireNumber;
2641 }
2642 ipOff:
2643 if ( usedWireNumber < 4 )
2644 {
2645 err = -2;
2646 break;
2647 }
2648
2649 //...Check condition...
2650 double change = chi2Old - chi2;
2651 if ( fabs( change ) < convergence ) break;
2652 if ( change < 0. )
2653 {
2654# ifdef TRKRECO_DEBUG_DETAIL
2655 std::cout << "chi2Old, chi2=" << chi2Old << " " << chi2 << std::endl;
2656# endif
2657 // change to the old value.
2658 a += factor * da;
2659 _helix->a( a );
2660
2661 chi2 = 0.;
2662 for ( unsigned j = 0; j < 3; ++j ) dchi2da[j] = 0.;
2663 d2chi2d2a = SymMatrix( 3, 0 );
2664 usedWireNumber = 0;
2665
2666 //...Loop with hits...
2667 unsigned i = 0;
2668 while ( TMLink* l = _links[i++] )
2669 {
2670 if ( l->fit2D() == 0 ) continue;
2671 const TMDCWireHit& h = *l->hit();
2672
2673 //...Cal. closest points...
2674 if ( approach2D( *l ) != 0 ) continue;
2675 double dPhi = l->dPhi();
2676 const HepPoint3D& onTrack = l->positionOnTrack();
2677 const HepPoint3D& onWire = l->positionOnWire();
2678 HepPoint3D onTrack2( onTrack.x(), onTrack.y(), 0. );
2679 HepPoint3D onWire2( onWire.x(), onWire.y(), 0. );
2680
2681 //...Obtain drift distance and its error...
2682 unsigned leftRight = WireHitRight;
2683 if ( onWire2.cross( onTrack2 ).z() < 0. ) leftRight = WireHitLeft;
2684 double distance = h.drift( leftRight );
2685 double eDistance = h.dDrift( leftRight );
2686 double eDistance2 = eDistance * eDistance;
2687
2688 //...Residual...
2689 HepVector3D v = onTrack2 - onWire2;
2690 double vmag = v.mag();
2691 double dDistance = vmag - distance;
2692
2693 //...dxda...
2694 if ( this->dxda2D( *l, dPhi, dxda, dyda ) != 0 ) continue;
2695
2696 //...Chi2 related...
2697 dDda = ( vmag > 0. ) ? ( v.x() * dxda + v.y() * dyda ) / vmag : Vector( 3, 0 );
2698 if ( vmag <= 0.0 )
2699 {
2700 std::cout << " in fit2D " << onTrack << ", " << onWire;
2701 h.dump();
2702 continue;
2703 }
2704 dchi2da += ( dDistance / eDistance2 ) * dDda;
2705 d2chi2d2a += vT_times_v( dDda ) / eDistance2;
2706 double pChi2 = dDistance * dDistance / eDistance2;
2707 chi2 += pChi2;
2708
2709 //...Store results...
2710 l->update( onTrack2, onWire2, leftRight, pChi2 );
2711 ++usedWireNumber;
2712 }
2713 if ( ipFlag != 0 )
2714 {
2715 double kappa = _helix->a()[2];
2716 double phi0 = _helix->a()[1];
2717 HepPoint3D xc( _helix->center() );
2718 HepPoint3D onWire( 0., 0., 0. );
2719 HepPoint3D xt( _helix->x() );
2720 HepVector3D v0( xt - xc );
2721 HepVector3D v1( onWire - xc );
2722 double vCrs = v0.x() * v1.y() - v0.y() * v1.x();
2723 double vDot = v0.x() * v1.x() + v0.y() * v1.y();
2724 double dPhi = atan2( vCrs, vDot );
2725 HepPoint3D onTrack( _helix->x( dPhi ).x(), _helix->x( dPhi ).y(), 0. );
2726 double distance = ipDistance;
2727 double eDistance = ipError;
2728 double eDistance2 = eDistance * eDistance;
2729
2730 HepVector3D v = onTrack - onWire;
2731 double vmag = v.mag();
2732 double dDistance = vmag - distance;
2733
2734 if ( this->dxda2D( dPhi, dxda, dyda ) != 0 ) goto ipOff2;
2735
2736 dDda = ( vmag > 0. ) ? ( v.x() * dxda + v.y() * dyda ) / vmag : Vector( 3, 0 );
2737 if ( vmag <= 0.0 ) { goto ipOff2; }
2738 dchi2da += ( dDistance / eDistance2 ) * dDda;
2739 d2chi2d2a += vT_times_v( dDda ) / eDistance2;
2740 double pChi2 = dDistance * dDistance / eDistance2;
2741 chi2 += pChi2;
2742
2743 ++usedWireNumber;
2744 }
2745 ipOff2:
2746 if ( usedWireNumber < 4 )
2747 {
2748 err = -2;
2749 break;
2750 }
2751 // break;
2752 factor *= 0.75;
2753# ifdef TRKRECO_DEBUG_DETAIL
2754 std::cout << "factor = " << factor << std::endl;
2755 std::cout << "chi2 = " << chi2 << std::endl;
2756# endif
2757 if ( factor < 0.01 ) break;
2758 }
2759
2760 chi2Old = chi2;
2761
2762 //...Cal. helix parameters for next loop...
2763 f = solve( d2chi2d2a, dchi2da );
2764 da[0] = f[0];
2765 da[1] = f[1];
2766 da[2] = f[2];
2767 da[3] = 0.;
2768 da[4] = 0.;
2769
2770 a -= factor * da;
2771 _helix->a( a );
2772 ++nTrial;
2773 }
2774 if ( err ) { return err; }
2775
2776 //...Cal. error matrix...
2777 SymMatrix Ea( 5, 0 );
2778 unsigned dim = 3;
2779 SymMatrix Eb = d2chi2d2a;
2780 SymMatrix Ec = Eb.inverse( err );
2781 Ea[0][0] = Ec[0][0];
2782 Ea[0][1] = Ec[0][1];
2783 Ea[0][2] = Ec[0][2];
2784 Ea[1][1] = Ec[1][1];
2785 Ea[1][2] = Ec[1][2];
2786 Ea[2][2] = Ec[2][2];
2787
2788 //...Store information...
2789 if ( !err )
2790 {
2791 _helix->a( a );
2792 _helix->Ea( Ea );
2793 }
2794 else { err = -2; }
2795
2796 _ndf = usedWireNumber - dim;
2797 _chi2 = chi2;
2798
2799 //_fitted = true;
2800
2801# define JJJTEST 0
2802# if JJJTEST
2803 double gmaxX = -9999., gminX = 9999.;
2804 double gmaxY = -9999., gminY = 9999.;
2805 FILE * gnuplot, *data;
2806 double step = 200.;
2807 double dStep = 2. * M_PI / step;
2808 for ( int i = 0, size = _links.length(); i < size; ++i )
2809 {
2810 TMLink* l = _links[i];
2811 double drift = l->hit()->distance( 0 );
2812 char name[100] = "dat";
2813 char counter[10] = "";
2814 sprintf( counter, "%02d", i );
2815 strcat( name, counter );
2816 if ( ( data = fopen( name, "w" ) ) != NULL )
2817 {
2818 for ( int ii = 0; ii < step; ++ii )
2819 {
2820 double X =
2821 l->wire()->xyPosition().x() + drift * cos( dStep * static_cast<double>( ii ) );
2822 double Y =
2823 l->wire()->xyPosition().y() + drift * sin( dStep * static_cast<double>( ii ) );
2824 fprintf( data, "%lf, %lf\n", X, Y );
2825 if ( gmaxX < X ) gmaxX = X;
2826 if ( gminX > X ) gminX = X;
2827 if ( gmaxY < Y ) gmaxY = Y;
2828 if ( gminY > Y ) gminY = Y;
2829 }
2830 fclose( data );
2831 }
2832 }
2833 step = 300.;
2834 dStep = 2. * M_PI / step;
2835 if ( ( data = fopen( "datc", "w" ) ) != NULL )
2836 {
2837 for ( int ii = 0; ii < step; ++ii )
2838 {
2839 double X =
2840 _helix->center().x() + _helix->radius() * cos( dStep * static_cast<double>( ii ) );
2841 double Y =
2842 _helix->center().y() + _helix->radius() * sin( dStep * static_cast<double>( ii ) );
2843 fprintf( data, "%lf, %lf\n", X, Y );
2844 }
2845 fclose( data );
2846 }
2847 if ( ( gnuplot = popen( "gnuplot", "w" ) ) != NULL )
2848 {
2849 fprintf( gnuplot, "set size 0.721,1.0 \n" );
2850 fprintf( gnuplot, "set xrange [%f:%f] \n", gminX, gmaxX );
2851 fprintf( gnuplot, "set yrange [%f:%f] \n", gminY, gmaxY );
2852 if ( _links.length() == 4 )
2853 {
2854 fprintf( gnuplot, "plot \"datc\" with line, \"dat00\" with line, \"dat01\" with line, "
2855 "\"dat02\" with line, "
2856 "\"dat03\" with line \n" );
2857 }
2858 else if ( _links.length() == 5 )
2859 {
2860 fprintf( gnuplot, "plot \"datc\" with line, \"dat00\" with line, \"dat01\" with line, "
2861 "\"dat02\" with line, "
2862 "\"dat03\" with line, \"dat04\" with line \n" );
2863 }
2864 else if ( _links.length() == 6 )
2865 {
2866 fprintf( gnuplot, "plot \"datc\" with line, \"dat00\" with line, \"dat01\" with line, "
2867 "\"dat02\" with line, "
2868 "\"dat03\" with line, \"dat04\" with line, \"dat05\" with line \n" );
2869 }
2870 else if ( _links.length() == 7 )
2871 {
2872 fprintf( gnuplot, "plot \"datc\" with line, \"dat00\" with line, \"dat01\" with line, "
2873 "\"dat02\" with line, "
2874 "\"dat03\" with line, \"dat04\" with line, \"dat05\" with line, "
2875 "\"dat06\" with line \n" );
2876 }
2877 else if ( _links.length() == 8 )
2878 {
2879 fprintf( gnuplot, "plot \"datc\" with line, \"dat00\" with line, \"dat01\" with line, "
2880 "\"dat02\" with line, \"dat03\" with "
2881 "line, \"dat04\" with line, \"dat05\" with line, \"dat06\" with line, "
2882 "\"dat07\" with line \n" );
2883 }
2884 else if ( _links.length() == 9 )
2885 {
2886 fprintf( gnuplot, "plot \"datc\" with line, \"dat00\" with line, \"dat01\" with line, "
2887 "\"dat02\" with line, "
2888 "\"dat03\" with line, \"dat04\" with line, \"dat05\" with line, "
2889 "\"dat06\" with line, \"dat07\" "
2890 "with line, \"dat08\" with line \n" );
2891 }
2892 else if ( _links.length() == 10 )
2893 {
2894 fprintf( gnuplot, "plot \"datc\" with line, \"dat00\" with line, \"dat01\" with line, "
2895 "\"dat02\" with line, "
2896 "\"dat03\" with line, \"dat04\" with line, \"dat05\" with line, "
2897 "\"dat06\" with line, \"dat07\" "
2898 "with line, \"dat08\" with line, \"dat09\" with line \n" );
2899 }
2900 else if ( _links.length() >= 11 )
2901 {
2902 fprintf( gnuplot,
2903 "plot \"datc\" with line, \"dat00\" with line, \"dat01\" with line, \"dat02\" "
2904 "with line, "
2905 "\"dat03\" with line, \"dat04\" with line, \"dat05\" with line, \"dat06\" with "
2906 "line, \"dat07\" "
2907 "with line, \"dat08\" with line, \"dat09\" with line, \"dat10\" with line \n" );
2908 }
2909 fflush( gnuplot );
2910 char tmp[8];
2911 gets( tmp );
2912 pclose( gnuplot );
2913 }
2914# endif // JJJTEST
2915
2916 return err;
2917}
2918
2919int TTrack::approach2D( TMLink& l ) const {
2920
2921 const TMDCWire& w = *l.wire();
2922 double kappa = _helix->a()[2];
2923 double phi0 = _helix->a()[1];
2924 HepPoint3D xc( _helix->center() );
2925 HepPoint3D xw( w.xyPosition() );
2926 HepPoint3D xt( _helix->x() );
2927 HepVector3D v0( xt - xc );
2928 HepVector3D v1( xw - xc );
2929
2930 double vCrs = v0.x() * v1.y() - v0.y() * v1.x();
2931 double vDot = v0.x() * v1.x() + v0.y() * v1.y();
2932 double dPhi = atan2( vCrs, vDot );
2933
2934 xt = _helix->x( dPhi );
2935 xt.setZ( 0. );
2936 l.positionOnTrack( xt );
2937 xw.setZ( 0. );
2938 l.positionOnWire( xw );
2939 l.dPhi( dPhi );
2940 return 0;
2941}
2942
2943int TTrack::dxda2D( const TMLink& link, double dPhi, Vector& dxda, Vector& dyda ) const {
2944
2945 //...Setup...
2946 double kappa = ( _helix->a() )[2];
2947 if ( kappa == 0. )
2948 {
2949 std::cout << "Error(?) : kappa == 0 in dxda2D of TrkReco." << std::endl;
2950 return 1;
2951 }
2952 const TMDCWire& w = *link.wire();
2953 double dRho = ( _helix->a() )[0];
2954 double phi0 = ( _helix->a() )[1];
2955 // double rho = Helix::ConstantAlpha / kappa;
2956 // double rho = 333.564095 / kappa;
2957 const double Bz = -1000 * getPmgnIMF()->getReferField();
2958 double rho = 333.564095 / ( kappa * Bz );
2959
2960 double sinPhi0 = sin( phi0 );
2961 double cosPhi0 = cos( phi0 );
2962 double sinPhi0dPhi = sin( phi0 + dPhi );
2963 double cosPhi0dPhi = cos( phi0 + dPhi );
2964 Vector dphida( 3 );
2965
2966 HepPoint3D d = _helix->center() - w.xyPosition();
2967 d.setZ( 0. );
2968 double dmag2 = d.mag2();
2969
2970 if ( dmag2 == 0. )
2971 {
2972 std::cout << "Error(?) : Distance(center-xyPosition) == 0 in dxda2D of TrkReco."
2973 << std::endl;
2974 return 1;
2975 }
2976
2977 dphida[0] = ( sinPhi0 * d.x() - cosPhi0 * d.y() ) / dmag2;
2978 dphida[1] = ( dRho + rho ) * ( cosPhi0 * d.x() + sinPhi0 * d.y() ) / dmag2 - 1.;
2979 dphida[2] = ( -rho / kappa ) * ( sinPhi0 * d.x() - cosPhi0 * d.y() ) / dmag2;
2980
2981 dxda[0] = cosPhi0 + rho * sinPhi0dPhi * dphida[0];
2982 dxda[1] = -( dRho + rho ) * sinPhi0 + rho * sinPhi0dPhi * ( 1. + dphida[1] );
2983 dxda[2] = -rho / kappa * ( cosPhi0 - cosPhi0dPhi ) + rho * sinPhi0dPhi * dphida[2];
2984
2985 dyda[0] = sinPhi0 - rho * cosPhi0dPhi * dphida[0];
2986 dyda[1] = ( dRho + rho ) * cosPhi0 - rho * cosPhi0dPhi * ( 1. + dphida[1] );
2987 dyda[2] = -rho / kappa * ( sinPhi0 - sinPhi0dPhi ) - rho * cosPhi0dPhi * dphida[2];
2988
2989 return 0;
2990}
2991
2992int TTrack::dxda2D( double dPhi, Vector& dxda, Vector& dyda ) const {
2993
2994 //...Setup...
2995 double kappa = ( _helix->a() )[2];
2996 if ( kappa == 0. )
2997 {
2998 std::cout << "Error(?) : kappa == 0 in dxda2D of TrkReco." << std::endl;
2999 return 1;
3000 }
3001 double dRho = ( _helix->a() )[0];
3002 double phi0 = ( _helix->a() )[1];
3003 // double rho = Helix::ConstantAlpha / kappa;
3004 // double rho = 333.564095 / kappa;
3005 const double Bz = -1000 * getPmgnIMF()->getReferField();
3006 double rho = 333.564095 / ( kappa * Bz );
3007
3008 double sinPhi0 = sin( phi0 );
3009 double cosPhi0 = cos( phi0 );
3010 double sinPhi0dPhi = sin( phi0 + dPhi );
3011 double cosPhi0dPhi = cos( phi0 + dPhi );
3012 Vector dphida( 3 );
3013
3014 HepPoint3D d = _helix->center(); // d = center - (0,0,0)
3015 d.setZ( 0. );
3016 double dmag2 = d.mag2();
3017
3018 if ( dmag2 == 0. )
3019 {
3020 std::cout << "Error(?) : Distance(center-xyPosition) == 0 in dxda2D of TrkReco."
3021 << std::endl;
3022 return 1;
3023 }
3024
3025 dphida[0] = ( sinPhi0 * d.x() - cosPhi0 * d.y() ) / dmag2;
3026 dphida[1] = ( dRho + rho ) * ( cosPhi0 * d.x() + sinPhi0 * d.y() ) / dmag2 - 1.;
3027 dphida[2] = ( -rho / kappa ) * ( sinPhi0 * d.x() - cosPhi0 * d.y() ) / dmag2;
3028
3029 dxda[0] = cosPhi0 + rho * sinPhi0dPhi * dphida[0];
3030 dxda[1] = -( dRho + rho ) * sinPhi0 + rho * sinPhi0dPhi * ( 1. + dphida[1] );
3031 dxda[2] = -rho / kappa * ( cosPhi0 - cosPhi0dPhi ) + rho * sinPhi0dPhi * dphida[2];
3032
3033 dyda[0] = sinPhi0 - rho * cosPhi0dPhi * dphida[0];
3034 dyda[1] = ( dRho + rho ) * cosPhi0 - rho * cosPhi0dPhi * ( 1. + dphida[1] );
3035 dyda[2] = -rho / kappa * ( sinPhi0 - sinPhi0dPhi ) - rho * cosPhi0dPhi * dphida[2];
3036
3037 return 0;
3038}
3039#endif
3040
3041unsigned TTrack::defineType( void ) const {
3042 bool highPt = true;
3043 if ( pt() < 0.2 ) highPt = false; // Bes
3044 bool fromIP = true;
3045 if ( fabs( impact() ) > 8.3 ) fromIP = false;
3046
3047 if ( fromIP && highPt ) return _type = TrackTypeNormal;
3048 else if ( fromIP && ( !highPt ) ) return _type = TrackTypeCurl;
3049
3050 if ( ( fabs( radius() ) * 2. + fabs( impact() ) ) < 81. ) // Bes
3051 return _type = TrackTypeCircle;
3052 return _type = TrackTypeCosmic;
3053}
3054
3055std::string TrackType( unsigned type ) {
3056 switch ( type )
3057 {
3058 case TrackTypeUndefined: return std::string( "undefined" );
3059 case TrackTypeNormal: return std::string( "normal" );
3060 case TrackTypeCurl: return std::string( "curl " );
3061 case TrackTypeCircle: return std::string( "circle" );
3062 case TrackTypeCosmic: return std::string( "cosmic" );
3063 }
3064 return std::string( "unknown " );
3065}
3066
3067#ifdef OPTNK
3068# define t_dot( a, b ) ( a[0] * b[0] + a[1] * b[1] + a[2] * b[2] )
3069# define t_dot2( a, b ) ( a[0] * b[0] + a[1] * b[1] )
3070# define t_print( a, b ) \
3071 std::cout << b << " = " << a[0] << " " << a[1] << " " << a[2] << std::endl;
3072#endif
3073
3074// #define TRKRECO_DEBUG
3075#ifdef TRKRECO_DEBUG
3076// #include "tuple/BelleTupleManager.h"
3077// BelleHistogram * h_nTrial;
3078#endif
3079
3080int TTrack::approach( TMLink& link, bool doSagCorrection ) const {
3081 //...Cal. dPhi to rotate...
3082 const TMDCWire& w = *link.wire();
3083 double wp[3];
3084 w.xyPosition( wp );
3085 double wb[3];
3086 w.backwardPosition( wb );
3087 double v[3];
3088 v[0] = w.direction().x();
3089 v[1] = w.direction().y();
3090 v[2] = w.direction().z();
3091 //...Sag correction...
3092 if ( doSagCorrection )
3093 {
3094 HepVector3D dir = w.direction();
3095 HepPoint3D xw( wp[0], wp[1], wp[2] );
3096 HepPoint3D wireBackwardPosition( wb[0], wb[1], wb[2] );
3097 w.wirePosition( link.positionOnTrack().z(), xw, wireBackwardPosition, dir );
3098 v[0] = dir.x();
3099 v[1] = dir.y();
3100 v[2] = dir.z();
3101 wp[0] = xw.x();
3102 wp[1] = xw.y();
3103 wp[2] = xw.z();
3104 wb[0] = wireBackwardPosition.x();
3105 wb[1] = wireBackwardPosition.y();
3106 wb[2] = wireBackwardPosition.z();
3107 }
3108 //...Cal. dPhi to rotate...
3109 const HepPoint3D& xc = _helix->center();
3110 double xt[3];
3111 _helix->x( 0., xt );
3112 double x0 = -xc.x();
3113 double y0 = -xc.y();
3114 double x1 = wp[0] + x0;
3115 double y1 = wp[1] + y0;
3116 x0 += xt[0];
3117 y0 += xt[1];
3118 double dPhi = atan2( x0 * y1 - y0 * x1, x0 * x1 + y0 * y1 );
3119 //...Setup...
3120 double kappa = _helix->kappa();
3121 double phi0 = _helix->phi0();
3122// yzhang 2012-08-30
3123//...Axial case...
3124// if (w.axial()) {
3125// link.positionOnTrack(_helix->x(dPhi));
3126// HepPoint3D x(wp[0], wp[1], wp[2]);
3127// x.setZ(link.positionOnTrack().z());
3128// link.positionOnWire(x);
3129// link.dPhi(dPhi);
3130// return 0;
3131// }
3132#ifdef TRKRECO_DEBUG
3133 double firstdfdphi = 0.;
3134 static bool first = true;
3135 if ( first )
3136 {
3137 // extern BelleTupleManager * BASF_Histogram;
3138 // BelleTupleManager * m = BASF_Histogram;
3139 // h_nTrial = m->histogram("TTrack::approach nTrial", 100, 0., 100.);
3140 }
3141#endif
3142
3143 //...Stereo case...
3144 // double rho = Helix::ConstantAlpha / kappa;
3145 // double rho = 333.564095 / kappa;
3146 // yzhang 2012-05-03 delete
3147 // double rho = 333.564095/ kappa;
3148 // yzhang add
3149 Gaudi::svcLocator()->service( "MagneticFieldSvc", m_pmgnIMF );
3150 if ( m_pmgnIMF == NULL )
3151 {
3152 std::cout << __FILE__ << " " << __LINE__ << " Unable to open Magnetic field service"
3153 << std::endl;
3154 }
3155 const double Bz = -1000 * getPmgnIMF()->getReferField();
3156 double rho = 333.564095 / ( kappa * Bz );
3157 // zhangy
3158 double tanLambda = _helix->tanl();
3159 static Vector x( 3 );
3160 double t_x[3];
3161 double t_dXdPhi[3];
3162 const double convergence = 1.0e-5;
3163 double l;
3164 unsigned nTrial = 0;
3165 while ( nTrial < 100 )
3166 {
3167
3168 x = link.positionOnTrack( _helix->x( dPhi ) );
3169 t_x[0] = x[0];
3170 t_x[1] = x[1];
3171 t_x[2] = x[2];
3172
3173 l = v[0] * t_x[0] + v[1] * t_x[1] + v[2] * t_x[2] - v[0] * wb[0] - v[1] * wb[1] -
3174 v[2] * wb[2];
3175
3176 double rcosPhi = rho * cos( phi0 + dPhi );
3177 double rsinPhi = rho * sin( phi0 + dPhi );
3178 t_dXdPhi[0] = rsinPhi;
3179 t_dXdPhi[1] = -rcosPhi;
3180 t_dXdPhi[2] = -rho * tanLambda;
3181
3182 //...f = d(Distance) / d phi...
3183 double t_d2Xd2Phi[2];
3184 t_d2Xd2Phi[0] = rcosPhi;
3185 t_d2Xd2Phi[1] = rsinPhi;
3186
3187 //...iw new...
3188 double n[3];
3189 n[0] = t_x[0] - wb[0];
3190 n[1] = t_x[1] - wb[1];
3191 n[2] = t_x[2] - wb[2];
3192
3193 double a[3];
3194 a[0] = n[0] - l * v[0];
3195 a[1] = n[1] - l * v[1];
3196 a[2] = n[2] - l * v[2];
3197 double dfdphi = a[0] * t_dXdPhi[0] + a[1] * t_dXdPhi[1] + a[2] * t_dXdPhi[2];
3198
3199#ifdef TRKRECO_DEBUG
3200 if ( nTrial == 0 )
3201 {
3202 // break;
3203 firstdfdphi = dfdphi;
3204 }
3205
3206 //...Check bad case...
3207 if ( nTrial > 3 )
3208 {
3209 std::cout << "TTrack::approach ... " << w.name() << " "
3210 << "dfdphi(0)=" << firstdfdphi << ",(" << nTrial << ")=" << dfdphi << endl;
3211 }
3212#endif
3213
3214 //...Is it converged?...
3215 if ( fabs( dfdphi ) < convergence ) break;
3216
3217 double dv = v[0] * t_dXdPhi[0] + v[1] * t_dXdPhi[1] + v[2] * t_dXdPhi[2];
3218 double t0 =
3219 t_dXdPhi[0] * t_dXdPhi[0] + t_dXdPhi[1] * t_dXdPhi[1] + t_dXdPhi[2] * t_dXdPhi[2];
3220 double d2fd2phi = t0 - dv * dv + a[0] * t_d2Xd2Phi[0] + a[1] * t_d2Xd2Phi[1];
3221 // + a[2] * t_d2Xd2Phi[2];
3222
3223 dPhi -= dfdphi / d2fd2phi;
3224
3225 // cout << "nTrial=" << nTrial << endl;
3226 // cout << "iw f,df,dphi=" << dfdphi << "," << d2fd2phi << "," << dPhi << endl;
3227
3228 ++nTrial;
3229 }
3230
3231 //...Cal. positions...
3232 link.positionOnWire( HepPoint3D( wb[0] + l * v[0], wb[1] + l * v[1], wb[2] + l * v[2] ) );
3233 link.dPhi( dPhi );
3234
3235#ifdef TRKRECO_DEBUG
3236// h_nTrial->accumulate((float) nTrial + .5);
3237#endif
3238
3239 return nTrial;
3240}
3241
3242/*
3243void
3244TTrack::relationClusterWithWire(){
3245//----------------------------------------------------------
3246// Note: Selection of associating cluster to TMLink is done.
3247// if appropriate cluster was found ,
3248// relation to cluster is set to the TMLink.
3249//----------------------------------------------------------
3250
3251 int j = 0;
3252 while( TMLink *l = _links[j] ) {
3253
3254 // initialization
3255 int flag(-1);
3256 l->setusecathode(0);
3257 double dPhi = l->dPhi();
3258 l->setZphiBeforeCathode( _helix->x(dPhi).z() );
3259
3260 int k1 = 0;
3261 while( TMDCCatHit *c = _catHits[k1] ){
3262
3263 // Matching of layer
3264 if( c->layerID() != l->hit()->wire()->layerId() ) {
3265 k1++; continue;
3266 }
3267
3268 // Matching of wire hit
3269 AList<TMDCWireHit>& cwire = c->candwire();
3270 AList<TMDCClust>& cclust = c->candclust();
3271
3272 if( cwire.length() == 0 || cclust.length() == 0 ){
3273 k1++; continue;
3274 }
3275
3276 int k2 = 0;
3277 while( TMDCWireHit *cw = cwire[k2] ){
3278 if( cw == l->hit() ){ flag = 1; break; }
3279 k2++;
3280 }
3281
3282 if( flag == -1 ){
3283 k1++; continue;
3284 }
3285
3286 // Selection of cluster
3287 flag = -1;
3288 float cand_sector = CathodeSectorId( cwire[k2]->wire()->id());
3289
3290 // --- debug ---
3291 // std::cout << "cand_sector = " << cand_sector << std::endl;
3292 //--- debug ---
3293
3294 cand_sector = cand_sector - (int(cand_sector/8))*8;
3295
3296 // --- debug ---
3297 // std::cout << "cand_sector(local) = " << cand_sector << std::endl;
3298 // --- debug ---
3299
3300 int count_at_boundary(0);
3301 float tmaxpad_at_boundary(0.);
3302 TMDCClust * cc_at_boundary = NULL;
3303
3304 k2 = 0;
3305 while( TMDCClust *cc = cclust[k2] ){
3306
3307 unsigned cathit_sector = cc->maxpad()->sectorId();
3308
3309 if( cand_sector == float(cathit_sector) ) {
3310 l->setusecathode(1); l->setmclust(cc); break ;
3311 } else if( abs( cand_sector - float(cathit_sector)) == 0.5 ) {
3312
3313 count_at_boundary += 1;
3314 l->setusecathode(1);
3315 l->setmclust(cc);
3316
3317 if( count_at_boundary == 1 ) {
3318 cc_at_boundary = cc;
3319 tmaxpad_at_boundary = cc->tmaxpad();
3320 }
3321
3322 //.. if 2 candidates exist at sector boundary,
3323 // choose the best matching of T between cluster and wire
3324 if( count_at_boundary == 2 ) {
3325 if( fabs( cc->tmaxpad() - l->hit()->reccdc()->m_tdc )
3326 > fabs( tmaxpad_at_boundary - l->hit()->reccdc()->m_tdc ) )
3327 l->setmclust(cc_at_boundary);
3328 break;
3329 }
3330
3331 }
3332
3333 k2++;
3334 }
3335 k1++;
3336 } // TMDCCatHit loop
3337
3338 j++;
3339 } // TMLink loop
3340
3341}
3342
3343// addition by matsu ( 1999/07/05 )
3344void TTrack::relationClusterWithLayer( int SysCorr ){
3345
3346 //... selection of cathode cluster
3347 for(unsigned it0=0;it0<3;it0++){
3348
3349 AList<TMLink> catLink = SameLayer( cores(), it0 );
3350
3351 unsigned n(catLink.length());
3352 if( !n ) continue;
3353
3354 int BestID(-1);
3355
3356 double tmpz(DBL_MAX);
3357 for(unsigned it1=0;it1<n;it1++){
3358 TMLink & l = * catLink[it1];
3359
3360 if( l.getmclust() && l.usecathode() != 0 ){
3361
3362 l.setusecathode(2);
3363
3364 double Zdiff(l.positionOnTrack().z());
3365 if(SysCorr ) Zdiff -= l.getmclust()->zclustWithSysCorr();
3366 else Zdiff -= l.getmclust()->zclust();
3367 if( fabs(Zdiff) < tmpz ){
3368 tmpz = fabs(Zdiff); BestID = it1;
3369 }
3370 }
3371 }
3372
3373 if( BestID == -1 ) continue;
3374
3375 for(unsigned it1=0;it1<n;it1++){
3376 TMLink & l = * catLink[it1];
3377 if( l.getmclust() && l.usecathode() != 0 ){
3378 if( it1 == BestID ) {
3379 l.setusecathode(3); break;
3380 }
3381 }
3382 }
3383 }
3384}
3385*/
3386
3387int TTrack::szPosition( TMLink& link ) const {
3388 const TMDCWireHit& h = *link.hit();
3389 HepVector3D X = 0.5 * ( h.wire()->forwardPosition() + h.wire()->backwardPosition() );
3390 // double theta = atan2(X.y(), X.x());
3391 // HepVector3D lr(h.distance(WireHitLeft) * sin(theta),
3392 // - h.distance(WireHitLeft) * cos(theta),
3393 // 0.);
3394
3395 HepVector3D xx = HepVector3D( X.x(), X.y(), 0. );
3396 HepPoint3D center = _helix->center();
3397 HepVector3D yy = center - xx;
3398 HepVector3D ww = HepVector3D( yy.x(), yy.y(), 0. );
3399 double wwmag2 = ww.mag2();
3400 double wwmag = sqrt( wwmag2 );
3401 HepVector3D lr( h.drift( WireHitLeft ) / wwmag * ww.x(),
3402 h.drift( WireHitLeft ) / wwmag * ww.y(), 0. );
3403
3404#ifdef TRKRECO_DEBUG_DETAIL
3405 std::cout << "old lr " << lr << endl;
3406 std::cout << "old X " << X << endl;
3407 std::cout << "link.leftRight " << link.leftRight() << endl;
3408#endif
3409 //...Check left or right...
3410 // // change to bes3 .. test..
3411 // if (link.leftRight() == WireHitRight) {
3412 // lr = - lr;
3413 //}
3414 // if (link.leftRight() == WireHitLeft) lr = - lr;
3415 // else if (link.leftRight() == 2) lr = ORIGIN;
3416
3417 // yzhang 2012-05-07
3418 const double Bz = -1000 * getPmgnIMF()->getReferField();
3419#ifdef TRKRECO_DEBUG_DETAIL
3420 std::cout << "charge " << _charge << " Bz " << Bz << endl;
3421#endif
3422 if ( _charge * Bz > 0 )
3423 { // yzhang 2012-05-07
3424 if ( link.leftRight() == WireHitRight )
3425 {
3426 lr = -lr; // right
3427 }
3428 }
3429 else
3430 {
3431 if ( link.leftRight() == WireHitLeft )
3432 {
3433 lr = -lr; // left
3434 }
3435 }
3436 if ( link.leftRight() == 2 ) lr = ORIGIN;
3437 // zhangy
3438
3439 X += lr;
3440
3441 //...Prepare vectors...
3442 // HepPoint3D center = _helix->center();
3443 HepPoint3D tmp( -9999., -9999., 0. );
3444 HepVector3D x = HepVector3D( X.x(), X.y(), 0. );
3445 HepVector3D w = x - center;
3446 // //modified the next sentence because the direction are different from belle.
3447 HepVector3D V = h.wire()->direction();
3448 // // to bes3
3449 // // HepVector3D V = - h.wire()->direction();
3450
3451 HepVector3D v = HepVector3D( V.x(), V.y(), 0. );
3452 double vmag2 = v.mag2();
3453 double vmag = sqrt( vmag2 );
3454
3455 double r = _helix->curv();
3456 double wv = w.dot( v );
3457 // //zsl for bes3
3458 // wv = abs(wv);
3459 double d2 = wv * wv - vmag2 * ( w.mag2() - r * r );
3460#ifdef TRKRECO_DEBUG_DETAIL
3461 std::cout << "lr " << lr << endl;
3462 std::cout << "forwardPosition " << h.wire()->forwardPosition() << endl;
3463 std::cout << "backwardPosition " << h.wire()->backwardPosition() << endl;
3464 std::cout << "X " << X << endl;
3465 std::cout << "center " << center << endl;
3466 std::cout << "xx " << xx << endl;
3467 std::cout << "ww " << ww << endl;
3468 std::cout << "TTrack::wire direction:" << h.wire()->direction() << endl;
3469 std::cout << "x " << x << endl;
3470 std::cout << "w " << w << endl;
3471 std::cout << "sz,Track::vmag:" << vmag << ", helix_r:" << r << ", wv:" << wv
3472 << ", d:" << sqrt( d2 ) << endl;
3473#endif
3474
3475 //...No crossing in R/Phi plane... This is too tight...
3476
3477 if ( d2 < 0. )
3478 {
3479 link.position( tmp );
3480
3481#ifdef TRKRECO_DEBUG
3482 std::cout << "TTrack !!! stereo: 0. > d2 = " << d2 << " " << link.leftRight() << std::endl;
3483#endif
3484 return -1;
3485 }
3486 double d = sqrt( d2 );
3487
3488 //...Cal. length to crossing points...
3489 double l[2];
3490 l[0] = ( -wv + d ) / vmag2;
3491 l[1] = ( -wv - d ) / vmag2;
3492
3493 //...Cal. z of crossing points...
3494 bool ok[2];
3495 ok[0] = true;
3496 ok[1] = true;
3497 double z[2];
3498 z[0] = X.z() + l[0] * V.z();
3499 z[1] = X.z() + l[1] * V.z();
3500#ifdef TRKRECO_DEBUG_DETAIL
3501 std::cout << "X.z():" << X.z() << endl;
3502 std::cout << "szPosition::z(0) " << z[0] << " z(1)" << z[1] << " leftRight "
3503 << link.leftRight() << endl;
3504 std::cout << "szPosition::wire backwardPosition and forwardPosition:"
3505 << h.wire()->backwardPosition().z() << "," << h.wire()->forwardPosition().z()
3506 << endl;
3507 std::cout << " l0, l1 = " << l[0] << ", " << l[1] << std::endl;
3508 std::cout << " z0, z1 = " << z[0] << ", " << z[1] << std::endl;
3509 std::cout << " backward = " << h.wire()->backwardPosition().z() << std::endl;
3510 std::cout << " forward = " << h.wire()->forwardPosition().z() << std::endl;
3511#endif
3512
3513 //...Check z position... //yzhang change 2012-05-03
3514 // modified because Belle backward and forward are different from BESIII
3515
3516 /*
3517 if (link.leftRight() == 2) {
3518 if (z[0] > h.wire()->backwardPosition().z()+20.
3519 || z[0] < h.wire()->forwardPosition().z()-20.) ok[0] = false;
3520 if (z[1] > h.wire()->backwardPosition().z()+20.
3521 || z[1] < h.wire()->forwardPosition().z()-20.) ok[1] = false;
3522 }
3523 else {
3524 if (z[0] > h.wire()->backwardPosition().z()
3525 || z[0] < h.wire()->forwardPosition().z() ) ok[0] = false;
3526 if (z[1] > h.wire()->backwardPosition().z()
3527 || z[1] < h.wire()->forwardPosition().z() ) ok[1] = false;
3528 }
3529 if ((! ok[0]) && (! ok[1])) {
3530 link.position(tmp);
3531 return -2;
3532 }*/
3533 // belle...
3534 if ( link.leftRight() == 2 )
3535 {
3536 if ( z[0] < h.wire()->backwardPosition().z() - 20. ||
3537 z[0] > h.wire()->forwardPosition().z() + 20. )
3538 ok[0] = false;
3539 if ( z[1] < h.wire()->backwardPosition().z() - 20. ||
3540 z[1] > h.wire()->forwardPosition().z() + 20. )
3541 ok[1] = false;
3542 }
3543 else
3544 {
3545 if ( z[0] < h.wire()->backwardPosition().z() || z[0] > h.wire()->forwardPosition().z() )
3546 ok[0] = false;
3547 if ( z[1] < h.wire()->backwardPosition().z() || z[1] > h.wire()->forwardPosition().z() )
3548 ok[1] = false;
3549 }
3550 if ( ( !ok[0] ) && ( !ok[1] ) )
3551 {
3552 link.position( tmp );
3553 return -2;
3554 }
3555
3556 //...Cal. xy position of crossing points...
3557 HepVector3D p[2];
3558 p[0] = x + l[0] * v;
3559 p[1] = x + l[1] * v;
3560 /* if (_charge * (center.x() * p[0].y() - center.y() * p[0].x()) < 0.) //liuqg,
3561 cosmic... ok[0] = false; if (_charge * (center.x() * p[1].y() - center.y() * p[1].x()) <
3562 0.) ok[1] = false; if ((! ok[0]) && (! ok[1])){
3563 // double tmp1 = _charge * (center.x() * p[0].y() - center.y() * p[0].x());
3564 // double tmp2 = _charge * (center.x() * p[1].y() - center.y() * p[1].x()) ;
3565 // if (link.leftRight() == 2) std::cout<<tmp1<<" "<<tmp2<<std::endl;
3566 link.position(tmp);
3567 return -3;
3568 }
3569 */
3570 //...Which one is the best?... Study needed...
3571 unsigned best = 0;
3572 if ( ok[1] ) best = 1;
3573
3574 //...Cal. arc length...
3575 double cosdPhi = -center.dot( ( p[best] - center ).unit() ) / center.mag();
3576 double dPhi;
3577 if ( fabs( cosdPhi ) <= 1.0 ) { dPhi = acos( cosdPhi ); }
3578 else if ( cosdPhi > 1.0 ) { dPhi = 0.0; }
3579 else { dPhi = M_PI; }
3580
3581 //...Finish...
3582 tmp.setX( r * dPhi );
3583 tmp.setY( z[best] );
3584 link.position( tmp );
3585
3586 return 0;
3587}
3588
3589int TTrack::szPosition( const TSegment& segment, TMLink& a ) const {
3590
3591 //...Pick up a wire which represents segment position...
3592 AList<TMLink> links = segment.cores();
3593 unsigned n = links.length();
3594 // std::cout<<" szPosition TTrack:: segment.cores().length():"<<n<<endl;
3595 // for (unsigned i = 1; i < n; i++) {
3596 // std::cout<<"drift distance"<<links[i]->hit()->drift()<<endl;
3597 // }
3598 if ( !n ) return -1;
3599 TMLink* minL = links[0];
3600 float minDist = links[0]->drift();
3601#ifdef TRKRECO_DEBUG_DETAIL
3602 std::cout << "minDist szPosition TTrack:" << minDist << endl;
3603#endif
3604 for ( unsigned i = 1; i < n; i++ )
3605 {
3606 if ( links[i]->hit()->drift() < minDist )
3607 {
3608 minDist = links[i]->drift();
3609#ifdef TRKRECO_DEBUG_DETAIL
3610 std::cout << "minDist szPosition TTrack:" << minDist << endl;
3611#endif
3612 minL = links[i];
3613 }
3614 }
3615
3616 //...sz calculation...
3617 a.position( minL->position() );
3618 a.leftRight( 2 );
3619 a.hit( minL->hit() );
3620 int err = szPosition( a );
3621#ifdef TRKRECO_DEBUG_DETAIL
3622 std::cout << "err of szPosition TTrack:" << err << endl;
3623#endif
3624 if ( err ) return -2;
3625 return 0;
3626}
3627
3628int TTrack::szPosition( const HepPoint3D& p, HepPoint3D& sz ) const {
3629
3630 //...Cal. arc length...
3631 HepPoint3D center = _helix->center();
3632 HepPoint3D xy = p;
3633 xy.setZ( 0. );
3634 double cosdPhi = -center.dot( ( xy - center ).unit() ) / center.mag();
3635 double dPhi;
3636 if ( fabs( cosdPhi ) <= 1.0 ) { dPhi = acos( cosdPhi ); }
3637 else if ( cosdPhi > 1.0 ) { dPhi = 0.0; }
3638 else { dPhi = M_PI; }
3639
3640 //...Finish...
3641 sz.setX( _helix->curv() * dPhi );
3642 sz.setY( p.z() );
3643 sz.setZ( 0. );
3644
3645 return 0;
3646}
3647
3648void TTrack::assign( unsigned hitMask ) {
3649 hitMask |= WireHitUsed;
3650
3651 unsigned n = _links.length();
3652 for ( unsigned i = 0; i < n; i++ )
3653 {
3654 TMLink* l = _links[i];
3655 const TMDCWireHit* h = l->hit();
3656#ifdef TRKRECO_DEBUG
3657 if ( h->track() )
3658 {
3659 std::cout << "TTrack::assign !!! hit(" << h->wire()->name();
3660 std::cout << ") already assigned" << std::endl;
3661 }
3662#endif
3663
3664 //...This function will be moved to TTrackManager...
3665 h->track( (TTrack*)this );
3666 h->state( h->state() | hitMask );
3667 }
3668}
3669
3670IBesMagFieldSvc* TTrack::getPmgnIMF() const {
3671 if ( !m_pmgnIMF )
3672 {
3673 StatusCode scmgn = Gaudi::svcLocator()->service( "MagneticFieldSvc", m_pmgnIMF );
3674 if ( scmgn != StatusCode::SUCCESS )
3675 { std::cout << "Unable to open Magnetic field service" << std::endl; }
3676 }
3677 return m_pmgnIMF;
3678}
3679
3680Helix Track2Helix( const MdcTrk_localz& t ) {
3681 HepVector a( 5 );
3682 Hep3Vector p( t.pivot[0], t.pivot[1], t.pivot[2] );
3683 HepSymMatrix er( 5, 0 );
3684 a( 1 ) = t.helix[0];
3685 a( 2 ) = t.helix[1];
3686 a( 3 ) = t.helix[2];
3687 a( 4 ) = t.helix[3];
3688 a( 5 ) = t.helix[4];
3689 er( 1, 1 ) = t.error[0];
3690 er( 2, 1 ) = t.error[1];
3691 er( 2, 2 ) = t.error[2];
3692 er( 3, 1 ) = t.error[3];
3693 er( 3, 2 ) = t.error[4];
3694 er( 3, 3 ) = t.error[5];
3695 er( 4, 1 ) = t.error[6];
3696 er( 4, 2 ) = t.error[7];
3697 er( 4, 3 ) = t.error[8];
3698 er( 4, 4 ) = t.error[9];
3699 er( 5, 1 ) = t.error[10];
3700 er( 5, 2 ) = t.error[11];
3701 er( 5, 3 ) = t.error[12];
3702 er( 5, 4 ) = t.error[13];
3703 er( 5, 5 ) = t.error[14];
3704 return Helix( p, a, er );
3705}
3706
3707Helix Track2Helix( const MdcRec_trk& t ) {
3708 HepVector a( 5 );
3709 Hep3Vector p( t.pivot[0], t.pivot[1], t.pivot[2] );
3710 HepSymMatrix er( 5, 0 );
3711 a( 1 ) = t.helix[0];
3712 a( 2 ) = t.helix[1];
3713 a( 3 ) = t.helix[2];
3714 a( 4 ) = t.helix[3];
3715 a( 5 ) = t.helix[4];
3716 er( 1, 1 ) = t.error[0];
3717 er( 2, 1 ) = t.error[1];
3718 er( 2, 2 ) = t.error[2];
3719 er( 3, 1 ) = t.error[3];
3720 er( 3, 2 ) = t.error[4];
3721 er( 3, 3 ) = t.error[5];
3722 er( 4, 1 ) = t.error[6];
3723 er( 4, 2 ) = t.error[7];
3724 er( 4, 3 ) = t.error[8];
3725 er( 4, 4 ) = t.error[9];
3726 er( 5, 1 ) = t.error[10];
3727 er( 5, 2 ) = t.error[11];
3728 er( 5, 3 ) = t.error[12];
3729 er( 5, 4 ) = t.error[13];
3730 er( 5, 5 ) = t.error[14];
3731 return Helix( p, a, er );
3732}
3733
3734Helix Track2Helix( const Mdst_trk_fit& t ) {
3735 HepVector a( 5 );
3736 Hep3Vector p( t.pivot_x, t.pivot_y, t.pivot_z );
3737 HepSymMatrix er( 5, 0 );
3738 a( 1 ) = t.helix[0];
3739 a( 2 ) = t.helix[1];
3740 a( 3 ) = t.helix[2];
3741 a( 4 ) = t.helix[3];
3742 a( 5 ) = t.helix[4];
3743 er( 1, 1 ) = t.error[0];
3744 er( 2, 1 ) = t.error[1];
3745 er( 2, 2 ) = t.error[2];
3746 er( 3, 1 ) = t.error[3];
3747 er( 3, 2 ) = t.error[4];
3748 er( 3, 3 ) = t.error[5];
3749 er( 4, 1 ) = t.error[6];
3750 er( 4, 2 ) = t.error[7];
3751 er( 4, 3 ) = t.error[8];
3752 er( 4, 4 ) = t.error[9];
3753 er( 5, 1 ) = t.error[10];
3754 er( 5, 2 ) = t.error[11];
3755 er( 5, 3 ) = t.error[12];
3756 er( 5, 4 ) = t.error[13];
3757 er( 5, 5 ) = t.error[14];
3758 return Helix( p, a, er );
3759}
3760
3761std::string TrackKinematics( const Helix& h ) {
3762 static const HepPoint3D IP( 0., 0., 0. );
3763 Helix hIp = h;
3764 hIp.pivot( IP );
3765
3766 float chrg = hIp.a()[2] / fabs( hIp.a()[2] );
3767 std::string s;
3768 if ( chrg > 0. ) s = "+";
3769 else s = "-";
3770
3771 float x[4];
3772 x[0] = fabs( hIp.a()[0] );
3773 x[1] = hIp.a()[3];
3774 x[2] = 1. / fabs( hIp.a()[2] );
3775 x[3] = ( 1. / fabs( hIp.a()[2] ) ) * hIp.a()[4];
3776
3777 if ( ( x[0] < 2. ) && ( fabs( x[1] ) < 4. ) ) s += "i ";
3778 else s += " ";
3779
3780 for ( unsigned i = 0; i < 4; i++ )
3781 {
3782 if ( i ) s += " ";
3783
3784 if ( x[i] < 0. ) s += "-";
3785 else s += " ";
3786
3787 int y = int( fabs( x[i] ) );
3788 s += std::to_string( y ) + ".";
3789 float z = fabs( x[i] );
3790 for ( unsigned j = 0; j < 3; j++ )
3791 {
3792 z *= 10.;
3793 y = ( int( z ) % 10 );
3794 s += std::to_string( y );
3795 }
3796 }
3797
3798 return s;
3799}
3800
3801std::string TrackStatus( const TTrack& t ) {
3802 return TrackStatus( t.finder(), t.type(), t.quality(), t.fitting(), 0, 0 );
3803}
3804
3805std::string TrackStatus( const MdcRec_trk& c ) {
3806 // const reccdc_trk_add & a =
3807 // * (reccdc_trk_add *) BsGetEnt(RECMDC_TRK_ADD, c.m_ID, BBS_No_Index);
3808 // return TrackStatus(a);
3809 return std::string( "" );
3810}
3811
3812std::string TrackStatus( const MdcRec_trk_add& a ) {
3813 // return TrackStatus(a.decision,
3814 // a.kind,
3815 // a.quality,
3816 // a.stat,
3817 // a.mother,
3818 // a.daughter);
3819 return std::string( "" );
3820}
3821
3822std::string TrackStatus( unsigned md, unsigned mk, unsigned mq, unsigned ms, unsigned mm,
3823 unsigned ma ) {
3824
3825 std::string f;
3826 if ( md & TrackOldConformalFinder ) f += "o";
3827 if ( md & TrackFastFinder ) f += "f";
3828 if ( md & TrackSlowFinder ) f += "s";
3829 if ( md & TrackCurlFinder ) f += "c";
3830 if ( md & TrackTrackManager ) f += "t";
3831 if ( f == "" ) f = "?";
3832
3833 std::string k;
3834 if ( mk & TrackTypeNormal ) k += "Norm";
3835 if ( mk & TrackTypeCurl ) k += "Curl";
3836 if ( mk & TrackTypeCircle ) k += "Circ";
3837 if ( mk & TrackTypeIncomingCosmic ) k += "Inco";
3838 if ( mk & TrackTypeOutgoingCosmic ) k += "Outc";
3839 if ( mk & TrackTypeKink ) k += "Kink";
3840 if ( mk & TrackTypeSVDOnly ) k += "Svd";
3841 if ( k == "" ) k = "?";
3842
3843 std::string b;
3844 if ( mq & TrackQualityOutsideCurler ) b += "Curlback";
3845 if ( mq & TrackQualityAfterKink ) b += "Afterkink";
3846 if ( mq & TrackQualityCosmic ) b += "Cosmic";
3847 if ( mq & TrackQuality2D ) b += "2D";
3848 if ( b == "" ) b = "ok";
3849
3850 std::string s;
3851 if ( ms & TrackFitGlobal ) s += "HFit";
3852 if ( ms & TrackFitCosmic ) s += "CFit";
3853 if ( ms & TrackFitCdcKalman ) s += "CKal";
3854 if ( ms & TrackFitSvdCdcKalman ) s += "SKal";
3855 if ( s == "" ) s = "?";
3856
3857 int m = mm;
3858 if ( m ) --m;
3859
3860 int d = ma;
3861 if ( d ) --d;
3862
3863 std::string p = " ";
3864 return f + p + k + p + b + p + s + p + std::to_string( m ) + p + std::to_string( d );
3865}
3866
3867std::string TrackInformation( const TTrack& t ) {
3868 const AList<TMLink> cores = t.cores();
3869 unsigned n = cores.length();
3870 unsigned nS = NStereoHits( cores );
3871 unsigned nA = n - nS;
3872 std::string p;
3873 if ( !PositiveDefinite( t.helix() ) ) p = " negative";
3874 if ( HelixHasNan( t.helix() ) ) p += " NaN";
3875 return TrackInformation( nA, nS, n, t.chi2() ) + p;
3876}
3877
3878std::string TrackInformation( const MdcRec_trk& r ) {
3879 std::string p;
3880 if ( PositiveDefinite( Track2Helix( r ) ) ) p = " posi";
3881 else p = " nega";
3882 if ( HelixHasNan( Track2Helix( r ) ) ) p += " with NaN";
3883 return TrackInformation( r.nhits - r.nster, r.nster, r.nhits, r.chiSq ) + p;
3884}
3885
3886std::string TrackInformation( unsigned nA, unsigned nS, unsigned n, float chisq ) {
3887 std::string s;
3888
3889 s += "a" + std::to_string( int( nA ) );
3890 s += " s" + std::to_string( int( nS ) );
3891 s += " n" + std::to_string( int( n ) );
3892 // s += " ndf" + std::string(int(r.m_ndf));
3893 float x = chisq;
3894
3895 if ( x < 0. ) s += " -";
3896 else s += " ";
3897
3898 int y = int( fabs( x ) );
3899 s += std::to_string( y ) + ".";
3900 float z = fabs( x );
3901 for ( unsigned j = 0; j < 1; j++ )
3902 {
3903 z *= 10.;
3904 y = ( int( z ) % 10 );
3905 s += std::to_string( y );
3906 }
3907
3908 return s;
3909}
3910
3911std::string TrackLayerUsage( const TTrack& t ) {
3912 unsigned n[11];
3913 NHitsSuperLayer( t.links(), n );
3914 std::string nh;
3915 for ( unsigned i = 0; i < 11; i++ )
3916 {
3917 nh += std::to_string( n[i] );
3918 if ( i % 2 ) nh += "-";
3919 else if ( i < 10 ) nh += ",";
3920 }
3921 return nh;
3922}
3923
3924bool PositiveDefinite( const Helix& h ) {
3925 const SymMatrix& e = h.Ea();
3926 SymMatrix e2 = e.sub( 1, 2 );
3927 SymMatrix e3 = e.sub( 1, 3 );
3928 SymMatrix e4 = e.sub( 1, 4 );
3929
3930 bool positive = true;
3931 if ( e[0][0] <= 0. ) positive = false;
3932 else if ( e2.determinant() <= 0. ) positive = false;
3933 else if ( e3.determinant() <= 0. ) positive = false;
3934 else if ( e4.determinant() <= 0. ) positive = false;
3935 else if ( e.determinant() <= 0. ) positive = false;
3936 return positive;
3937}
3938
3939bool HelixHasNan( const Helix& h ) {
3940 const Vector& a = h.a();
3941 for ( unsigned i = 0; i < 5; i++ )
3942 // maqm if (isnan(a[i]))
3943 if ( std::isnan( a[i] ) ) return true;
3944 const SymMatrix& Ea = h.Ea();
3945 for ( unsigned i = 0; i < 5; i++ )
3946 for ( unsigned j = 0; j <= i; j++ )
3947 // maqm if (isnan(Ea[i][j]))
3948 if ( std::isnan( Ea[i][j] ) ) return true;
3949 return false;
3950}
3951
3952Helix Track2Helix( const Gen_hepevt& t ) {
3953 float charge = 1;
3954 if ( t.idhep == 11 ) charge = -1;
3955 else if ( t.idhep == -11 ) charge = 1;
3956 else if ( t.idhep == 13 ) charge = -1;
3957 else if ( t.idhep == -13 ) charge = 1;
3958 else if ( t.idhep == 211 ) charge = 1;
3959 else if ( t.idhep == -211 ) charge = -1;
3960 else if ( t.idhep == 321 ) charge = 1;
3961 else if ( t.idhep == -321 ) charge = -1;
3962 else if ( t.idhep == 2212 ) charge = 1;
3963 else if ( t.idhep == -2212 ) charge = -1;
3964 else
3965 {
3966 std::cout << "Track2Helix(gen_hepevt) !!! charge of id=";
3967 std::cout << t.idhep << " is unknown" << std::endl;
3968 }
3969
3970 Hep3Vector mom( t.P[0], t.P[1], t.P[2] );
3971 Hep3Vector pos( t.V[0] / 10., t.V[1] / 10., t.V[2] / 10. );
3972 return Helix( pos, mom, charge );
3973}
HepVector3D
HepGeom::Vector3D< double > HepVector3D
Definition Analysis/VertexFit/src/BField.cxx:13
test
std::string test
Definition CalibModel.cxx:44
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
sprintf
sprintf(cut, "kal_costheta0_em>-0.93&&kal_costheta0_em<0.93&&kal_pxy0_em>=0.05+%d*0.1&&kal_" "pxy0_em<0.15+%d*0.1&&NGch>=2", j, j)
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
n
const Int_t n
Definition DataBase/tau_mode.c:57
data
TTree * data
Definition DataBase/tau_mode.c:5
x
Double_t x[10]
Definition DataBase/tau_mode.c:49
e2
Double_t e2
Definition DataBase/tau_mode.c:7
w
double w
Definition EvtPhokharaDef.hh:28
unit
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in !Latex Output unit
Definition FoamA.h:90
mk
const double mk
Definition Gam4pikp.cxx:33
s
XmlRpcServer s
Definition HelloServer.cpp:11
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
zero
const DifNumber zero
Matrix
CLHEP::HepMatrix Matrix
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:43
SymMatrix
CLHEP::HepSymMatrix SymMatrix
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:39
Vector
CLHEP::HepVector Vector
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:38
intersection
int intersection(const HepPoint3D &c1, double r1, const HepPoint3D &c2, double r2, double eps, HepPoint3D &x1, HepPoint3D &x2)
Circle utilities.
Definition TMDCUtil.cxx:93
ORIGIN
const HepPoint3D ORIGIN
Constants.
Definition TMDCUtil.cxx:47
WireHitLeft
#define WireHitLeft
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:21
WireHitRight
#define WireHitRight
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:22
WireHitUsed
#define WireHitUsed
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:47
NHitsSuperLayer
void NHitsSuperLayer(const AList< TMLink > &links, unsigned nHits[11])
returns # of hits per super layer.
Definition TMLink.cxx:163
NStereoHits
unsigned NStereoHits(const AList< TMLink > &links)
returns # of stereo hits.
Definition TMLink.cxx:273
TrackFitCosmic
#define TrackFitCosmic
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:53
TrackFitCdcKalman
#define TrackFitCdcKalman
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:54
TrackFitSvdCdcKalman
#define TrackFitSvdCdcKalman
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:55
TrackQualityAfterKink
#define TrackQualityAfterKink
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:45
TrackTypeCircle
#define TrackTypeCircle
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:36
TrackFastFinder
#define TrackFastFinder
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:24
TrackCurlFinder
#define TrackCurlFinder
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:26
TrackTypeNormal
#define TrackTypeNormal
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:34
TrackQuality2D
#define TrackQuality2D
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:47
TrackQualityCosmic
#define TrackQualityCosmic
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:46
TrackTypeUndefined
#define TrackTypeUndefined
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:33
TrackTypeKink
#define TrackTypeKink
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:40
TrackTrackManager
#define TrackTrackManager
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:27
TrackSlowFinder
#define TrackSlowFinder
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:25
TrackTypeCosmic
#define TrackTypeCosmic
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:37
TrackTypeCurl
#define TrackTypeCurl
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:35
TrackTypeIncomingCosmic
#define TrackTypeIncomingCosmic
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:38
TrackOldConformalFinder
#define TrackOldConformalFinder
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:23
TrackTypeSVDOnly
#define TrackTypeSVDOnly
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:41
TrackQualityOutsideCurler
#define TrackQualityOutsideCurler
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:44
TrackTypeOutgoingCosmic
#define TrackTypeOutgoingCosmic
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:39
TrackFitGlobal
#define TrackFitGlobal
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:52
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
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition Reconstruction/TrkReco/include/TrkReco/TTrack.h:102
n2
int n2
Definition SD0Tag.cxx:59
n1
int n1
Definition SD0Tag.cxx:58
M_PI
#define M_PI
Definition TConstant.h:4
Track2Helix
Helix Track2Helix(const MdcTrk_localz &t)
Definition TTrack.cxx:3680
SortByPt
int SortByPt(const void *av, const void *bv)
Utility functions.
Definition TTrack.cxx:2514
TrackInformation
std::string TrackInformation(const TTrack &t)
Definition TTrack.cxx:3867
TrackStatus
std::string TrackStatus(const TTrack &t)
returns string of track status.
Definition TTrack.cxx:3801
TrackLayerUsage
std::string TrackLayerUsage(const TTrack &t)
Definition TTrack.cxx:3911
PositiveDefinite
bool PositiveDefinite(const Helix &h)
Error matrix validity.
Definition TTrack.cxx:3924
TrackKinematics
std::string TrackKinematics(const Helix &h)
Definition TTrack.cxx:3761
HelixHasNan
bool HelixHasNan(const Helix &h)
Helix parameter validity.
Definition TTrack.cxx:3939
Gen_hepevt
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcTables/HepevtTables.h:13
Helix::center
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:197
Helix::Ea
const HepSymMatrix & Ea(void) const
returns error matrix.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:219
Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition Reconstruction/TrackUtil/src/Helix.cxx:120
Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:217
Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:199
IBesMagFieldSvc::getReferField
virtual double getReferField()=0
MdcRec_trk_add
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcTables/MdcTables.h:763
MdcRec_trk
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcTables/MdcTables.h:374
MdcRec_trk::chiSq
float chiSq
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcTables/MdcTables.h:426
MdcRec_trk::nhits
int nhits
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcTables/MdcTables.h:429
MdcRec_trk::nster
int nster
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcTables/MdcTables.h:430
MdcTrk_localz
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcTables/TrkTables.h:176
Mdst_trk_fit
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcTables/MdstTables.h:373
TCircle
A class to represent a circle in tracking.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TCircle.h:42
TCircle::center
const HepPoint3D & center(void) const
returns position of center.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TCircle.h:106
THelixFitter
A class to fit a TTrackBase object to a helix.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/THelixFitter.h:47
TMDCWireHit::dump
void dump(const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
dumps debug information.
Definition TMDCWireHit.cxx:62
TMDCWireHit::track
const TTrack *const track(void) const
assigns a pointer to a TTrack.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:240
TMDCWireHit::dDrift
float dDrift(unsigned) const
returns drift distance error.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:223
TMDCWireHit::drift
float drift(unsigned) const
returns drift distance.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:218
TMDCWireHit::wire
const TMDCWire *const wire(void) const
returns a pointer to a TMDCWire.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:212
TMDCWire::layerId
unsigned layerId(void) const
returns layer id.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:207
TMDCWire::direction
const HepVector3D & direction(void) const
returns direction vector of the wire.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:264
TMDCWire::xyPosition
const HepPoint3D & xyPosition(void) const
returns middle position of a wire. z componet is 0.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:255
TMDCWire::forwardPosition
const HepPoint3D & forwardPosition(void) const
returns position in forward endplate.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:242
TMDCWire::backwardPosition
const HepPoint3D & backwardPosition(void) const
returns position in backward endplate.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:244
TMDCWire::name
std::string name(void) const
returns name.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:308
TRunge
A class to represent a track in tracking.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TRunge.h:65
TSegment
A class to relate TMDCWireHit and TTrack objects.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TSegment.h:43
TTrackBase::distance
virtual double distance(const TMLink &) const
returns distance to a position of TMLink in TMLink space.
Definition TTrackBase.cxx:86
TTrackBase::dump
virtual void dump(const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
dumps debug information.
Definition TTrackBase.cxx:58
TTrackBase::links
const AList< TMLink > & links(unsigned mask=0) const
Definition TTrackBase.cxx:289
TTrackBase::fitter
const TMFitter *const fitter(void) const
returns a pointer to a default fitter.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrackBase.h:236
TTrackBase::nCores
unsigned nCores(unsigned mask=0) const
returns # of masked TMLinks for fit. 'mask' will be applied if mask is not 0.
Definition TTrackBase.cxx:314
TTrackBase::TTrackBase
TTrackBase()
Constructor.
Definition TTrackBase.cxx:41
TTrackBase::cores
const AList< TMLink > & cores(unsigned mask=0) const
returns a list of masked TMLinks for fit. 'mask' will be applied if mask is not 0.
Definition TTrackBase.cxx:308
TTrack
A class to represent a track in tracking.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:128
TTrack::assign
void assign(unsigned maskForWireHit)
assigns wire hits to this track.
Definition TTrack.cxx:3648
TTrack::helix
const Helix & helix(void) const
returns helix parameter.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:434
TTrack::dump
void dump(const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
dumps debug information.
Definition TTrack.cxx:204
TTrack::stereoHitForCurl
int stereoHitForCurl(TMLink &link, AList< HepPoint3D > &arcZList) const
TTrack::name
const std::string & name(void) const
returns/sets name.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:459
TTrack::center
TPoint2D center(void) const
returns position of helix center.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:493
TTrack::quality
unsigned quality(void) const
sets/returns quality.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:502
TTrack::TTrack
TTrack()
Default constructor.
Definition TTrack.cxx:190
TTrack::fit2D
int fit2D(unsigned=0, double=0.1, double=0.015)
fits itself. Error was happened if return value is not zero.
Definition TTrack.cxx:2524
TTrack::radius
double radius(void) const
returns signed radius.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:487
TTrack::HelCyl
int HelCyl(double rhole, double rcyl, double zb, double zf, double epsl, double &phi, HepPoint3D &xp) const
returns a cathode hit list.
Definition TTrack.cxx:283
TTrack::impact
double impact(void) const
returns signed impact parameter to the origin.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:483
TTrack::approach
int approach(TMLink &) const
Definition TTrack.cxx:265
TTrack::refine2D
void refine2D(AList< TMLink > &list, float maxSigma)
fits itself with cathode hits.
Definition TTrack.cxx:267
TTrack::approach2D
int approach2D(TMLink &) const
Definition TTrack.cxx:2919
TTrack::movePivot
void movePivot(void)
moves pivot to the inner most hit.
Definition TTrack.cxx:230
TTrack::szPosition
int szPosition(TMLink &link) const
calculates arc length and z for a stereo hit.
Definition TTrack.cxx:3387
TTrack::charge
double charge(void) const
returns charge.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:455
TTrack::~TTrack
virtual ~TTrack()
Destructor.
Definition TTrack.cxx:202
TTrack::deleteListForCurl
void deleteListForCurl(AList< HepPoint3D > &l1, AList< HepPoint3D > &l2) const
TTrack::p
Hep3Vector p(void) const
returns momentum.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:477
rb::R
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition TUtil.h:22
t
int t()
Definition t.c:1