BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
RecMucTrack.cxx
Go to the documentation of this file.
1//$id$
2//
3//$log$
4
5/*
6 * 2004/03/08 Zhengyun You Peking University
7 *
8 * 2004/09/12 Zhengyun You Peking University
9 * transplanted to Gaudi framework
10 */
11
12using namespace std;
13
14#include "geomdefs.hh"
15#include <CLHEP/Geometry/Point3D.h>
16#include <CLHEP/Vector/ThreeVector.h>
17#include <iostream>
18#include <vector>
19
20#include "MucGeomSvc/MucConstant.h"
21#include "MucRecEvent/MucRec2DRoad.h"
22#include "MucRecEvent/MucRec3DRoad.h"
23#include "MucRecEvent/MucTrackParameter.h"
24#include "MucRecEvent/RecMucTrack.h"
25
26// Constructor.
27RecMucTrack::RecMucTrack()
28 : // m_trackId(-1), //--------------->
29 m_ExtTrackID( -1 )
30 , m_MdcPos( 0.0, 0.0, 0.0 )
31 , m_MdcMomentum( 0.0, 0.0, 0.0 )
32 , m_MucPos( 0.0, 0.0, 0.0 )
33 , m_MucPosSigma( 0.0, 0.0, 0.0 )
34 , m_MucMomentum( 0.0, 0.0, 0.0 )
35 , m_CurrentPos( 0.0, 0.0, 0.0 )
36 , m_CurrentDir( 0.0, 0.0, 0.0 )
37 , m_CurrentInsct( 0.0, 0.0, 0.0 )
38 , m_Good3DLine( 0 )
39 , m_pHits( 0 )
40 , m_pExpectedHits( 0 )
41 , m_Intersections( 0 )
42 , m_Directions( 0 ) {
43 // initialize m_IntersectionInner/Outer.
44 for ( int igap = 0; igap < 9; igap++ )
45 {
46 m_IntersectionInner[igap].set( -9999, -9999, -9999 );
47 m_IntersectionOuter[igap].set( -9999, -9999, -9999 );
48 }
49 m_id = 0;
50 m_status = -1;
51 m_type = -1;
52
53 m_numHits = 0;
54 m_startPart = -1;
55 m_endPart = -1;
56 m_brLastLayer = -1;
57 m_ecLastLayer = -1;
58 m_brFirstLayer = -1;
59 m_ecFirstLayer = -1;
60 m_ecPart = -1;
61 m_numLayers = 0;
62 m_maxHitsInLayer = 0;
63 m_depth = -99;
64 m_dof = 0;
65 m_chi2 = 0.0;
66 m_rms = 0.0;
67 m_deltaPhi = 0.0;
68
69 m_xPosSigma = 0.0;
70 m_yPosSigma = 0.0;
71 m_zPosSigma = 0.0;
72
73 m_changeUnit = false;
74 m_recmode = 0;
75 // added by LI Chunhua
76 m_kalrechi2 = 0.;
77 m_kaldof = 0;
78 m_kaldepth = -99;
79 m_kalbrLastLayer = -1;
80 m_kalecLastLayer = -1;
81}
82
83// Assignment constructor.
84RecMucTrack& RecMucTrack::operator=( const RecMucTrack& orig ) {
85 // Assignment operator.
86 if ( this != &orig )
87 { // Watch out for self-assignment!
88 // m_trackId = orig.m_trackId; //--------------->
89 m_ExtTrackID = orig.m_ExtTrackID;
90 m_MdcPos = orig.m_MdcPos;
91 m_MdcMomentum = orig.m_MdcMomentum;
92 m_MucPos = orig.m_MucPos;
93 m_MucPosSigma = orig.m_MucPosSigma;
94 m_MucMomentum = orig.m_MucMomentum;
95 m_CurrentPos = orig.m_CurrentPos;
96 m_CurrentDir = orig.m_CurrentDir;
97 m_CurrentInsct = orig.m_CurrentInsct;
98 m_id = orig.m_id;
99 m_status = orig.m_status;
100 m_type = orig.m_type;
101 m_numHits = orig.m_numHits; //--------------->
102 m_startPart = orig.m_startPart;
103 m_endPart = orig.m_endPart;
104 m_brLastLayer = orig.m_brLastLayer;
105 m_ecLastLayer = orig.m_ecLastLayer;
106 m_numLayers = orig.m_numLayers;
107 m_maxHitsInLayer = orig.m_maxHitsInLayer;
108 m_depth = orig.m_depth;
109 m_dof = orig.m_dof;
110 m_chi2 = orig.m_chi2;
111 m_rms = orig.m_rms;
112 m_deltaPhi = orig.m_deltaPhi;
113 m_pHits = orig.m_pHits;
114 m_pExpectedHits = orig.m_pExpectedHits;
115 m_Intersections = orig.m_Intersections;
116 m_Directions = orig.m_Directions;
117 m_xPosSigma = orig.m_xPosSigma;
118 m_yPosSigma = orig.m_yPosSigma;
119 m_zPosSigma = orig.m_zPosSigma;
120 m_changeUnit = orig.m_changeUnit;
121 m_recmode = orig.m_recmode;
122 //*******
123 m_kalrechi2 = orig.m_kalrechi2;
124 m_kaldof = orig.m_kaldof;
125 m_kaldepth = orig.m_kaldepth;
126 m_kalbrLastLayer = orig.m_kalbrLastLayer;
127 m_kalecLastLayer = orig.m_kalecLastLayer;
128 }
129
130 return *this;
131}
132
133// Copy constructor.
134RecMucTrack::RecMucTrack( const RecMucTrack& source )
135 : // m_trackId (source.m_trackId), //--------------->
136 m_ExtTrackID( source.m_ExtTrackID )
137 , m_MdcPos( source.m_MdcPos )
138 , m_MdcMomentum( source.m_MdcMomentum )
139 , m_MucPos( source.m_MucPos )
140 , m_MucPosSigma( source.m_MucPosSigma )
141 , m_MucMomentum( source.m_MucMomentum )
142 , m_CurrentPos( source.m_CurrentPos )
143 , m_CurrentDir( source.m_CurrentDir )
144 , m_CurrentInsct( source.m_CurrentInsct )
145 , m_pHits( source.m_pHits )
146 , m_pExpectedHits( source.m_pExpectedHits )
147 , m_Intersections( source.m_Intersections )
148 , m_Directions( source.m_Directions ) {
149 m_id = source.m_id;
150 m_status = source.m_status;
151 m_type = source.m_type;
152 m_numHits = source.m_numHits; //--------------->
153 m_startPart = source.m_startPart;
154 m_endPart = source.m_endPart;
155 m_brLastLayer = source.m_brLastLayer;
156 m_ecLastLayer = source.m_ecLastLayer;
157 m_numLayers = source.m_numLayers;
158 m_maxHitsInLayer = source.m_maxHitsInLayer;
159 m_depth = source.m_depth;
160 m_dof = source.m_dof;
161 m_chi2 = source.m_chi2;
162 m_rms = source.m_rms;
163 m_deltaPhi = source.m_deltaPhi;
164 m_xPosSigma = source.m_xPosSigma;
165 m_yPosSigma = source.m_yPosSigma;
166 m_zPosSigma = source.m_zPosSigma;
167 m_changeUnit = source.m_changeUnit;
168 m_recmode = source.m_recmode;
169 //******
170 m_kalrechi2 = source.m_kalrechi2;
171 m_kaldof = source.m_kaldof;
172 m_kaldepth = source.m_kaldepth;
173 m_kalbrLastLayer = source.m_kalbrLastLayer;
174 m_kalecLastLayer = source.m_kalecLastLayer;
175}
176
177RecMucTrack::RecMucTrack( const DstMucTrack& dstTrack ) : DstMucTrack( dstTrack ) {
178 SetDefault();
179}
180
181RecMucTrack& RecMucTrack::operator=( const DstMucTrack& dstTrack ) {
182 SetDefault();
183 DstMucTrack::operator=( dstTrack );
184 return *this;
185}
186
187void RecMucTrack::SetDefault() {
188 m_brFirstLayer = -99;
189 m_ecFirstLayer = -99;
190 m_Good3DPart = -99;
191}
192
193// Destructor.
194RecMucTrack::~RecMucTrack() {
195 for ( int i = 0; i < m_pExpectedHits.size(); i++ ) delete m_pExpectedHits[i];
196}
197
198// Set the index for this track.
199void RecMucTrack::setTrackId( const int trackId ) {
200 if ( trackId >= 0 ) { m_trackId = trackId; }
201}
202
203// now, what we get is cm, so it' multiplied by 10 to convert to mm
204// and GeV -> MeV
205void RecMucTrack::SetMdcPos( const float x, const float y, const float z ) {
206 m_MdcPos.set( x * 10, y * 10, z * 10 );
207}
208
209void RecMucTrack::SetMdcMomentum( const float px, const float py, const float pz ) {
210 m_MdcMomentum.set( px * 1000, py * 1000, pz * 1000 );
211}
212
213void RecMucTrack::SetMucPos( const float x, const float y, const float z ) {
214 m_MucPos.set( x * 10, y * 10, z * 10 );
215 m_xPos = x * 10;
216 m_yPos = y * 10;
217 m_zPos = z * 10;
218}
219
220void RecMucTrack::SetMucPosSigma( const float x, const float y, const float z ) {
221 m_MucPosSigma.set( x * 10, y * 10, z * 10 );
222 m_xPosSigma = x * 10;
223 m_yPosSigma = y * 10;
224 m_zPosSigma = z * 10;
225}
226
227void RecMucTrack::SetMucMomentum( const float px, const float py, const float pz ) {
228 m_MucMomentum.set( px * 1000, py * 1000, pz * 1000 );
229 m_px = px * 1000;
230 m_py = py * 1000;
231 m_pz = pz * 1000;
232}
233
234void RecMucTrack::SetExtMucPos( const float x, const float y, const float z ) {
235 m_ExtMucPos.set( x * 10, y * 10, z * 10 );
236}
237
238void RecMucTrack::SetExtMucMomentum( const float px, const float py, const float pz ) {
239 m_ExtMucMomentum.set( px * 1000, py * 1000, pz * 1000 );
240}
241
242void RecMucTrack::SetCurrentPos( const float x, const float y, const float z ) {
243 m_CurrentPos.set( x * 10, y * 10, z * 10 );
244}
245
246void RecMucTrack::SetCurrentDir( const float x, const float y, const float z ) {
247 m_CurrentDir.set( x * 1000, y * 1000, z * 1000 ); // unnecessary, because it's dir, not mom
248}
249
250void RecMucTrack::SetCurrentInsct( const float x, const float y, const float z ) {
251 m_CurrentInsct.set( x, y,
252 z ); /////this is intenal function, it receive "mm", so need not convert.
253}
254
255// set corresponding monte carlo track pointer.
256// RecMucTrack::SetMCTrack(const BesPersTrack* mcTrack);
257//{
258// m_MCTrack = mcTrack;
259//}
260
261void RecMucTrack::SetExtTrack( RecExtTrack* extTrack ) {
262 m_ExtTrack = extTrack;
263 if ( m_ExtTrack ) m_ExtTrackID = extTrack->GetTrackId();
264}
265
266// reload setVecHits
267
268void RecMucTrack::setVecHits( vector<MucRecHit*>& pHits ) {
269 for ( int i = 0; i < pHits.size(); i++ )
270 m_vecHits.push_back( ( pHits[i]->GetID() ).get_value() );
271}
272
273void RecMucTrack::setExpHits( vector<MucRecHit*>& pHits ) {
274 for ( int i = 0; i < pHits.size(); i++ )
275 m_expHits.push_back( ( pHits[i]->GetID() ).get_value() );
276}
277
278void RecMucTrack::GetMdcExtTrack( Hep3Vector mdcStartPos, Hep3Vector mdcStartMomentum,
279 int charge, Hep3Vector& mucStartPos,
280 Hep3Vector& mucStartMomentum ) {
281
282 // cm->mm; GeV->MeV//
283 mdcStartPos *= 10;
284 mdcStartMomentum *= 1000;
285 ////////////////////
286
287 Hep3Vector pt = mdcStartMomentum;
288 pt.setZ( 0 );
289 // cout << "pt " << pt.mag() << endl;
290 double radius = ( pt.mag() * 1e6 / kvC * 1e3 ) / ( fabs( charge * kMagnetField ) );
291 // double startPhi = startP.phi();
292 double deltaPhi = -1.0 * ( charge / abs( charge ) ) * kDeltaPhi;
293 double deltaZ = ( mdcStartMomentum.z() * 1e6 / kvC * 1e3 ) * kDeltaPhi /
294 ( abs( charge ) * kMagnetField );
295
296 // cout << "r = " << radius << endl;
297 mucStartPos = mdcStartPos;
298 mucStartMomentum = mdcStartMomentum;
299 double phi;
300 int iter = 0;
301 do {
302 phi = mucStartMomentum.getPhi() + deltaPhi;
303 mucStartPos.setX( mucStartPos.x() + radius * kDeltaPhi * cos( phi ) );
304 mucStartPos.setY( mucStartPos.y() + radius * kDeltaPhi * sin( phi ) );
305 mucStartPos.setZ( mucStartPos.z() + deltaZ );
306
307 mucStartMomentum.setPhi( mucStartMomentum.phi() + deltaPhi );
308 iter++;
309 // cout << endP << " " << mucStartPos << endl;
310 } while ( IsInsideSuperConductor( mucStartPos ) && iter < kMdcExtIterMax );
311
312 // mm->cm; MeV->GeV//
313 mucStartPos /= 10;
314 mucStartMomentum /= 1000;
315 ////////////////////
316}
317
318bool RecMucTrack::IsInsideSuperConductor( Hep3Vector pos ) {
319 if ( pos.mag() < kSuperConductorR && fabs( pos.z() ) < kSuperConductorZ ) { return true; }
320 else { return false; }
321}
322
323// Attach the given hit to this track.
324// Assume that this hit has been verified to be consistent with the road.
325void RecMucTrack::AttachHit( MucRecHit* hit ) {
326 // cout << "Muc2DRoad::AttachHit-I0 hit = " << hit << endl;
327
328 if ( !hit )
329 {
330 cout << "RecMucTrack::AttachHit-E1 null hit pointer!" << endl;
331 return;
332 }
333
334 int part = hit->Part();
335 int gap = hit->Gap();
336 if ( ( gap < 0 ) || ( gap >= (int)MucID::getGapNum( part ) ) )
337 {
338 // The gap number of the hit is out of range.
339 cout << "Muc2DRoad::AttachHit(MucRecHit*), bad gap number = " << gap << endl;
340 return;
341 }
342
343 // Attach the hit to the road.
344 m_pHits.push_back( hit );
345
346 // m_HitDistance[gap] = dX;
347
348 // Now recalculate the total number of hits and the max. number of
349 // hits per gap.
350 // CountHits();
351}
352
353// Where does the trajectory of this track intersect a specific gap?
354void RecMucTrack::Project( const int& part, const int& gap, float& x, float& y, float& z,
355 int& seg ) {
356 seg = -1;
357 x = 0.0;
358 y = 0.0;
359 z = 0.0;
360
361 if ( ( gap < 0 ) || ( gap >= (int)MucID::getGapNum( part ) ) )
362 {
363 cout << "Muc2DRoad::Project(), invalid gap = " << gap << endl;
364 return;
365 }
366
367 HepPoint3D pos = m_CurrentPos;
368 Hep3Vector dir = m_CurrentDir;
369
370 vector<HepPoint3D> insctCon =
371 MucGeoGeneral::Instance()->FindIntersections( part, gap, pos, dir );
372 if ( insctCon.size() == 0 ) { return; }
373
374 HepPoint3D intersection = insctCon[0];
375
376 x = intersection.x();
377 y = intersection.y();
378 z = intersection.z();
379
380 // cout << " x " << x << " y " << y << " z " << z << endl;
381
382 float phi;
383 if ( ( x * x + y * y + z * z ) < kMinor ) { return; }
384 else
385 {
386 if ( part == 1 )
387 {
388 phi = acos( x / sqrt( x * x + y * y ) );
389 if ( y < 0 ) phi = 2.0 * kPi - phi;
390 phi = fmod( ( phi + kPi / 8.0 ), 2.0 * kPi );
391 seg = int( phi / ( kPi / 4.0 ) );
392 }
393 else
394 {
395 if ( x >= 0.0 )
396 {
397 if ( y >= 0.0 ) { seg = 0; }
398 else { seg = 3; }
399 }
400 else
401 {
402 if ( y >= 0.0 ) { seg = 1; }
403 else { seg = 2; }
404 }
405 }
406 }
407}
408
409float RecMucTrack::GetHitDistance( const MucRecHit* hit ) {
410 if ( !hit )
411 {
412 cout << "RecMucTrack:GetHitDistance-E1 null hit pointer!" << endl;
413 return kInfinity;
414 }
415
416 int part = hit->Part();
417 int gap = hit->Gap();
418 if ( ( gap < 0 ) || ( gap >= (int)MucID::getGapNum( part ) ) )
419 {
420 cout << "RecMucTrack::GetHitDistance() bad gap number = " << gap << endl;
421 return kInfinity;
422 }
423
424 HepPoint3D posHit = hit->GetCenterPos();
425 HepPoint3D posHitLocal = hit->GetGap()->TransformToGap( posHit );
426 HepPoint3D posInsctLocal = hit->GetGap()->TransformToGap( m_CurrentInsct );
427 int orient = hit->GetGap()->Orient();
428
429 float distance = -9990;
430 if ( orient == 1 ) distance = fabs( posInsctLocal.x() - posHitLocal.x() );
431 if ( orient == 0 ) distance = fabs( posInsctLocal.y() - posHitLocal.y() );
432
433 return distance;
434}
435
436// not abs value
437float RecMucTrack::GetHitDistance2( const MucRecHit* hit ) {
438 if ( !hit )
439 {
440 cout << "RecMucTrack:GetHitDistance-E1 null hit pointer!" << endl;
441 return kInfinity;
442 }
443
444 int part = hit->Part();
445 int gap = hit->Gap();
446 if ( ( gap < 0 ) || ( gap >= (int)MucID::getGapNum( part ) ) )
447 {
448 cout << "RecMucTrack::GetHitDistance() bad gap number = " << gap << endl;
449 return kInfinity;
450 }
451
452 HepPoint3D posHit = hit->GetCenterPos();
453 HepPoint3D posHitLocal = hit->GetGap()->TransformToGap( posHit );
454 HepPoint3D posInsctLocal = hit->GetGap()->TransformToGap( m_CurrentInsct );
455 int orient = hit->GetGap()->Orient();
456
457 float distance = -9990;
458 if ( orient == 1 ) distance = ( posInsctLocal.x() - posHitLocal.x() );
459 if ( orient == 0 ) distance = ( posInsctLocal.y() - posHitLocal.y() );
460
461 // cout<<"========in RecMucTrack: line insct: "<<posInsctLocal.x()<<" "<<posInsctLocal.y()<<"
462 // ->
463 // "<<posHitLocal.x()<<" "<<posHitLocal.y()<<endl;
464
465 return distance;
466}
467
468/// Calculate intersection from all hits attached on this gap.
469Hep3Vector RecMucTrack::CalculateInsct( const int part, const int seg, const int gap ) {
470 MucGeoGap* gapPtr = MucGeoGeneral::Instance()->GetGap( part, seg, gap );
471 vector<int> hitSeq;
472 for ( int i = 0; i < (int)m_pHits.size(); i++ )
473 {
474 MucRecHit* aHit = m_pHits[i];
475 if ( aHit->Part() == part && aHit->Seg() == seg && aHit->Gap() == gap )
476 { hitSeq.push_back( i ); }
477 }
478 int nHitInGap = hitSeq.size();
479 // cout << "nHitInGap " << nHitInGap << endl;
480
481 HepPoint3D insctLocal = gapPtr->TransformToGap( m_CurrentInsct );
482 // HepPoint3D newInsct(0.0, 0.0, 0.0);
483 HepPoint3D newInsctLocal = insctLocal;
484
485 vector<float> x;
486 for ( int i = 0; i < nHitInGap; i++ ) x.push_back( 0 );
487 float xInsct = 0, xNewInsct = 0;
488
489 int orient = gapPtr->Orient();
490 if ( orient == 1 ) xInsct = insctLocal.x();
491 if ( orient == 0 ) xInsct = insctLocal.y();
492
493 for ( int i = 0; i < nHitInGap; i++ )
494 {
495 float xStrip, yStrip, zStrip;
496 ( m_pHits[hitSeq[i]] )->GetStrip()->GetCenterPos( xStrip, yStrip, zStrip );
497 if ( orient == 1 ) x[i] = xStrip;
498 if ( orient == 0 ) x[i] = yStrip;
499 }
500 // cout << "local Old" << insctLocal << endl;
501
502 // if == 0, no direction change;
503 if ( nHitInGap > 0 )
504 {
505 xNewInsct = xInsct * kInsctWeight;
506
507 // float minDist = kInfinity;
508 for ( int i = 0; i < nHitInGap; i++ )
509 {
510 // if(fabs(x[i] - xInsct) < minDist) {
511 // xNewInsct = x[i]; //}
512 xNewInsct += x[i] * ( ( 1.0 - kInsctWeight ) / nHitInGap );
513 }
514
515 if ( orient == 1 )
516 {
517 newInsctLocal.setX( xNewInsct );
518 newInsctLocal.setY( insctLocal.y() );
519 }
520 if ( orient == 0 )
521 {
522 newInsctLocal.setX( insctLocal.x() );
523 newInsctLocal.setY( xNewInsct );
524 }
525 }
526
527 m_CurrentInsct = gapPtr->TransformToGlobal( newInsctLocal );
528 // cout << "local New" << newInsctLocal << endl;
529 // cout << "global New" << m_CurrentInsct << endl;
530
531 return m_CurrentInsct;
532}
533
534void RecMucTrack::AttachInsct( Hep3Vector insct ) { m_Intersections.push_back( insct ); }
535
536void RecMucTrack::AttachDirection( Hep3Vector dir ) { m_Directions.push_back( dir ); }
537
538void RecMucTrack::CorrectDir() {
539 // cout << "Before CorrectDir(), fCurrentDir " << m_CurrentDir << endl;
540 m_CurrentDir = m_CurrentInsct - m_CurrentPos;
541 m_CurrentDir.setMag( 1.0 );
542 // cout << "After CorrectDir(), fCurrentDir " << m_CurrentDir << endl;
543}
544
545void RecMucTrack::CorrectPos() { m_CurrentPos = m_CurrentInsct; }
546
547void RecMucTrack::ComputeLastGap() {
548 int lastGap1, lastGap2;
549 int firstGap1, firstGap2;
550 lastGap1 = lastGap2 = -1;
551 firstGap1 = firstGap2 = 99;
552 vector<MucRecHit*>::const_iterator iHit;
553 iHit = m_pHits.begin();
554 for ( ; iHit != m_pHits.end(); iHit++ )
555 {
556 if ( *iHit )
557 { // Check for a null pointer.
558 int part = ( *iHit )->Part();
559 int gap = ( *iHit )->Gap();
560
561 if ( part == 1 )
562 {
563 if ( gap > lastGap1 ) lastGap1 = gap;
564 if ( gap < firstGap1 ) firstGap1 = gap;
565 }
566 else
567 {
568 if ( gap > lastGap2 )
569 {
570 m_ecPart = part;
571 m_endPart = part;
572 lastGap2 = gap;
573 }
574 if ( gap < firstGap2 ) firstGap2 = gap;
575 }
576 }
577 }
578
579 m_brLastLayer = lastGap1;
580 if ( firstGap1 == 99 ) m_brFirstLayer = -1;
581 else if ( firstGap1 >= 0 && firstGap1 < 9 ) m_brFirstLayer = firstGap1;
582
583 m_ecLastLayer = lastGap2;
584 if ( firstGap2 == 99 ) m_ecFirstLayer = -1;
585 else if ( firstGap2 >= 0 && firstGap2 < 8 ) m_ecFirstLayer = firstGap2;
586
587 // cout<<"MucTrack, br: "<<m_brFirstLayer<<", "<<m_brLastLayer<<"\tec: "<<m_ecFirstLayer<<",
588 // "<<m_ecLastLayer<<endl;
589}
590
591void RecMucTrack::ComputeDepth( int method ) {
592 if ( m_numLayers == 0 )
593 {
594 m_depth = m_depth_3 = -99;
595 return;
596 }
597 else if ( m_numLayers == 1 && ( m_brLastLayer == 0 || m_ecLastLayer == 0 ) )
598 {
599 m_depth = m_depth_3 = 0;
600 return;
601 }
602
603 m_depth_3 = 0.0;
604
605 float brThickness = 0.0;
606 float ecThickness = 0.0;
607 float deltaPhi = 0.0;
608 int betweenSeg = 0;
609
610 float phi = m_MucMomentum.phi();
611 float theta = m_MucMomentum.theta();
612
613 vector<MucRecHit*>::const_iterator iHit;
614 int ngap = MucID::getGapMax();
615 // cout<<"ngap:\t"<<ngap<<endl;
616
617 int Seg[9];
618 int part, seg, gap, strip;
619 for ( int gap = 0; gap < ngap; gap++ ) { Seg[gap] = -1; }
620 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ )
621 {
622 if ( *iHit )
623 { // Check for a null pointer.
624 part = ( *iHit )->Part();
625 seg = ( *iHit )->Seg();
626 gap = ( *iHit )->Gap();
627 strip = ( *iHit )->Strip();
628 if ( part == 1 ) Seg[gap] = seg;
629 }
630 }
631
632 int segTrackBr = -1;
633 for ( int gap = 0; gap <= brLastLayer(); gap++ )
634 {
635 if ( Seg[gap] != -1 ) segTrackBr = Seg[gap];
636 }
637 // BR, the innermost layer is RPC module
638 for ( int gap = 0; gap <= brLastLayer(); gap++ )
639 {
640 float thickness = 0.0;
641 if ( Seg[gap] != -1 && Seg[brLastLayer() - 1] != -1 && Seg[gap] != Seg[brLastLayer() - 1] )
642 betweenSeg = 1;
643 thickness = MucGeoGeneral::Instance()->GetGap( 1, 0, gap )->GetIronThickness();
644 // cout<<"RecMucTrack gap="<<gap<<" brlastgap="<<brLastLayer()<<" "<<thickness<<endl;
645 if ( sin( m_MucMomentum.theta() ) != 0 ) thickness /= sin( m_MucMomentum.theta() );
646 else
647 ; // cout<<"RecMucTrack::ComputeDepth, In Barrel,theta=0?"<<endl;
648 deltaPhi = m_MucMomentum.phi() - segTrackBr * kPi / 4;
649 if ( Seg[gap] == -1 && betweenSeg == 1 )
650 {
651 thickness += 40; // gap width
652 }
653 if ( cos( deltaPhi ) != 0 ) thickness /= cos( deltaPhi );
654 else
655 ; // cout<<"RecMucTrack::ComputeDepth, In Barrel,Cos(phi)=0?"<<endl;
656 if ( deltaPhi == 0 && Seg[brLastLayer() - 1] == 2 ) thickness = 0;
657 // cout<<"in muctrack "<<thickness<<" "<<brThickness<<" theta "<<m_MucMomentum.theta()<<"
658 // phi="<<deltaPhi<<"
659 // "<<m_MucMomentum.phi()<<endl;
660
661 if ( brFirstLayer() < brLastLayer() ) brThickness += thickness;
662 else if ( brFirstLayer() == brLastLayer() )
663 { // only one gap
664 if ( m_MucMomentum.mag() > 1000 || brFirstLayer() < 2 )
665 brThickness += thickness; // high momentum or only one or two gap
666 // cout<<"mom="<<m_MucMomentum.mag()<<" "<<brFirstLayer()<<" "<<brThickness<<"
667 // "<<thickness<<endl;
668 }
669 else cout << "in RecMucTrack: Wrong Gap Info" << endl;
670
671 // cout<<brThickness<<endl;
672 }
673
674 // cout<<"eclastgap= "<<ecLastLayer()<<" ecfirstgap= "<<ecFirstLayer()<<endl;
675
676 // EC, the innermost layer is Iron
677 // for (int gap = ecFirstLayer(); gap!=-1&&gap <= ecLastLayer(); gap++) {
678 for ( int gap = 0; gap != -1 && gap <= ecLastLayer(); gap++ )
679 { ecThickness += MucGeoGeneral::Instance()->GetGap( 0, 0, gap )->GetIronThickness(); }
680
681 if ( cos( theta ) != 0.0 ) ecThickness /= cos( theta );
682 else
683 ; // cout << "RecMucTrack::ComputeDepth, In EndCap, Track theta = 90.0 ? " << endl;
684 ecThickness = fabs( ecThickness );
685
686 // cout<<"eclastgap= "<<ecLastLayer()<<" ecthickness="<<ecThickness<<endl;
687
688 if ( method == 2 )
689 {
690 // barrel first
691 if ( ( m_Good3DLine == 1 ) && ( m_Good3DPart == 1 ) )
692 {
693 if ( m_IntersectionInner[0].x() != -9999 )
694 {
695 for ( int gap = 1; gap <= brLastLayer(); gap++ ) // Inner[gap1]-Outer[gap0] ...
696 {
697 if ( m_IntersectionInner[gap].x() != -9999 &&
698 m_IntersectionInner[gap - 1].x() != -9999 )
699 m_depth_3 += ( m_IntersectionInner[gap] - m_IntersectionOuter[gap - 1] ).mag();
700 }
701 // may be pass some gap in endcap!
702 m_depth_3 += ecThickness;
703 }
704 else m_depth_3 = brThickness + ecThickness;
705 }
706 if ( ( m_Good3DLine == 1 ) && ( m_Good3DPart != 1 ) )
707 {
708 for ( int gap = 1; gap <= ecLastLayer(); gap++ ) // Inner[gap1]-Outer[gap0] ...
709 {
710 if ( m_IntersectionInner[gap].x() != -9999 &&
711 m_IntersectionInner[gap - 1].x() != -9999 )
712 m_depth_3 += ( m_IntersectionInner[gap] - m_IntersectionOuter[gap - 1] ).mag();
713 }
714 // may be pass some gap in barrel!
715 m_depth_3 += brThickness;
716 if ( cos( theta ) != 0.0 )
717 m_depth_3 += 40 / cos( theta ); // there is 1 absorber before first gap in endcap!
718 }
719 if ( m_Good3DLine == 0 ) m_depth_3 = brThickness + ecThickness;
720 if ( m_depth_3 > 2000 )
721 m_depth_3 = brThickness + ecThickness; // unreasonable depth! so use previous method!
722 }
723 else // method == 1 linefit
724 { m_depth_3 = brThickness + ecThickness; }
725
726 double offset = 50.0;
727 // if(GetTotalHits() > 0) m_depth_3 += offset;
728
729 m_depth = m_depth_3;
730
731 // if(m_depth<0||m_depth>2000) m_depth = 0;
732 if ( m_depth > 2000 ) m_depth = -99;
733 // cout<<"depth= "<<m_depth<<endl;
734}
735
736void RecMucTrack::ComputeDepth() {
737 m_depth = -99.0;
738 // Part 1 first.
739 float brThickness = 0.0;
740 for ( int gap = 0; gap <= brLastLayer(); gap++ )
741 { brThickness += MucGeoGeneral::Instance()->GetGap( 1, 0, gap )->GetIronThickness(); }
742
743 // second alg
744 float brThickness_2 = 0.0;
745 for ( int gap = 0; gap <= brLastLayer(); gap++ )
746 {
747 if ( HasHitInGap( 1, gap ) )
748 { brThickness_2 += MucGeoGeneral::Instance()->GetGap( 1, 0, gap )->GetIronThickness(); }
749 }
750 // third alg
751 float brThickness_3 = 0.0;
752 vector<MucRecHit*>::const_iterator iHit;
753 int ngap = MucID::getGapMax();
754 int Seg[9];
755 int part, seg, gap, strip;
756 for ( int gap = 0; gap < ngap; gap++ ) { Seg[gap] = -1; }
757 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ )
758 {
759 if ( *iHit )
760 { // Check for a null pointer.
761 part = ( *iHit )->Part();
762 seg = ( *iHit )->Seg();
763 gap = ( *iHit )->Gap();
764 strip = ( *iHit )->Strip();
765 if ( part == 1 ) Seg[gap] = seg;
766 }
767 }
768
769 float deltaPhi_3 = 0.0;
770 int betweenSeg = 0;
771
772 int segTrackBr = -1;
773 for ( int gap = 0; gap <= brLastLayer(); gap++ )
774 {
775 if ( Seg[gap] != -1 ) segTrackBr = Seg[gap];
776 }
777
778 for ( int gap = 0; gap <= brLastLayer(); gap++ )
779 {
780 float thickness = 0.0;
781 if ( Seg[gap] != -1 && Seg[brLastLayer() - 1] != -1 && Seg[gap] != Seg[brLastLayer() - 1] )
782 betweenSeg = 1;
783 thickness = MucGeoGeneral::Instance()->GetGap( 1, 0, gap )->GetIronThickness();
784 if ( sin( m_MucMomentum.theta() ) != 0 ) thickness /= sin( m_MucMomentum.theta() );
785 else cout << "RecMucTrack::ComputeDepth, In Barrel,theta=0?" << endl;
786 // if(Seg[gap] != -1) deltaPhi_3 = m_MucMomentum.phi() - Seg[gap]*kPi/4;
787 deltaPhi_3 = m_MucMomentum.phi() -
788 segTrackBr * kPi / 4; // some times, no hits in a gap, but a good track exist!
789 if ( Seg[gap] == -1 && betweenSeg == 1 )
790 {
791 cout << "between segment" << endl;
792 thickness += 40; // gap width
793 }
794
795 if ( cos( deltaPhi_3 ) != 0 ) thickness /= cos( deltaPhi_3 );
796 else cout << "RecMucTrack::ComputeDepth, In Barrel,Cos(phi)=0?" << endl;
797
798 if ( deltaPhi_3 == 0 && Seg[brLastLayer() - 1] == 2 ) thickness = 0;
799 // cout<<"in muctrack "<<thickness<<" "<<brThickness_3<<" theta "<<m_MucMomentum.theta()<<"
800 // phi="<<deltaPhi_3<<"
801 // "<<m_MucMomentum.phi()<<endl;
802 brThickness_3 += thickness;
803 }
804
805 // cout<<"in RecMucTrack: compare thickness "<<brThickness<<" "<<brThickness_2<<"
806 // "<<brThickness_3<<endl;
807
808 float phi = m_MucMomentum.phi();
809 float deltaPhi = phi - kPi / 4 * (int)( phi / ( kPi / 4 ) );
810 if ( deltaPhi > kPi / 8 ) deltaPhi = kPi / 4 - deltaPhi;
811 float theta = m_MucMomentum.theta();
812 // cout << "br LastGap " << brLastLayer() << " Thick " << brThickness
813 // << " 1/sin(theta) " << 1/sin(theta) << " 1/cos(deltaPhi) " << 1/cos(deltaPhi) << endl;
814
815 if ( sin( theta ) != 0.0 ) brThickness /= sin( theta );
816 else cout << "RecMucTrack::ComputeDepth, In Barrel, Track theta = 0.0 ? " << endl;
817
818 brThickness /= cos( deltaPhi );
819 brThickness = fabs( brThickness );
820 // cout << "br Depth " << brThickness << endl;
821
822 // EC, then
823 float ecThickness = 0.0;
824 for ( int gap = 0; gap <= ecLastLayer(); gap++ )
825 { ecThickness += MucGeoGeneral::Instance()->GetGap( 0, 0, gap )->GetIronThickness(); }
826 // cout << "ec LastGap " << ecLastLayer() << " Thick " << ecThickness
827 // << " 1/cos(theta) " << 1/cos(theta) << endl;
828
829 if ( cos( theta ) != 0.0 ) ecThickness /= cos( theta );
830 else
831 ; // cout << "RecMucTrack::ComputeDepth, In EndCap, Track theta = 90.0 ? " << endl;
832 ecThickness = fabs( ecThickness );
833 // cout << "ec Depth " << ecThickness << endl;
834
835 m_depth = brThickness + ecThickness;
836 m_depth_3 = brThickness_3 + ecThickness;
837 cout << "Depth " << m_depth << " Depth3 = " << m_depth_3 << endl;
838 m_depth = m_depth_3;
839 double offset = 50.0;
840 if ( GetTotalHits() > 0 )
841 m_depth +=
842 offset; // since brThickness on gap 0 is zero, give an offset for track arriving muc.
843}
844
845void RecMucTrack::ComputeDistanceMatch() {
846 bool firstHitFound = false;
847 MucGeoGap* firstGap = 0;
848 vector<MucRecHit*>::const_iterator iHit;
849 vector<MucRecHit*> hitsGap0;
850 float stripLocal[3] = { 0.0, 0.0, 0.0 };
851 float stripGlobal[3] = { 0.0, 0.0, 0.0 };
852 int nStrip = 0;
853
854 int part, seg, gap, strip;
855 int barrel_gap0_exist = 0;
856
857 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ )
858 {
859 if ( *iHit )
860 { // Check for a null pointer.
861 part = ( *iHit )->Part();
862 seg = ( *iHit )->Seg();
863 gap = ( *iHit )->Gap();
864 strip = ( *iHit )->Strip();
865 if ( !firstHitFound && gap == 0 )
866 {
867 firstGap = MucGeoGeneral::Instance()->GetGap( part, seg, gap );
868 firstHitFound = true;
869 }
870 if ( firstGap && part == firstGap->Part() && seg == firstGap->Seg() &&
871 gap == firstGap->Gap() )
872 {
873 // cout<<"in RecMucTrack "<<part<<" "<<seg<<" "<<gap<<" "<<strip<<endl;
874 HepPoint3D posHit = ( *iHit )->GetCenterPos();
875 HepPoint3D posHitLocal = ( *iHit )->GetGap()->TransformToGap( posHit );
876 if ( part == 1 && gap == 0 ) barrel_gap0_exist = 1; // exist
877
878 stripLocal[0] += posHitLocal.x();
879 stripLocal[1] += posHitLocal.y();
880 stripLocal[2] += posHitLocal.z();
881
882 stripGlobal[0] += posHit.x(); // to calc phi of this strip
883 stripGlobal[1] += posHit.y();
884 stripGlobal[2] += posHit.z();
885
886 nStrip++;
887 }
888 }
889 }
890
891 // cout<<"in RecMucTrack: extpos "<<m_ExtMucPos<<" mucpos "<< m_MucPos<<endl;
892
893 int apart = -1, aseg = -1;
894 int nHits = FindSegWithMaxHits( apart, aseg );
895 MucGeoGap* fakefirstGap = 0; // maybe not exist!
896 if ( apart == -1 && aseg == -1 ) { m_Dist_muc_ext.set( 0, 0, 0 ); }
897 else
898 {
899 fakefirstGap = MucGeoGeneral::Instance()->GetGap( apart, aseg, 0 );
900 HepPoint3D fextLocal = fakefirstGap->TransformToGap( m_ExtMucPos );
901 HepPoint3D fmucLocal = fakefirstGap->TransformToGap( m_MucPos );
902 float dist_x = fextLocal.x() - fmucLocal.x();
903 float dist_y = fextLocal.y() - fmucLocal.y();
904 float dist_z = fextLocal.z() - fmucLocal.z();
905 m_Dist_muc_ext.set( dist_x, dist_y, dist_z );
906 // cout<<"in RecMucTrack dist = "<<dist_x<<" "<<dist_y<<" "<<dist_z<<endl;
907 if ( fakefirstGap->Orient() == 0 )
908 { // part 0,2
909 // cout<<"in RecMucTrack "<< extLocal.y()<<" "<<stripLocal[1]<<endl;
910 }
911 else
912 {
913 // cout<<"in RecMucTrack1 "<< extLocal.x()<<" "<<stripLocal[0]<<endl;
914 }
915 }
916
917 float distance = -9990;
918
919 if ( nStrip == 0 || !firstGap )
920 {
921 distance = -9990;
922 // cout << "no hits on first gap" << endl;
923 }
924 else
925 {
926 for ( int k = 0; k < 3; k++ ) stripLocal[k] /= nStrip;
927 for ( int k = 0; k < 3; k++ ) stripGlobal[k] /= nStrip;
928
929 m_StripPhi.set( stripGlobal[0], stripGlobal[1], stripGlobal[2] );
930
931 HepPoint3D extLocal = firstGap->TransformToGap( m_ExtMucPos );
932 // extmucpos may be can be used to identify mu/pion
933
934 // cout<<"in RecMucTrack mom "<<m_MucMomentum.x()<<" "<<m_MucMomentum.y()<<"
935 // "<<m_MucMomentum.z()<<endl; cout<<"in RecMucTrack extpos "<<m_ExtMucPos.x()<<"
936 // "<<m_ExtMucPos.y()<<" "<<m_ExtMucPos.z()<<endl; cout<<"in RecMucTrack extpos2
937 // "<<ExtMucPos_2.x()<<" "<<ExtMucPos_2.y()<<" "<<ExtMucPos_2.z()<<endl; cout<<"in
938 // RecMucTrack extloc "<<extLocal.x()<<" "<<extLocal.y()<<" "<<extLocal.z()<<endl;
939 if ( firstGap->Orient() == 0 )
940 { // part 0,2
941 distance = extLocal.y() - stripLocal[1];
942 // cout<<"in RecMucTrack "<< extLocal.y()<<" "<<stripLocal[1]<<endl;
943 }
944 else
945 {
946 distance = extLocal.x() - stripLocal[0];
947 // cout<<"in RecMucTrack1 "<< extLocal.x()<<" "<<stripLocal[0]<<endl;
948 }
949 }
950
951 m_distance = distance;
952
953 // use m_chi2 temporary
954 // m_chi2 = distance;
955 // cout<<"in RecMucTrack distance= "<<m_distance<<" n= "<<nStrip<<endl;
956}
957
958void RecMucTrack::Extend() {
959 float x = m_ExtMucPos.x(), y = m_ExtMucPos.y(), z = m_ExtMucPos.z();
960 // cout<<"in MucTrac, MucPos= "<<x<<" "<<y<<" "<<z<<endl;
961 float step = 0.005;
962 float myMucR = 1720.0;
963 float myMucZ = 2110.0;
964 int count_ = 0; // added by LI Chunhua to avoid loop infinitely
965 while ( ( ( fabs( x ) <= myMucR ) && ( fabs( y ) <= myMucR ) &&
966 ( ( fabs( x - y ) / sqrt( 2. ) ) <= myMucR ) &&
967 ( ( fabs( x + y ) / sqrt( 2. ) ) <= myMucR ) && ( fabs( z ) <= myMucZ ) ) &&
968 count_ < 1000 )
969 {
970 x += step * m_MucMomentum.x();
971 y += step * m_MucMomentum.y();
972 z += step * m_MucMomentum.z();
973 count_++;
974 // cout << "in RecMucTrack+ x " << x << " y " << y << " z " << z <<" mom "<<
975 // m_MucMomentum.x()<< "
976 // "<<m_MucMomentum.y()<<" "<<m_MucMomentum.z()<<endl;
977 }
978
979 int count = 0;
980 while ( !( ( fabs( x ) <= myMucR ) && ( fabs( y ) <= myMucR ) &&
981 ( ( fabs( x - y ) / sqrt( 2. ) ) <= myMucR ) &&
982 ( ( fabs( x + y ) / sqrt( 2. ) ) <= myMucR ) && ( fabs( z ) <= myMucZ ) ) &&
983 count < 1000 )
984 {
985 x -= step * m_MucMomentum.x();
986 y -= step * m_MucMomentum.y();
987 z -= step * m_MucMomentum.z();
988 count++;
989 // cout << "in RecMucTrack- x " << x << " y " << y << " z " << z <<" mom "<<
990 // m_MucMomentum.x()<< "
991 // "<<m_MucMomentum.y()<<" "<<m_MucMomentum.z()<<endl;
992 }
993
994 if ( count < 1000 && count_ < 1000 )
995 {
996 if ( fabs( x ) < 2600 && fabs( y ) < 2600 && fabs( z ) < 2800 )
997 {
998 m_ExtMucPos.set( x, y, z );
999 m_MucPos.set( x, y, z );
1000 m_xPos = x;
1001 m_yPos = y;
1002 m_zPos = z;
1003 }
1004 }
1005}
1006
1007void RecMucTrack::LineFit( int fittingMethod ) {
1008 Extend();
1009
1010 int part = -1, seg = -1;
1011 int nHits = FindSegWithMaxHits( part, seg );
1012 // cout << "RecMucTrack::ComputeDepth(), part " << part << " seg " << seg
1013 // << " contains most hits " << nHits << endl;
1014 if ( part < 0 || seg < 0 )
1015 {
1016 m_depth = 0;
1017 // cout << "No hit in this track" << endl;
1018 return;
1019 }
1020
1021 float startPos[3] = { 0.0, 0.0, 0.0 };
1022 MucRec2DRoad road2D[2];
1023 vector<MucRecHit*>::const_iterator iHit;
1024 for ( int orient = 0; orient <= 1; orient++ )
1025 {
1026 road2D[orient] = MucRec2DRoad( part, seg, orient, startPos[0], startPos[1], startPos[2],
1027 fittingMethod );
1028 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ )
1029 {
1030 if ( *iHit )
1031 { // Check for a null pointer.
1032 int hitPart = ( *iHit )->Part();
1033 int hitSeg = ( *iHit )->Seg();
1034 int hitOrient = ( *iHit )->GetGap()->Orient();
1035 if ( hitPart == part && hitSeg == seg && hitOrient == orient )
1036 { road2D[orient].AttachHit( *iHit ); }
1037 }
1038 }
1039 // cout << "Orient " << orient
1040 // << " Slope " << road2D[orient].GetSlope()
1041 // << " Intercept " << road2D[orient].GetIntercept()
1042 // << " LastGap " << road2D[orient].GetLastGap()
1043 // << endl;
1044 }
1045 MucRec3DRoad road3D( &road2D[0], &road2D[1] );
1046 float vx, vy, x0, y0;
1047 if ( road3D.GetPart() == 1 )
1048 {
1049 road3D.TransformPhiRToXY( road2D[1].GetSlope(), road2D[0].GetSlope(),
1050 road2D[1].GetIntercept(), road2D[0].GetIntercept(), vx, vy, x0,
1051 y0 );
1052 }
1053 else
1054 {
1055 vx = road2D[1].GetSlope();
1056 x0 = road2D[1].GetIntercept();
1057 vy = road2D[0].GetSlope();
1058 y0 = road2D[0].GetIntercept();
1059 }
1060
1061 // cout << "road3D Last Gap " << road3D.GetLastGap() << endl;
1062 // cout << "vx " << vx << " x0 " << x0 << " vy " << vy << " y0 " << y0 << endl;
1063
1064 // if number of gaps with hits >= 2 on both orient, change muc pos and momentum
1065 if ( road2D[0].GetNGapsWithHits() >= 2 && road2D[1].GetNGapsWithHits() >= 2 )
1066 {
1067
1068 // good 3d road!!!
1069 m_Good3DLine = 1;
1070 m_status = 1;
1071
1072 m_Good3DPart = road3D.GetPart();
1073 road3D.SetSimpleFitParams( vx, x0, vy, y0 );
1074 // road3D.PrintHitsInfo();
1075
1076 float startx = 0.0, starty = 0.0, startz = 0.0;
1077 float startxSigma = 0.0, startySigma = 0.0, startzSigma = 0.0;
1078 float x1 = 0.0, y1 = 0.0, z1 = 0.0, x2 = 0.0, y2 = 0.0, z2 = 0.0;
1079 int gap = 0;
1080
1081 if ( fittingMethod == 2 )
1082 { // if choose quadratic fitting method!
1083 for ( int igap = 0; igap < 9; igap++ )
1084 {
1085 road3D.Project( igap, x1, y1, z1, x2, y2, z2 );
1086 m_IntersectionInner[igap].set( x1, y1, z1 );
1087 m_IntersectionOuter[igap].set( x2, y2, z2 );
1088 // cout<<"3dproject sur "<<x1<<" "<<y1<<" "<<z1<<" "<<x2<<" "<<y2<<" "<<z2<<endl;
1089 }
1090 }
1091
1092 do {
1093 // road3D.Project(gap, startx, starty, startz);
1094 road3D.ProjectWithSigma( gap, startx, starty, startz, startxSigma, startySigma,
1095 startzSigma );
1096 gap++;
1097 } while ( ( startx * startx + starty * starty + startz * startz ) < kMinor &&
1098 gap < (int)MucID::getGapNum( part ) );
1099
1100 if ( fabs( startx ) < 2600 && fabs( starty ) < 2600 && fabs( startz ) < 2800 )
1101 {
1102 Hep3Vector MucPos_self;
1103 MucPos_self.set( startx, starty, startz );
1104 float dist = ( MucPos_self - m_MucPos ).mag();
1105
1106 if ( dist < 1000 )
1107 { // (mm) maybe the fit is bad, if the dist is too big.
1108 m_MucPos.set( startx, starty, startz );
1109 m_xPos = startx;
1110 m_yPos = starty;
1111 m_zPos = startz;
1112 }
1113 }
1114 m_MucPosSigma.set( startxSigma, startySigma, startzSigma );
1115 m_xPosSigma = startxSigma;
1116 m_yPosSigma = startySigma;
1117 m_zPosSigma = startzSigma;
1118
1119 // cout<<"in RecMucTrack gap= "<<gap<<" start= "<<startx<<" "<<starty<<" "<<startz<<"
1120 // "<<endl;
1121 float momentum = m_MucMomentum.mag();
1122 float zDir = 1.0;
1123 if ( m_MucMomentum.z() != 0.0 ) zDir = m_MucMomentum.z() / fabs( m_MucMomentum.z() );
1124
1125 float px = road3D.GetSlopeZX() * zDir;
1126 float py = road3D.GetSlopeZY() * zDir;
1127 float pz = zDir;
1128 float segPhi = 0.25 * kPi * seg;
1129 // if vr is opposite to original, but both are very small, e.x. -0.01 -> + 0.01, or you can
1130 // say, pt~p, change it to guarantee vx, vy is correct.
1131 if ( part == 1 && px * cos( segPhi ) + py * sin( segPhi ) < 0.0 )
1132 { // when in barrel, pt dot segPhi < 0, angle between them > 90deg
1133 px *= -1.0;
1134 py *= -1.0;
1135 pz *= -1.0;
1136 }
1137
1138 // if(sqrt(px*px+py*py)>0.01)m_MucMomentum.set(px, py, pz);
1139 m_MucMomentum.set( px, py, pz );
1140 m_MucMomentum.setMag( momentum );
1141
1142 m_px = px;
1143 m_py = py;
1144 m_pz = pz;
1145 Hep3Vector phi_mdc, phi_muc;
1146 phi_mdc.set( m_MdcMomentum.x(), m_MdcMomentum.y(), 0 );
1147 phi_muc.set( px, py, 0 );
1148 double deltaPhi = phi_mdc.angle( phi_muc );
1149
1150 m_deltaPhi = deltaPhi;
1151 m_dof = road3D.GetDegreesOfFreedom();
1152 m_chi2 = road3D.GetReducedChiSquare();
1153 if ( m_chi2 < 0 || m_chi2 > 1000 ) m_chi2 = 0;
1154 m_rms = 0.0; // road3D.GetReducedChiSquare(); // ??? in MucRecRoad3D
1155
1156 if ( road2D[0].GetNGapsWithHits() >= 3 && road2D[1].GetNGapsWithHits() >= 3 )
1157 { // need 3+ gap to do refit
1158 // calc distance between each hit with this track.
1159 // use GetHitDistance(), so we should provide m_CurrentInsct.
1160 float xInsct = 0.0, yInsct = 0.0, zInsct = 0.0, sigmax = 0.0, sigmay = 0.0, sigmaz = 0.0;
1161 float quad_x1 = 0.0, quad_y1 = 0.0, quad_z1 = 0.0, quad_x2 = 0.0, quad_y2 = 0.0,
1162 quad_z2 = 0.0;
1163 for ( int ihit = 0; ihit < m_pHits.size(); ihit++ )
1164 {
1165 int igap = ( m_pHits[ihit] )->Gap();
1166 road3D.ProjectNoCurrentGap( igap, xInsct, yInsct, zInsct, sigmax, sigmay, sigmaz,
1167 quad_x1, quad_y1, quad_z1, quad_x2, quad_y2, quad_z2 );
1168
1169 SetCurrentInsct( xInsct, yInsct, zInsct );
1170
1171 float dist_hit_track = GetHitDistance2( m_pHits[ihit] );
1172 m_distHits.push_back( dist_hit_track );
1173
1174 // cout<<"===========in RecMucTrack: line dist: "<<dist_hit_track<<endl;
1175
1176 if ( fittingMethod == 2 ) // ------calc local insct in a gap with quad line.
1177 {
1178
1179 Hep3Vector center = m_pHits[ihit]->GetCenterPos();
1180 int iPart = m_pHits[ihit]->Part();
1181 int iSeg = m_pHits[ihit]->Seg();
1182 int iGap = m_pHits[ihit]->Gap();
1183 float xHit, yHit, zHit = 0.;
1184 if ( iPart == 1 )
1185 {
1186 if ( iGap % 2 == 1 )
1187 {
1188 xHit = center.z();
1189 yHit = sqrt( center.x() * center.x() + center.y() * center.y() );
1190 if ( iSeg == 2 )
1191 yHit = center.y(); // deal with seg2 seperately! because there is a hole here.
1192 // 2006.11.9
1193 }
1194 else
1195 {
1196 xHit = center.x();
1197 yHit = center.y();
1198 }
1199 }
1200 else
1201 {
1202 if ( iGap % 2 == 0 )
1203 {
1204 xHit = center.z();
1205 yHit = center.y();
1206 }
1207 else
1208 {
1209 xHit = center.z();
1210 yHit = center.x();
1211 }
1212 }
1213
1214 float distance1 = fabs( xHit - quad_x1 ) / ( xHit - quad_x1 ) *
1215 sqrt( ( xHit - quad_x1 ) * ( xHit - quad_x1 ) +
1216 ( yHit - quad_y1 ) * ( yHit - quad_y1 ) );
1217 float distance2 = fabs( xHit - quad_x2 ) / ( xHit - quad_x2 ) *
1218 sqrt( ( xHit - quad_x2 ) * ( xHit - quad_x2 ) +
1219 ( yHit - quad_y2 ) * ( yHit - quad_y2 ) );
1220
1221 float dist_quad = distance1;
1222 if ( fabs( distance1 ) > fabs( distance2 ) ) dist_quad = distance2;
1223
1224 // cout<<"============in RecMucTrack: quad xyz: "<<quad_x1<<" "<<quad_y1<<"
1225 // "<<quad_z1<<" "<<quad_x2<<"
1226 // "<<quad_y2<<endl; cout<<"============in RecMucTrack: hit pos: "<<xHit<<"
1227 // "<<yHit<<" "<<zHit<<endl; cout<<"============in RecMucTrack: dist1 =
1228 // "<<distance1<<" dist2 = "<<distance2<<" dist
1229 // ="<<dist_quad<<endl;
1230
1231 if ( quad_x1 == -9999 ) m_distHits_quad.push_back( -99 );
1232 else m_distHits_quad.push_back( dist_quad );
1233 }
1234 }
1235 }
1236 else
1237 {
1238 for ( int ihit = 0; ihit < m_pHits.size(); ihit++ )
1239 {
1240 m_distHits.push_back( -99 );
1241 m_distHits_quad.push_back( -99 );
1242 }
1243 }
1244 //////////////////////////////////////
1245
1246 ///////////////find expected strips for this 3D road
1247 ComputeLastGap();
1248 // HepPoint3D initPos(startx, starty,startz);
1249 HepPoint3D initPos( startx - m_MucMomentum.x() / momentum,
1250 starty - m_MucMomentum.y() / momentum,
1251 startz - m_MucMomentum.z() / momentum );
1252
1253 // for(int igap = 0; igap <= m_brLastLayer; igap++){
1254 for ( int igap = 0; igap < 9; igap++ )
1255 {
1256 vector<int> padID;
1257 vector<float> intersect_x;
1258 vector<float> intersect_y;
1259 vector<float> intersect_z;
1260
1261 // refit---------------------------------
1262 float zx, zy, x0, y0;
1263 int status = road3D.RefitNoCurrentGap( igap, zx, x0, zy, y0 );
1264 Hep3Vector mom_refit;
1265 if ( status == 0 )
1266 { // refit succeed
1267 float zDir = 1.0;
1268 if ( m_MucMomentum.z() != 0.0 ) zDir = m_MucMomentum.z() / fabs( m_MucMomentum.z() );
1269 float px_refit = zx * zDir;
1270 float py_refit = zy * zDir;
1271 float pz_refit = zDir;
1272 float segPhi = 0.25 * kPi * seg;
1273 // if vr is opposite to original, but both are very small, e.x. -0.01 -> + 0.01, or you
1274 // can say, pt~p, change it to guarantee vx, vy is correct.
1275 if ( part == 1 && px_refit * cos( segPhi ) + py_refit * sin( segPhi ) < 0.0 )
1276 { // when in barrel, pt dot segPhi < 0, angle between them > 90deg
1277 px_refit *= -1.0;
1278 py_refit *= -1.0;
1279 pz_refit *= -1.0;
1280 }
1281
1282 mom_refit.setX( px_refit );
1283 mom_refit.setY( py_refit );
1284 mom_refit.setZ( pz_refit );
1285 mom_refit.setMag( momentum );
1286 }
1287 else mom_refit = m_MucMomentum; // use the former momentum
1288
1289 HepPoint3D initPos_refit;
1290 if ( status == 0 )
1291 {
1292 float initPosx_refit, initPosy_refit, initPosz_refit;
1293 float sigmax_refit, sigmay_refit, sigmaz_refit;
1294 road3D.Project( 0, zx, x0, zy, y0, initPosx_refit, initPosy_refit, initPosz_refit );
1295 initPos_refit.setX( initPosx_refit - mom_refit.x() / momentum );
1296 initPos_refit.setY( initPosy_refit - mom_refit.y() / momentum );
1297 initPos_refit.setZ( initPosz_refit - mom_refit.z() / momentum );
1298 }
1299 else initPos_refit = initPos;
1300
1301 // cout<<"initPos: "<<initPos<<" "<<initPos_refit<<endl;
1302 // if((initPos - initPos_refit).mag()>100) cout<<"--------=====------to far"<<endl;
1303 // refit---------------------------------
1304 // cout<<"no gap: "<<igap<<" mom: "<<m_MucMomentum<<" refit: "<<mom_refit<<endl;
1305
1306 // vector<Identifier> MuId =
1307 // road3D.ProjectToStrip(1,igap,initPos,m_MucMomentum,padID,intersect_x,intersect_y,intersect_z);
1308 vector<Identifier> MuId = road3D.ProjectToStrip(
1309 1, igap, initPos_refit, mom_refit, padID, intersect_x, intersect_y, intersect_z );
1310
1311 // vector<Identifier> MuId =
1312 // road3D.ProjectToStrip(1,igap,initPos,mom_refit,padID,intersect_x,intersect_y,intersect_z);
1313
1314 vector<Identifier>::const_iterator mucid;
1315 int i = 0;
1316 for ( mucid = MuId.begin(); mucid != MuId.end(); mucid++, i++ )
1317 {
1318 MucRecHit* pHit = new MucRecHit( MucID::part( *mucid ), MucID::seg( *mucid ),
1319 MucID::gap( *mucid ), MucID::strip( *mucid ) );
1320 pHit->SetPadID( padID[i] );
1321 pHit->SetIntersectX( intersect_x[i] );
1322 pHit->SetIntersectY( intersect_y[i] );
1323 pHit->SetIntersectZ( intersect_z[i] );
1324
1325 m_pExpectedHits.push_back( pHit );
1326 }
1327 }
1328
1329 if ( m_endPart != -1 && m_ecLastLayer != -1 )
1330 {
1331 // for(int igap = 0; igap <= m_ecLastLayer; igap++){
1332 for ( int igap = 0; igap < 8; igap++ )
1333 {
1334 vector<int> padID;
1335 vector<float> intersect_x;
1336 vector<float> intersect_y;
1337 vector<float> intersect_z;
1338
1339 // refit---------------------------------
1340 float zx, zy, x0, y0;
1341 // int status = road3D.RefitNoCurrentGap(igap,zy,y0,zx,x0); //!!!!!different sequence
1342 int status = road3D.RefitNoCurrentGap( igap, zx, x0, zy, y0 );
1343 Hep3Vector mom_refit;
1344 if ( status == 0 )
1345 { // refit succeed
1346 float zDir = 1.0;
1347 if ( m_MucMomentum.z() != 0.0 ) zDir = m_MucMomentum.z() / fabs( m_MucMomentum.z() );
1348 float px_refit = zx * zDir;
1349 float py_refit = zy * zDir;
1350 float pz_refit = zDir;
1351 float segPhi = 0.25 * kPi * seg;
1352 // if vr is opposite to original, but both are very small, e.x. -0.01 -> + 0.01, or
1353 // you can say, pt~p, change it to guarantee vx, vy is correct.
1354 if ( part == 1 && px_refit * cos( segPhi ) + py_refit * sin( segPhi ) < 0.0 )
1355 { // when in barrel, pt dot segPhi < 0, angle between them > 90deg
1356 px_refit *= -1.0;
1357 py_refit *= -1.0;
1358 pz_refit *= -1.0;
1359 }
1360
1361 mom_refit.setX( px_refit );
1362 mom_refit.setY( py_refit );
1363 mom_refit.setZ( pz_refit );
1364 mom_refit.setMag( momentum );
1365 }
1366 else mom_refit = m_MucMomentum; // use the former momentum
1367
1368 HepPoint3D initPos_refit;
1369 if ( status == 0 )
1370 {
1371 float initPosx_refit, initPosy_refit, initPosz_refit;
1372 float sigmax_refit, sigmay_refit, sigmaz_refit;
1373 road3D.Project( 0, zx, x0, zy, y0, initPosx_refit, initPosy_refit, initPosz_refit );
1374 initPos_refit.setX( initPosx_refit - mom_refit.x() / momentum );
1375 initPos_refit.setY( initPosy_refit - mom_refit.y() / momentum );
1376 initPos_refit.setZ( initPosz_refit - mom_refit.z() / momentum );
1377 }
1378 else initPos_refit = initPos;
1379 // refit---------------------------------
1380 // cout<<"mom: "<<m_MucMomentum<<" "<<mom_refit<<" pos: "<<initPos<<"
1381 // "<<initPos_refit<<endl; vector<Identifier> MuId =
1382 // road3D.ProjectToStrip(m_endPart,igap,initPos,m_MucMomentum,padID,intersect_x,intersect_y,intersect_z);
1383
1384 vector<Identifier> MuId =
1385 road3D.ProjectToStrip( m_endPart, igap, initPos_refit, mom_refit, padID,
1386 intersect_x, intersect_y, intersect_z );
1387
1388 vector<Identifier>::const_iterator mucid;
1389 int i = 0;
1390 for ( mucid = MuId.begin(); mucid != MuId.end(); mucid++, i++ )
1391 {
1392 MucRecHit* pHit = new MucRecHit( MucID::part( *mucid ), MucID::seg( *mucid ),
1393 MucID::gap( *mucid ), MucID::strip( *mucid ) );
1394 pHit->SetPadID( padID[i] );
1395 pHit->SetIntersectX( intersect_x[i] );
1396 pHit->SetIntersectY( intersect_y[i] );
1397 pHit->SetIntersectZ( intersect_z[i] );
1398
1399 m_pExpectedHits.push_back( pHit );
1400 }
1401 }
1402 }
1403
1404 // cout<<"in RecMucTrack push back expected hits size= "<<m_pExpectedHits.size()<<" ?
1405 // "<<m_pHits.size()<<endl;
1406 for ( int i = 0; i < m_pExpectedHits.size(); i++ )
1407 {
1408 MucRecHit* ihit = m_pExpectedHits[i];
1409 // cout<<"expected Hits: "<<ihit->Part()<<" "<<ihit->Seg()<<" "<<ihit->Gap()<<"
1410 // "<<ihit->Strip()<<endl; cout<<"pad: "<<ihit->GetPadID()<<" pos:
1411 // "<<ihit->GetIntersectX()<<" "<<ihit->GetIntersectY()<<"
1412 // "<<ihit->GetIntersectZ()<<" "<<endl;
1413 }
1414
1415 for ( int i = 0; i < m_pHits.size(); i++ )
1416 {
1417 MucRecHit* ihit = m_pHits[i];
1418 // cout<<"attachedd Hits: "<<ihit->Part()<<" "<<ihit->Seg()<<" "<<ihit->Gap()<<"
1419 // "<<ihit->Strip()<<" isseed:
1420 // "<<ihit->HitIsSeed()<<endl;
1421 }
1422 // cout<<"1st: "<<brFirstLayer()<<" "<<ecFirstLayer()<<" last: "<<brLastLayer()<<"
1423 // "<<ecLastLayer()<<endl;
1424 }
1425 else
1426 {
1427 // calc distance between each hit with this track.
1428 // initialize distHits
1429 for ( int ihit = 0; ihit < m_pHits.size(); ihit++ )
1430 {
1431 m_distHits.push_back( -99 );
1432 m_distHits_quad.push_back( -99 );
1433 }
1434 }
1435 // cout << "Muc Pos: x " << m_MucPos.x() << " y " << m_MucPos.y() << " z " << m_MucPos.z() <<
1436 // endl;
1437}
1438
1439void RecMucTrack::ComputeNGapsWithHits() {
1440 int ngap = MucID::getGapMax();
1441 int gap, count = 0;
1442 vector<int> firedGap;
1443 for ( gap = 0; gap < ngap; gap++ ) firedGap.push_back( 0 );
1444
1445 vector<MucRecHit*>::const_iterator iHit;
1446 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ )
1447 {
1448 if ( *iHit )
1449 { // Check for a null pointer.
1450 gap = ( *iHit )->Gap();
1451 firedGap[gap] = 1;
1452 }
1453 }
1454
1455 for ( gap = 0; gap < ngap; gap++ ) { count += firedGap[gap]; }
1456
1457 m_numHits = m_pHits.size();
1458 m_numLayers = count;
1459 // cout<<"nLayer:\t"<<m_numLayers<<endl;
1460}
1461
1462int RecMucTrack::GetTotalHits() const { return m_pHits.size(); }
1463
1464int RecMucTrack::GetHitInGap( const int part, const int gap ) const {
1465 if ( part < 0 || part > 2 )
1466 {
1467 cout << "RecMucTrack::GetHitInGap(), invalid part " << part << endl;
1468 return 0;
1469 }
1470 if ( ( gap < 0 ) || ( gap >= (int)MucID::getGapNum( part ) ) )
1471 {
1472 cout << "RecMucTrack::GetHitInGap(), invalid gap " << gap << endl;
1473 return 0;
1474 }
1475
1476 vector<MucRecHit*>::const_iterator iHit;
1477 int hitsInGap = 0;
1478
1479 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ )
1480 {
1481
1482 if ( !( *iHit ) )
1483 {
1484 cout << "RecMucTrack::GetHitInGap() null hit pointer !" << endl;
1485 return 0;
1486 }
1487 else
1488 {
1489 if ( part == ( *iHit )->Part() && gap == ( *iHit )->Gap() ) hitsInGap++;
1490 }
1491 }
1492
1493 return hitsInGap;
1494}
1495
1496int RecMucTrack::GetHitInSeg( const int part, const int seg ) const {
1497 int num = GetHitInSegOrient( part, seg, 0 ) + GetHitInSegOrient( part, seg, 1 );
1498 return num;
1499}
1500
1501int RecMucTrack::GetHitInSegOrient( const int part, const int seg, const int orient ) const {
1502 if ( part < 0 || part > 2 )
1503 {
1504 cout << "RecMucTrack::GetHitInSeg(), invalid part " << part << endl;
1505 return 0;
1506 }
1507 if ( ( seg < 0 ) || ( seg >= (int)MucID::getSegNum( part ) ) )
1508 {
1509 cout << "RecMucTrack::GetHitInSeg(), invalid seg = " << seg << endl;
1510 return 0;
1511 }
1512
1513 vector<MucRecHit*>::const_iterator iHit;
1514 int hitsInSeg = 0;
1515
1516 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ )
1517 {
1518
1519 if ( !( *iHit ) )
1520 {
1521 cout << "RecMucTrack::GetHitInSeg(), null hit pointer !" << endl;
1522 return 0;
1523 }
1524 else
1525 {
1526 if ( part == ( *iHit )->Part() && seg == ( *iHit )->Seg() &&
1527 orient == ( *iHit )->GetGap()->Orient() )
1528 hitsInSeg++;
1529 }
1530 }
1531
1532 return hitsInSeg;
1533}
1534
1535int RecMucTrack::FindSegWithMaxHits( int& findPart, int& findSeg ) {
1536 int maxHits = 0;
1537 int hits = 0;
1538 for ( int part = 0; part < (int)MucID::getPartNum(); part++ )
1539 {
1540 for ( int seg = 0; seg < (int)MucID::getSegNum( part ); seg++ )
1541 {
1542 hits = GetHitInSeg( part, seg );
1543 if ( hits > maxHits )
1544 {
1545 maxHits = hits;
1546 findPart = part;
1547 findSeg = seg;
1548 }
1549 }
1550 }
1551
1552 return maxHits;
1553}
1554
1555void RecMucTrack::ComputeMaxHitsInGap() {
1556 int maxHits = 0;
1557 int hits = 0;
1558 for ( int part = 0; part < (int)MucID::getPartNum(); part++ )
1559 {
1560 for ( int gap = 0; gap < (int)MucID::getGapNum( part ); gap++ )
1561 {
1562 hits = GetHitInGap( part, gap );
1563 if ( hits > maxHits ) maxHits = hits;
1564 }
1565 }
1566
1567 m_maxHitsInLayer = maxHits;
1568}
1569
1570bool RecMucTrack::HasHit( const int part, const int seg, const int gap,
1571 const int strip ) const {
1572 bool found = false;
1573 vector<MucRecHit*>::const_iterator iHit;
1574
1575 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ )
1576 {
1577 if ( *iHit )
1578 { // Check for a null pointer.
1579 if ( ( *iHit )->Part() == part && ( *iHit )->Seg() == seg && ( *iHit )->Gap() == gap &&
1580 ( *iHit )->Strip() == strip )
1581 { found = true; }
1582 }
1583 }
1584
1585 return found;
1586}
1587
1588bool RecMucTrack::HasHitInGap( const int part, const int gap ) const {
1589 if ( ( gap < 0 ) || ( gap >= (int)MucID::getGapNum( part ) ) )
1590 {
1591 cout << "RecMucTrack::HasHitInGap-E2 invalid gap = " << gap << endl;
1592 return false;
1593 }
1594
1595 bool found = false;
1596 vector<MucRecHit*>::const_iterator iHit;
1597
1598 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ )
1599 {
1600 if ( *iHit )
1601 { // Check for a null pointer.
1602 if ( ( *iHit )->Part() == part && ( *iHit )->Gap() == gap ) { found = true; }
1603 }
1604 }
1605
1606 return found;
1607}
1608
1609int RecMucTrack::GetNSharedHits( const RecMucTrack* track2 ) const {
1610 if ( !track2 ) { return 0; }
1611
1612 int count = 0;
1613 vector<MucRecHit*>::const_iterator iHit1;
1614 vector<MucRecHit*>::const_iterator iHit2;
1615 MucRecHit * hit1, *hit2;
1616
1617 for ( iHit1 = m_pHits.begin(); iHit1 != m_pHits.end(); iHit1++ )
1618 {
1619 for ( iHit2 = track2->m_pHits.begin(); iHit2 != track2->m_pHits.end(); iHit2++ )
1620 {
1621 hit1 = ( *iHit1 );
1622 hit2 = ( *iHit2 );
1623
1624 if ( ( hit1 != 0 ) && ( hit2 != 0 ) )
1625 {
1626 if ( hit1->GetID() == hit2->GetID() ) { count++; }
1627 }
1628 }
1629 }
1630
1631 return count;
1632}
1633
1634MucRecHit* RecMucTrack::GetHit( const int part, const int gap ) const {
1635 if ( ( gap < 0 ) || ( gap >= (int)MucID::getGapNum( part ) ) )
1636 {
1637 cout << "RecMucTrack::Hit-E1 invalid gap = " << gap << endl;
1638 return 0;
1639 }
1640
1641 vector<MucRecHit*>::const_iterator iHit;
1642
1643 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ )
1644 {
1645 if ( *iHit )
1646 { // Check for a null pointer.
1647 if ( ( *iHit )->Part() == part && ( *iHit )->Gap() == gap ) { return ( *iHit ); }
1648 }
1649 }
1650
1651 return 0L;
1652}
1653
1654vector<MucRecHit*> RecMucTrack::GetHits() const { return m_pHits; }
1655
1656vector<MucRecHit*> RecMucTrack::GetExpectedHits() const { return m_pExpectedHits; }
1657
1658vector<long> RecMucTrack::GetHitIndices() const {
1659 vector<long> v;
1660
1661 vector<MucRecHit*>::const_iterator iHit;
1662 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ )
1663 {
1664 if ( *iHit )
1665 { // Check for a null pointer.
1666 long index = ( *iHit )->GetID();
1667 v.push_back( index );
1668 /*
1669 cout << " Muc2DRoad::HitIndices gap orientation twopack= "
1670 << (*iHit)->ChannelID().Plane() << " "
1671 << (*iHit)->ChannelID().Orient() << " "
1672 << (*iHit)->ChannelID().TwoPack() << endl;
1673 */
1674 }
1675 }
1676 return v;
1677}
1678
1679void RecMucTrack::ComputeTrackInfo( int fittingmethod ) {
1680 setVecHits( m_pHits );
1681 setExpHits( m_pExpectedHits );
1682 ComputeLastGap();
1683 ComputeNGapsWithHits();
1684 ComputeMaxHitsInGap();
1685 // ComputeDepth();
1686 ComputeDepth( fittingmethod );
1687 ComputeDistanceMatch();
1688 OutputUnitChange();
1689}
1690
1691void RecMucTrack::PrintHitsInfo() const {
1692 vector<MucRecHit*>::const_iterator iHit;
1693 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ )
1694 {
1695 if ( *iHit )
1696 { // Check for a null pointer.
1697 float xl, yl, zl;
1698 ( *iHit )->GetStrip()->GetCenterPos( xl, yl, zl );
1699 HepPoint3D vl( xl, yl, zl );
1700 HepPoint3D vg = ( *iHit )->GetGap()->TransformToGlobal( vl );
1701
1702 cout << " part " << ( *iHit )->Part() << " seg " << ( *iHit )->Seg() << " gap "
1703 << ( *iHit )->Gap() << " strip " << ( *iHit )->Strip() << " pos (" << vg.x() << ", "
1704 << vg.y() << ", " << vg.z() << ")" << endl;
1705 }
1706 }
1707}
1708
1709void RecMucTrack::OutputUnitChange() {
1710 if ( m_changeUnit == false )
1711 {
1712 m_depth /= 10;
1713
1714 m_xPos /= 10;
1715 m_yPos /= 10;
1716 m_zPos /= 10;
1717
1718 m_xPosSigma /= 10;
1719 m_yPosSigma /= 10;
1720 m_zPosSigma /= 10;
1721
1722 m_px /= 1000;
1723 m_py /= 1000;
1724 m_pz /= 1000;
1725
1726 m_distance /= 10;
1727 }
1728
1729 m_changeUnit = true;
1730}
momentum
**********INTEGER nmxhep !maximum number of particles DOUBLE PRECISION vhep INTEGER jdahep COMMON hepevt $ !serial number $ !number of particles $ !status code $ !particle ident KF $ !parent particles $ !childreen particles $ !four momentum
Definition BesBdkRc/src/fortran/HepEvt.h:23
x
Double_t x[10]
Definition DataBase/tau_mode.c:49
count
DOUBLE_PRECISION count[3]
Definition EvtPhokharaDef.hh:63
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:83
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
kSuperConductorR
const double kSuperConductorR
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucGeomSvc/MucConstant.h:11
kMinor
const double kMinor
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucGeomSvc/MucConstant.h:8
kPi
const double kPi
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucGeomSvc/MucConstant.h:6
kSuperConductorZ
const double kSuperConductorZ
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucGeomSvc/MucConstant.h:12
kMagnetField
const double kMagnetField
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucGeomSvc/MucConstant.h:13
kvC
const double kvC
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/MucTrackParameter.h:6
kMdcExtIterMax
const int kMdcExtIterMax
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/MucTrackParameter.h:7
kInsctWeight
const double kInsctWeight
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/MucTrackParameter.h:9
kDeltaPhi
const double kDeltaPhi
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/MucTrackParameter.h:5
intersection
int intersection(const HepPoint3D &c1, double r1, const HepPoint3D &c2, double r2, double eps, HepPoint3D &x1, HepPoint3D &x2)
Circle utilities.
Definition TMDCUtil.cxx:93
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
MucGeoGap::TransformToGlobal
HepPoint3D TransformToGlobal(const HepPoint3D pPoint) const
Transform a point from gap coordinate to global coordinate.
Definition MucGeoGap.cxx:575
MucGeoGap::Part
int Part() const
Get part identifier (0,2 for cap, 1 for barrel).
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucGeomSvc/MucGeoGap.h:76
MucGeoGap::TransformToGap
HepPoint3D TransformToGap(const HepPoint3D gPoint) const
Transform a point from global coordinate to gap coordinate.
Definition MucGeoGap.cxx:580
MucGeoGeneral::FindIntersections
vector< HepPoint3D > FindIntersections(const int part, const int gap, const HepPoint3D gPoint, const Hep3Vector gDirection)
Definition MucGeoGeneral.cxx:685
MucGeoGeneral::GetGap
MucGeoGap * GetGap(const int part, const int seg, const int gap) const
Get a pointer to the gap identified by (part,seg,gap).
Definition MucGeoGeneral.cxx:435
MucGeoGeneral::Instance
static MucGeoGeneral * Instance()
Get a pointer to the single instance of MucGeoGeneral.
Definition MucGeoGeneral.cxx:428
MucID::part
static int part(const Identifier &id)
Definition MucID.cxx:43
MucID::getGapMax
static value_type getGapMax()
Definition MucID.cxx:165
MucID::getPartNum
static value_type getPartNum()
Definition MucID.cxx:131
MucID::gap
static int gap(const Identifier &id)
Definition MucID.cxx:63
MucID::seg
static int seg(const Identifier &id)
Definition MucID.cxx:53
MucID::getSegNum
static value_type getSegNum(int part)
Definition MucID.cxx:134
MucID::strip
static int strip(const Identifier &id)
Definition MucID.cxx:73
MucID::getGapNum
static value_type getGapNum(int part)
Definition MucID.cxx:141
MucRec2DRoad::GetIntercept
float GetIntercept() const
Intercept of trajectory.
Definition MucRec2DRoad.cxx:1038
MucRec2DRoad::AttachHit
void AttachHit(MucRecHit *hit)
Attach the given hit to this road.
Definition MucRec2DRoad.cxx:99
MucRec2DRoad::GetSlope
float GetSlope() const
Slope of trajectory.
Definition MucRec2DRoad.cxx:1035
MucRec3DRoad::RefitNoCurrentGap
int RefitNoCurrentGap(const int &gap, float &slopeZX, float &interceptZX, float &slopeZY, float &interceptZY)
refit the 3D road without the assigned gap
Definition MucRec3DRoad.cxx:260
MucRec3DRoad::GetSlopeZX
float GetSlopeZX() const
Get Z-X dimension slope.
Definition MucRec3DRoad.cxx:233
MucRec3DRoad::GetSlopeZY
float GetSlopeZY() const
Get Z-Y dimension slope.
Definition MucRec3DRoad.cxx:239
MucRec3DRoad::TransformPhiRToXY
void TransformPhiRToXY(const float &vk, const float &vr, const float &k0, const float &r0, float &vx, float &vy, float &x0, float &y0) const
Where does the trajectory of this road intersect the reference plane?
Definition MucRec3DRoad.cxx:653
MucRec3DRoad::ProjectNoCurrentGap
void ProjectNoCurrentGap(const int &gap, float &x, float &y, float &z, float &sigmaX, float &sigmaY, float &sigmaZ, float &quad_x1, float &quad_y1, float &quad_z1, float &quad_x2, float &quad_y2, float &quad_z2)
Definition MucRec3DRoad.cxx:419
MucRec3DRoad::GetPart
int GetPart() const
In which part was this road found?
Definition MucRec3DRoad.cxx:144
MucRec3DRoad::ProjectToStrip
vector< Identifier > ProjectToStrip(const int &part, const int &gap, const HepPoint3D &gPoint, const Hep3Vector &gDirection)
Definition MucRec3DRoad.cxx:309
MucRec3DRoad::SetSimpleFitParams
void SetSimpleFitParams(const float &m_SlopeZX, const float &m_InterceptZX, const float &m_SlopeZY, const float &m_InterceptZY)
set the fit parameters : slope and intercept for XZ and YZ.
Definition MucRec3DRoad.cxx:129
MucRec3DRoad::GetDegreesOfFreedom
int GetDegreesOfFreedom() const
How many degrees of freedom in the trajectory fit?
Definition MucRec3DRoad.cxx:220
MucRec3DRoad::Project
void Project(const int &gap, float &x1, float &y1, float &z1, float &x2, float &y2, float &z2)
Where does the trajectory of this road intersect two surfaces of a specific gap?
Definition MucRec3DRoad.cxx:532
MucRec3DRoad::ProjectWithSigma
void ProjectWithSigma(const int &gap, float &x, float &y, float &z, float &sigmaX, float &sigmaY, float &sigmaZ)
Where does the trajectory of this road intersect a specific gap?
Definition MucRec3DRoad.cxx:337
MucRec3DRoad::GetReducedChiSquare
float GetReducedChiSquare() const
Chi-square parameter (per degree of freedom) of the trajectory fit.
Definition MucRec3DRoad.cxx:223
MucRecHit::GetGap
MucGeoGap * GetGap() const
Get geometry data for the gap containing this hit.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/MucRecHit.h:82
MucRecHit::GetCenterPos
Hep3Vector GetCenterPos() const
Get Center position of the strip (global coords).
Definition MucRecHit.cxx:93
MucRecHit::GetID
Identifier GetID() const
Get soft identifier of this hit.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/MucRecHit.h:67
RecMucTrack::GetTotalHits
int GetTotalHits() const
How many hits in all does this track contain?
Definition RecMucTrack.cxx:1462
RecMucTrack::LineFit
void LineFit(int fittingMethod)
Line fit with hits on a seg with max hits.
Definition RecMucTrack.cxx:1007
RecMucTrack::CorrectPos
void CorrectPos()
Correct current position of this trajectory.
Definition RecMucTrack.cxx:545
RecMucTrack::PrintHitsInfo
void PrintHitsInfo() const
Print Hits Infomation.
Definition RecMucTrack.cxx:1691
RecMucTrack::CorrectDir
void CorrectDir()
Correct direction of this trajectory.
Definition RecMucTrack.cxx:538
RecMucTrack::ComputeTrackInfo
void ComputeTrackInfo(int fittingmethod)
Get corresponding monte carlo track pointer.
Definition RecMucTrack.cxx:1679
RecMucTrack::GetNGapsWithHits
int GetNGapsWithHits() const
How many gaps provide hits attached to this track?
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/RecMucTrack.h:231
RecMucTrack::SetDefault
void SetDefault()
Definition RecMucTrack.cxx:187
RecMucTrack::SetMdcMomentum
void SetMdcMomentum(const float px, const float py, const float pz)
set start moment of the track in Mdc.
Definition RecMucTrack.cxx:209
RecMucTrack::setExpHits
void setExpHits(vector< MucRecHit * > &pHits)
RecMucTrack::SetMucPosSigma
void SetMucPosSigma(const float sigmax, const float sigmay, const float sigmaz)
Definition RecMucTrack.cxx:220
RecMucTrack::GetHit
MucRecHit * GetHit(const int part, const int gap) const
Get a pointer to the first hit attached in a particular gap.
Definition RecMucTrack.cxx:1634
RecMucTrack::IsInsideSuperConductor
bool IsInsideSuperConductor(Hep3Vector pos)
Definition RecMucTrack.cxx:318
RecMucTrack::GetHitInSeg
int GetHitInSeg(const int part, const int seg) const
How many hits does a segment contains.
Definition RecMucTrack.cxx:1496
RecMucTrack::~RecMucTrack
~RecMucTrack()
Destructor.
Definition RecMucTrack.cxx:194
RecMucTrack::AttachInsct
void AttachInsct(Hep3Vector insct)
Attach the intersection to this trajectory.
Definition RecMucTrack.cxx:534
RecMucTrack::brFirstLayer
int brFirstLayer() const
Which gap on Barrel is the first one with hits attached to this track?
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/RecMucTrack.h:186
RecMucTrack::SetExtMucPos
void SetExtMucPos(const float x, const float y, const float z)
Definition RecMucTrack.cxx:234
RecMucTrack::GetNSharedHits
int GetNSharedHits(const RecMucTrack *track) const
How many hits do two tracks share?
Definition RecMucTrack.cxx:1609
RecMucTrack::ComputeNGapsWithHits
void ComputeNGapsWithHits()
ComputeNGapsWithHits;.
Definition RecMucTrack.cxx:1439
RecMucTrack::SetMdcPos
void SetMdcPos(const float x, const float y, const float z)
set start position of the track in Mdc.
Definition RecMucTrack.cxx:205
RecMucTrack::SetCurrentDir
void SetCurrentDir(const float x, const float y, const float z)
set current direction of the trajectory.
Definition RecMucTrack.cxx:246
RecMucTrack::ComputeDepth
void ComputeDepth()
chi2 of line fit
Definition RecMucTrack.cxx:736
RecMucTrack::GetHitDistance
float GetHitDistance(const MucRecHit *hit)
Calculate the distance of the hit to the intersection in read direction.
Definition RecMucTrack.cxx:409
RecMucTrack::SetExtTrack
void SetExtTrack(RecExtTrack *extTrack)
set Ext track point.
Definition RecMucTrack.cxx:261
RecMucTrack::setTrackId
void setTrackId(const int trackId)
set the index for this track.
Definition RecMucTrack.cxx:199
RecMucTrack::HasHit
bool HasHit(const int part, const int seg, const int gap, const int strip) const
How many hits were attached in the gap with the most attached hits?
Definition RecMucTrack.cxx:1570
RecMucTrack::GetMdcExtTrack
void GetMdcExtTrack(Hep3Vector mdcStartPos, Hep3Vector mdcStartMomentum, int charge, Hep3Vector &mucStartPos, Hep3Vector &mucStartMomentum)
compute ext track myself from mdc.
Definition RecMucTrack.cxx:278
RecMucTrack::ComputeMaxHitsInGap
void ComputeMaxHitsInGap()
ComputeMaxHitsInGap;.
Definition RecMucTrack.cxx:1555
RecMucTrack::GetHitInGap
int GetHitInGap(const int part, const int gap) const
How many hits per gap does this track contain?
Definition RecMucTrack.cxx:1464
RecMucTrack::FindSegWithMaxHits
int FindSegWithMaxHits(int &part, int &seg)
Find the segment which contains most hits, return max hits number.
Definition RecMucTrack.cxx:1535
RecMucTrack::CalculateInsct
Hep3Vector CalculateInsct(const int part, const int seg, const int gap)
Calculate intersection from all hits attached on this gap.
Definition RecMucTrack.cxx:469
RecMucTrack::Project
void Project(const int &part, const int &gap, float &x, float &y, float &z, int &seg)
Where does the trajectory of this track intersect a specific gap?
Definition RecMucTrack.cxx:354
RecMucTrack::SetMucMomentum
void SetMucMomentum(const float px, const float py, const float pz)
set start moment of the track in Muc.
Definition RecMucTrack.cxx:227
RecMucTrack::setVecHits
void setVecHits(vector< MucRecHit * > &pHits)
reload setVecHits
RecMucTrack::GetHitIndices
vector< long > GetHitIndices() const
Get indices of all hits in the track.
Definition RecMucTrack.cxx:1658
RecMucTrack::GetHitInSegOrient
int GetHitInSegOrient(const int part, const int seg, const int orient) const
How many hits does a segment contains in one orient.
Definition RecMucTrack.cxx:1501
RecMucTrack::SetExtMucMomentum
void SetExtMucMomentum(const float px, const float py, const float pz)
set start moment of ext track in Muc.
Definition RecMucTrack.cxx:238
RecMucTrack::operator=
RecMucTrack & operator=(const RecMucTrack &orig)
Assignment constructor.
Definition RecMucTrack.cxx:84
RecMucTrack::GetHits
vector< MucRecHit * > GetHits() const
Get all hits on this track.
Definition RecMucTrack.cxx:1654
RecMucTrack::AttachHit
void AttachHit(MucRecHit *hit)
set corresponding monte carlo track pointer.
Definition RecMucTrack.cxx:325
RecMucTrack::SetCurrentInsct
void SetCurrentInsct(const float x, const float y, const float z)
set current intersection of the trajectory with strip plane.
Definition RecMucTrack.cxx:250
RecMucTrack::ComputeDistanceMatch
void ComputeDistanceMatch()
Compute distance match //2006.11.08.
Definition RecMucTrack.cxx:845
RecMucTrack::OutputUnitChange
void OutputUnitChange()
change unit
Definition RecMucTrack.cxx:1709
RecMucTrack::SetCurrentPos
void SetCurrentPos(const float x, const float y, const float z)
set current position of the trajectory.
Definition RecMucTrack.cxx:242
RecMucTrack::GetExpectedHits
vector< MucRecHit * > GetExpectedHits() const
Definition RecMucTrack.cxx:1656
RecMucTrack::SetMucPos
void SetMucPos(const float x, const float y, const float z)
set start position of the track in Muc. (after line fit and correction)
Definition RecMucTrack.cxx:213
RecMucTrack::GetHitDistance2
float GetHitDistance2(const MucRecHit *hit)
no abs value
Definition RecMucTrack.cxx:437
RecMucTrack::Extend
void Extend()
Extend mucpos and extmucpos to first layer of muc.
Definition RecMucTrack.cxx:958
RecMucTrack::HasHitInGap
bool HasHitInGap(const int part, const int gap) const
Does this track contain any hits in the given gap?
Definition RecMucTrack.cxx:1588
RecMucTrack::ComputeLastGap
void ComputeLastGap()
Comute last gap in barrel and endcap.
Definition RecMucTrack.cxx:547
RecMucTrack::RecMucTrack
RecMucTrack()
Constructor.
Definition RecMucTrack.cxx:27
RecMucTrack::AttachDirection
void AttachDirection(Hep3Vector dir)
Attach the direction to this trajectory.
Definition RecMucTrack.cxx:536