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