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