CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
HoughTrack.cxx
Go to the documentation of this file.
1#include "HoughTransAlg/HoughTrack.h"
2#include "TrkBase/TrkExchangePar.h"
3#include "TrkFitter/TrkCircleMaker.h"
4#include "TrkFitter/TrkHelixMaker.h"
5#include "TrkBase/TrkHitList.h"
6#include "MdcData/MdcRecoHitOnTrack.h"
7//#include "CgemData/CgemHitOnTrack.h"
8#include "HoughTransAlg/Hough.h"
9#include "Identifier/MdcID.h"
10TGraph* HoughTrack::m_cut1_cgem =NULL;
11TGraph* HoughTrack::m_cut2_cgem =NULL;
12TGraph* HoughTrack::m_cut1_ODC1 =NULL;
13TGraph* HoughTrack::m_cut2_ODC1 =NULL;
14TGraph* HoughTrack::m_cut1_ODC2 =NULL;
15TGraph* HoughTrack::m_cut2_ODC2 =NULL;
16TGraph* HoughTrack::m_cut[2][43] = {NULL};
17int HoughTrack::m_clearTrack=1;
18int HoughTrack::m_useCgemInGlobalFit=3;
19//int HoughTrack::m_maxLayer=0;
20HoughTrack::HoughTrack(int charge, double angle, double rho, double dAngle, double dRho, int trkID):Helix()
21{
22 bFieldZ(-bFieldZ());
23
24 m_trkID = trkID;
25 m_charge = (charge>0)?1:-1;
26 m_angle = angle;
27 m_rho = rho;
28 m_flag = 0;
29 m_dAngle = dAngle;
30 m_dRho = dRho;
31 m_dTanl = 0;
32 m_dDz = 0;
33 m_chi2 = 0;
34 m_circleFitStat = 0;
35
36 m_dr = 0;
37 m_phi0 = (rho>0)? angle:angle+M_PI;
38 m_phi0 = (m_charge<0)? m_phi0:m_phi0+M_PI;
39 if(m_phi0>2*M_PI)m_phi0-=2*M_PI;
40 if(m_phi0<0) m_phi0+=2*M_PI;
41 m_kappa = fabs(m_alpha*rho)*m_charge;
42 m_dz = 0.0;
43 m_tanl = 0.0;
44
45 HepPoint3D pivot(0,0,0);
46 HepVector a(5,0);
47 HepSymMatrix Ea(5,0);
48 a(2) = m_phi0;
49 a(3) = m_kappa;
50 set(pivot,a,Ea);
51
52 //Helix helix(pivot,a); helix.bFieldZ(-helix.bFieldZ());
53 //KalmanFit::Helix helix(pivot,a);
54 //m_helix = helix;
55 //m_hot.clear();
56 m_trkRecoTrk = NULL;
57
58 m_vecHitPnt.clear();
59 m_vecHitResidual.clear();
60 m_vecHitChi2.clear();
61
62 //m_dr = m_helix.dr();
63 //m_phi0 = m_helix.phi0();
64 //m_kappa = m_helix.kappa();
65 //m_dz = m_helix.dz();
66 //m_tanl = m_helix.tanl();
67
68 m_unused=0;
69 m_untried=0;
70 m_mcTrack = NULL;
71 m_cgemHitVector.clear();
72}
73
74HoughTrack::HoughTrack(int charge, const HepPoint3D & position, const Hep3Vector & momentum, int trkID):Helix(position,momentum,charge)
75{
76 bFieldZ(-bFieldZ());
77 //Helix helix(position, momentum, charge); helix.bFieldZ(-helix.bFieldZ());
78 //m_helix = helix;
79 m_trkID = trkID;
80 m_charge = charge;
81 //m_angle = m_helix.center().phi();
82 //m_rho = 1./m_helix.radius();
83 m_angle = center().phi();
84 m_rho = 1./radius();
85 m_flag = 0;
86 if(m_angle<0)m_angle=+2*M_PI;
87 m_circleFitStat = 0;
88
89 //m_dr = m_helix.dr();
90 //m_phi0 = m_helix.phi0();
91 //m_kappa = m_helix.kappa();
92 //m_dz = m_helix.dz();
93 //m_tanl = m_helix.tanl();
94 m_dr = dr();
95 m_phi0 = phi0();
96 m_kappa = kappa();
97 m_dz = dz();
98 m_tanl = tanl();
99
100 //m_hot.clear();
101 m_trkRecoTrk = NULL;
102
103 m_vecHitPnt.clear();
104 m_vecHitResidual.clear();
105 m_vecHitChi2.clear();
106
107 m_dAngle = 0;
108 m_dRho = 0;
109 m_dTanl = 0;
110 m_dDz = 0;
111 m_chi2 = 0;
112 m_mcTrack = NULL;
113 m_cgemHitVector.clear();
114}
115
116HoughTrack::HoughTrack(HepPoint3D& pivot, HepVector& a, int trkID):Helix(pivot,a)
117{
118 bFieldZ(-bFieldZ());
119 //Helix helix(pivot,a); helix.bFieldZ(-helix.bFieldZ());
120 //m_helix = helix;
121 m_trkID = trkID;
122 //m_charge = (m_helix.kappa()>0)? 1:-1;;
123 //m_angle = m_helix.center().phi();
124 //m_rho = 1./m_helix.radius();
125 m_charge = (kappa()>0)? 1:-1;;
126 m_angle = center().phi();
127 m_rho = 1./radius();
128 m_flag = 0;
129 m_circleFitStat = 0;
130 if(m_angle<0)m_angle=+2*M_PI;
131
132 //m_dr = m_helix.dr();
133 //m_phi0 = m_helix.phi0();
134 //m_kappa = m_helix.kappa();
135 //m_dz = m_helix.dz();
136 //m_tanl = m_helix.tanl();
137 m_dr = dr();
138 m_phi0 = phi0();
139 m_kappa = kappa();
140 m_dz = dz();
141 m_tanl = tanl();
142
143 //m_hot.clear();
144 m_trkRecoTrk = NULL;
145
146 m_vecHitPnt.clear();
147 m_vecHitResidual.clear();
148 m_vecHitChi2.clear();
149
150 m_dAngle = 0;
151 m_dRho = 0;
152 m_dTanl = 0;
153 m_dDz = 0;
154 m_chi2 = 0;
155 m_mcTrack = NULL;
156 m_cgemHitVector.clear();
157}
158
159
160/*
161 HoughTrack::HoughTrack(Helix& helix, int trkID):Helix(helix)
162 {
163 m_trkID = trkID;
164 m_helix = helix;
165 m_charge = (m_helix.kappa()>0)? 1:-1;;
166 m_angle = m_helix.center().phi();
167 m_rho = 1./m_helix.radius();
168 m_flag = 0;
169 if(m_angle<0)m_angle=+2*M_PI;
170
171 m_dr = m_helix.dr();
172 m_phi0 = m_helix.phi0();
173 m_kappa = m_helix.kappa();
174 m_dz = m_helix.dz();
175 m_tanl = m_helix.tanl();
176
177 m_hot.clear();
178 m_trkRecoTrk = NULL;
179
180 m_dAngle = 0;
181 m_dRho = 0;
182 m_dTanl = 0;
183 m_dDz = 0;
184 m_chi2 = 0;
185 }
186 */
187HoughTrack::HoughTrack()
188{
189 m_trkID = -1;
190 m_charge = 0;
191 m_angle = 0;
192 m_rho = 1e-6;
193 m_flag = 0;
194 HepVector a(5,0);
195 HepPoint3D pivot(0,0,0);
196 bFieldZ(-bFieldZ());
197 m_circleFitStat = 0;
198
199 m_dr = dr();
200 m_phi0 = phi0();
201 m_kappa = kappa();
202 m_dz = dz();
203 m_tanl = tanl();
204
205 //m_hot.clear();
206 m_trkRecoTrk = NULL;
207
208 m_vecHitPnt.clear();
209 m_vecHitResidual.clear();
210 m_vecHitChi2.clear();
211
212 m_dAngle = 0;
213 m_dRho = 0;
214 m_dTanl = 0;
215 m_dDz = 0;
216 m_chi2 = 0;
217 m_mcTrack = NULL;
218 m_cgemHitVector.clear();
219}
220
221///*
222HoughTrack::HoughTrack(const HoughTrack& other):
223 Helix(other),
224 m_trkID(other.m_trkID),
225 m_charge(other.m_charge),
226 m_angle(other.m_angle),
227 m_rho(other.m_rho),
228 m_flag(other.m_flag),
229
230 m_dAngle(other.m_dAngle),
231 m_dRho(other.m_dRho),
232 m_dTanl(other.m_dTanl),
233 m_dDz(other.m_dTanl),
234
235 m_dr(other.m_dr),
236 m_phi0(other.m_phi0),
237 m_kappa(other.m_kappa),
238 m_dz(other.m_dz),
239 m_tanl(other.m_tanl),
240
241 m_chi2(other.m_chi2),
242 m_trkRecoTrk(other.m_trkRecoTrk),
243
244 m_circleFitStat(other.m_circleFitStat),
245 m_vecHitPnt(other.m_vecHitPnt),
246 m_vecHitResidual(other.m_vecHitResidual),
247 m_vecHitChi2(other.m_vecHitChi2),
248 m_vecStereoHitPnt(other.m_vecStereoHitPnt),
249 m_vecStereoHitRes(other.m_vecStereoHitRes),
250 m_mcTrack(other.m_mcTrack),
251 m_cgemHitVector(other.m_cgemHitVector),
252 m_vecRecMdcHit(other.m_vecRecMdcHit)
253{}
254
255HoughTrack& HoughTrack::operator=(const HoughTrack& other)
256{
257 if (this == & other) return * this;
258
259 Helix::operator=(other);
260 m_trkID = other.m_trkID;
261 m_charge = other.m_charge;
262 m_angle = other.m_angle;
263 m_rho = other.m_rho;
264 m_flag = other.m_flag;
265
266 m_dAngle = other.m_dAngle;
267 m_dRho = other.m_dRho;
268 m_dTanl = other.m_dTanl;
269 m_dDz = other.m_dDz;
270
271 m_dr = other.m_dr;
272 m_phi0 = other.m_phi0;
273 m_kappa = other.m_kappa;
274 m_dz = other.m_dz;
275 m_tanl = other.m_tanl;
276
277 m_chi2 = other.m_chi2;
278 m_trkRecoTrk = other.m_trkRecoTrk;
279
280 m_circleFitStat=other.m_circleFitStat;
281 m_vecHitPnt = other.m_vecHitPnt;
282 m_vecHitResidual = other.m_vecHitResidual;
283 m_vecHitChi2 = other.m_vecHitChi2;
284
285 m_vecStereoHitPnt = other.m_vecStereoHitPnt;
286 m_vecStereoHitRes = other.m_vecStereoHitRes;
287 m_mcTrack = other.m_mcTrack;
288 m_cgemHitVector = other.m_cgemHitVector;
289 m_vecRecMdcHit = other.m_vecRecMdcHit;
290 return * this;
291}
292//*/
293/*
294 int HoughTrack::findHot(vector<HoughHit>& hitList, int flag)
295 {
296 int nHot(0);
297//m_chi2 = 0;
298for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
299HoughHit hit(*hitIt);
300//hit.print();
301if(judgeHalf(hit)<0)continue;
302if(hit.getFlag() == flag){
303double cut = 999.0;//FIXME
304double res = driftDistRes(hit);
305if(fabs(res)<cut){
306m_hot.push_back(hit);
307m_chi2 =+ res*res;
308nHot++;
309}
310}else{
311double cut = 999.0;//FIXME
312vector<HoughHit::S_Z> sz = hit.getSZ();
313vector<HoughHit::S_Z>::iterator iter = sz.begin();
314for(;iter!=sz.end();iter++){
315double s = iter->first;
316double z = iter->second;
317double res = z-(m_dz+s*m_tanl);
318if(fabs(res)<cut){
319m_hot.push_back(hit);
320break;
321nHot++;
322}
323}
324//if(calculateZ_S(hit)>0)m_hot.push_back(hit);
325}
326}
327return nHot;
328}
329//*/
330/*
331 int HoughTrack::findHot(vector<HoughHit>& hitList, vector<HoughHit>& hot, int flag)
332 {
333 int nHot(0);
334//m_chi2 = 0;
335for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
336HoughHit hit(*hitIt);
337//hit.print();
338if(judgeHalf(hit)<0)continue;
339if(hit.getFlag() == flag){
340double cut = 999.0;//FIXME
341double res = driftDistRes(hit);
342if(fabs(res)<cut){
343hot.push_back(hit);
344m_chi2 =+ res*res;
345nHot++;
346}
347}else{
348double cut = 999.0;//FIXME
349vector<HoughHit::S_Z> sz = hit.getSZ();
350vector<HoughHit::S_Z>::iterator iter = sz.begin();
351for(;iter!=sz.end();iter++){
352double s = iter->first;
353double z = iter->second;
354double res = z-(m_dz+s*m_tanl);
355if(fabs(res)<cut){
356hot.push_back(hit);
357break;
358nHot++;
359}
360}
361//if(calculateZ_S(hit)>0)m_hot.push_back(hit);
362}
363}
364return nHot;
365}
366//*/
367
368/*
369 int HoughTrack::findHot(vector<HoughHit>& hitList, vector<HoughHit*>& hot, int flag, int charge)
370 {
371
372 int nHot(0);
373 HoughHit* CgemXCluster_sel[3]={NULL,NULL,NULL};
374 double resid_CgemXCluster_sel[3]={99.,99.,99.};
375//m_chi2 = 0;
376for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
377//hit.print();
378if(judgeCharge(*hitIt)*charge<0) continue;// charge requirement
379if(hitIt->getFlag() == flag){
380//double cut = 999.0;//FIXME
381double cut1, cut2;
382int iLayer = hitIt->getLayer();
383int used = hitIt->getUse();
384XhitCutWindow(m_rho, iLayer, charge, cut1, cut2);
385double res = driftDistRes(*hitIt);
386//if(fabs(res)<cut)
387cout<<" win: "<<setw(5)<<cut1<<" ~ "<<setw(5)<<cut2;
388map<int,double>::iterator it_map;
389it_map=m_map_lay_d.find(iLayer);
390if(it_map==m_map_lay_d.end()||fabs(res)<fabs(it_map->second))
391{
392m_map_lay_d[iLayer]=res;
393m_map_lay_hit[iLayer]=&(*hitIt);
394}
395if(res>cut1&&res<cut2)
396{
397cout<<" selected!";
398//if(iLayer>3) {
399hot.push_back(&(*hitIt));
400m_vecHitPnt.push_back(&(*hitIt));
401m_vecHitResidual.push_back(res);
402m_hot.push_back((*hitIt));
403m_chi2 =+ res*res;
404m_setLayer.insert(iLayer);
405if(used==0) m_untried++;
406if(used==0||used==-1) m_unused++;
407if(judgeHalf((*hitIt))>0) {
408m_nHitFirstHalf++;
409if(used<=0) m_nHitUnused_FirstHalf++;
410}
411else {
412m_nHitSecondHalf++;
413if(used<=0) m_nHitUnused_SecondHalf++;
414}
415nHot++;
416//}
417//else {
418// only keep one CGEM cluster at each layer
419// if(fabs(resid_CgemXCluster_sel[iLayer])>fabs(res))
420// {
421// resid_CgemXCluster_sel[iLayer]=res;
422// CgemXCluster_sel[iLayer]=&(*hitIt);
423// }
424//}
425
426}// within window
427cout<<endl;
428}// only one flag
429//else{
430//double cut = 999.0;//FIXME
431//vector<HoughHit::S_Z> sz = hitIt->getSZ();
432//vector<HoughHit::S_Z>::iterator iter = sz.begin();
433//for(;iter!=sz.end();iter++){
434//double s = iter->first;
435//double z = iter->second;
436//double res = z-(m_dz+s*m_tanl);
437//if(fabs(res)<cut){
438//hot.push_back(&(*hitIt));
439//break;
440//nHot++;
441//}
442//}
443//if(calculateZ_S(hit)>0)m_hot.push_back(hit);
444//}
445//
446}// loop hits
447//////////
448////////// only keep one CGEM cluster at each layer
449////////for(int iLayer=0; iLayer<3; iLayer++)
450////////{
451//////// if(CgemXCluster_sel[iLayer]!=NULL) {
452//////// double res = resid_CgemXCluster_sel[iLayer];
453//////// HoughHit* hit = CgemXCluster_sel[iLayer];
454//////// int used = hit->getUse();
455//////// hot.push_back(hit);
456//////// m_vecHitPnt.push_back(hit);
457//////// m_vecHitResidual.push_back(res);
458//////// m_hot.push_back(*hit);
459//////// m_chi2 =+ res*res;
460//////// m_setLayer.insert(iLayer);
461//////// if(used==0) m_untried++;
462//////// if(used==0||used==-1) m_unused++;
463//////// if(judgeHalf(*hit)>0) {
464//////// m_nHitFirstHalf++;
465//////// if(used<=0) m_nHitUnused_FirstHalf++;
466//////// }
467//////// else {
468//////// m_nHitSecondHalf++;
469//////// if(used<=0) m_nHitUnused_SecondHalf++;
470//////// }
471//////// nHot++;
472//////// }
473////////}//
474m_maxGap=0;
475m_nGap=XGapSize(m_setLayer,m_maxGap);
476return nHot;
477}
478*/
479
480
481int HoughTrack::findXHot(vector<HoughHit*>& hitList, int charge)
482{
483
484 bool printFlag(false); //printFlag=true;
485 int nHot(0);
486 for(vector<HoughHit*>::iterator hitIt = hitList.begin(); hitIt != hitList.end(); hitIt++){
487 //cout<<judgeCharge(*(*hitIt))*charge<<" ";
488 //if(0 == flag){
489 if((*hitIt)->getFlag() != 0) continue;
490 if(judgeCharge(*hitIt)*charge<0) continue;// charge requirement
491 //if((*hitIt)->getPairHit()==NULL)continue;//FIXME
492 //if((*hitIt)->getPairHit()->getHalfCircle()>1)continue;//FIXME
493 if(printFlag)(*hitIt)->print();
494 //double cut = 999.0;//FIXME
495 //double cut1, cut2;
496 double cut[2] = {0};
497 int iLayer = (*hitIt)->getLayer();
498 int used = (*hitIt)->getUse();
499 XhitCutWindow(m_rho, iLayer, charge, cut[0], cut[1]);
500 double res = driftDistRes(*hitIt);
501 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
502 double res2 = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
503 //if(fabs(res)<cut)
504 //cout<<"("<<(*hitIt)->getLayer()<<", "<<(*hitIt)->getWire()<<") ";
505 if(printFlag)cout<<"res:"<<setw(12)<<res;
506 //cout<<setw(12)<<res2;
507 if(printFlag)cout<<" win: "<<setw(5)<<cut[0]<<" ~ "<<setw(5)<<cut[1];
508 //cout<<endl;
509 map<int,double>::iterator it_map;
510 it_map=m_map_lay_d.find(iLayer);
511 if(it_map==m_map_lay_d.end()||fabs(res)<fabs(it_map->second))
512 {
513 m_map_lay_d[iLayer]=res;
514 m_map_lay_hit[iLayer]=(*hitIt);
515 }
516 if(res>cut[0]&&res<cut[1])
517 {
518 if(printFlag)cout<<" selected!";
519 //if(iLayer>3)
520 m_vecHitPnt.push_back((*hitIt));
521 m_vecHitResidual.push_back(res);
522 //m_hot.push_back((*(*hitIt)));
523 //m_chi2 =+ res*res;
524 m_setLayer.insert(iLayer);
525 if(used==0) m_untried++;
526 if(used==0||used==-1) m_unused++;
527 if(judgeHalf(*hitIt)>0)
528 {
529 m_nHitFirstHalf++;
530 if(used<=0) m_nHitUnused_FirstHalf++;
531 }
532 else
533 {
534 m_nHitSecondHalf++;
535 if(used<=0) m_nHitUnused_SecondHalf++;
536 }
537 nHot++;
538 }// within window
539 if(printFlag)cout<<endl;
540 //}// X-view
541 /*
542 else {
543 if((*hitIt)->getFlag() == 0) continue;
544 //if((*hitIt)->getPairHit()==NULL)continue;//FIXME
545 //if((*hitIt)->getPairHit()->getHalfCircle()!=1)continue;//FIXME
546 int layer = (*hitIt)->getLayer();
547 int wire = (*hitIt)->getWire();
548 double Pt = fabs(pt());
549 double cut[2] = {0};
550 if(layer<3){
551 if(m_cut[0][layer+3]!=NULL)cut[0]=m_cut[0][layer+3]->Eval(Pt);
552 if(m_cut[1][layer+3]!=NULL)cut[1]=m_cut[1][layer+3]->Eval(Pt);
553 }else{
554 if(m_cut[0][layer]!=NULL)cut[0]=m_cut[0][layer]->Eval(Pt);
555 if(m_cut[1][layer]!=NULL)cut[1]=m_cut[1][layer]->Eval(Pt);
556 }
557 //cout<<"("<<setw(2)<<layer<<","<<setw(3)<<wire<<") ";
558 //cout<<"("<<setw(12)<<cut[0]<<" , "<<setw(13)<<cut[1]<<") ";
559 //cout<<Pt<<endl;
560 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
561 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
562 bool isNewHot(true);
563 for(vector<HoughHit*>::iterator hotIt = m_vecStereoHitPnt.begin(); hotIt != m_vecStereoHitPnt.end(); hotIt++){
564 if((*hitIt)->getHitID() == (*hotIt)->getHitID()){
565 isNewHot = false;
566 break;
567 }
568 }
569 if(!isNewHot){
570 //cout<<" OK!";
571 //cout<<endl;
572 continue;
573 }
574 if(cut[0]<res&&res<cut[1]){
575 //HepPoint3D szz(minS,zTrk,zHit);
576 //(*hitIt)->setHitPosition(szz);
577 m_vecStereoHitPnt.push_back((*hitIt));
578 //m_vecStereoHitRes.push_back(minRes);
579 //m_hot.push_back(*(*hitIt));
580 nHot++;
581 //cout<<" OK!";
582 //(*hitIt)->print();
583 }
584 }// V-view
585 // */
586 }// loop hits
587
588 m_maxGap=0;
589 m_nGap=XGapSize(m_setLayer,m_maxGap);
590 return nHot;
591}
592
593
594
595//for first half return 1, for second half return -1;
596int HoughTrack::judgeHalf(HoughHit* hit)
597{
598 double xHit = hit->getHitPosition().x();
599 double yHit = hit->getHitPosition().y();
600 //double xCenter = m_helix.center().x();
601 //double yCenter = m_helix.center().y();
602 double xCenter = center().x();
603 double yCenter = center().y();
604 if(xHit*yCenter > xCenter*yHit)return m_charge;
605 else if(xHit*yCenter < xCenter*yHit)return -m_charge;
606 else return 0;
607}
608
609int HoughTrack::judgeCharge(HoughHit* hit)
610{
611 double xHit = hit->getHitPosition().x();
612 double yHit = hit->getHitPosition().y();
613 //double xCenter = m_helix.center().x();
614 //double yCenter = m_helix.center().y();
615 double xCenter = center().x();
616 double yCenter = center().y();
617 double leftOrRight = xHit*yCenter - xCenter*yHit;
618 if(leftOrRight>0) return 1;
619 else return -1;
620}
621
622int HoughTrack::judgeCharge(double xHit, double yHit)
623{
624 //xHit = hit->getHitPosition().x();
625 //yHit = hit->getHitPosition().y();
626 //double xCenter = m_helix.center().x();
627 //double yCenter = m_helix.center().y();
628 double xCenter = center().x();
629 double yCenter = center().y();
630 double leftOrRight = xHit*yCenter - xCenter*yHit;
631 if(leftOrRight>0) return 1;
632 else return -1;
633}
634
635
636double HoughTrack::driftDistRes(HoughHit* hit)// for X-view hits only
637{
638 double residual(9999.);
639 HepPoint3D hitPoint = hit->getHitPosition();
640 if(hit->getHitType() == HoughHit::MDCHIT || hit->getHitType() == HoughHit::MDCMCHIT){
641 //HepPoint3D circleCenter = m_helix.center();
642 HepPoint3D circleCenter = center();
643 //double distance = circleCenter.perp(hitPoint);
644 double distance = (circleCenter-hitPoint).perp();
645 //hit->print();
646 //cout<<distance<<endl
647 //<<sqrt((circleCenter-hitPoint).x()*(circleCenter-hitPoint).x()+(circleCenter-hitPoint).y()*(circleCenter-hitPoint).y())<<endl
648 //<<sqrt((circleCenter.x()-hitPoint.x())*(circleCenter.x()-hitPoint.x())+(circleCenter.y()-hitPoint.y())*(circleCenter.y()-hitPoint.y()))<<endl;
649 //double Rc = m_helix.radius();
650 double Rc = fabs(radius());
651 double driftDist(0);
652 //if(hit->getHitType() == HoughHit::MDCMCHIT) driftDist = 0;
653 if(hit->getHitType() == HoughHit::MDCHIT) driftDist = hit->getDriftDist();
654 //residual = fabs(distance - Rc) - driftDist;
655 residual = driftDist-fabs(distance - Rc);
656 //cout<<"driftDist, distance, Rc = "<<driftDist<<", "<<distance<<", "<<Rc<<endl;
657 //double distance = (m_helix.center()).perp(hit->getHitPosition());
658 //res = fabs(distance - m_helix.radius()) - hit->getDriftDist();
659 }else if(hit->getHitType() == HoughHit::CGEMHIT || hit->getHitType() == HoughHit::CGEMMCHIT){
660 double rCgem = hitPoint.perp();
661 //double phiTrkFlt = m_helix.IntersectCylinder(rCgem);
662 double phiTrkFlt = IntersectCylinder(rCgem);
663 phiTrkFlt=judgeHalf(hit)*phiTrkFlt;
664 //HepPoint3D crossPoint = m_helix.x(phiTrkFlt);
665 HepPoint3D crossPoint = x(phiTrkFlt);
666 double phiCrossPoint = crossPoint.phi();
667 double phiMeasure = hitPoint.phi();
668 residual = phiMeasure - phiCrossPoint;
669 while(residual<-M_PI)residual += 2*M_PI;
670 while(residual> M_PI)residual -= 2*M_PI;
671 residual = rCgem*residual;
672 }
673 hit->setResidual(residual);
674 return residual;
675}
676
677//void HoughTrack::dropHot(HoughHit& hit)
678//{
679//for(HitVector_Iterator hitIt = m_hot.begin(); hitIt != m_hot.end();hitIt++){
680//if(hitIt->getHitID() == hit.getHitID())m_hot.erase(hitIt);
681//}
682//}
683
684//void HoughTrack::dropHot(int hitID)
685//{
686//for(HitVector_Iterator hitIt = m_hot.begin(); hitIt != m_hot.end();hitIt++){
687//if(hitIt->getHitID() == hitID)m_hot.erase(hitIt);
688//}
689//}
690
691/*
692 TrkErrCode HoughTrack::fitHelix(const MdcDetector* mdcDetector, TrkContextEv* trkContextEv, double bunchT0, vector<MdcHit*> mdcHitCol)
693 {
694 double d0 = -m_dr;
695 double phi0 = m_phi0+M_PI/2;
696 double omega = m_kappa/fabs(alpha());
697 double z0 = m_dz;
698 double tanl = m_tanl;
699 while(phi0>M_PI)phi0-=2*M_PI;
700 while(phi0<-M_PI)phi0+=2*M_PI;
701
702 TrkExchangePar helixTrk(d0,phi0,omega,z0,tanl);
703 TrkHelixMaker helixfactory;
704 double chisum =0.;
705 TrkRecoTrk* trkRecoTrk = helixfactory.makeTrack(helixTrk, chisum, *trkContextEv, bunchT0);
706 bool permitFlips = true;
707 bool lPickHits = true;
708 helixfactory.setFlipAndDrop(*trkRecoTrk, permitFlips, lPickHits);
709 TrkHitList* trkHitList = trkRecoTrk->hits();
710
711//--- convert hits
712double ddCut=1;
713for(HitVector_Iterator hitIt = m_hot.begin(); hitIt != m_hot.end();hitIt++){
714if(hitIt->getFlag()!=0)continue;
715if(hitIt->getHitType()==HoughHit::MDCHIT){
716//const MdcDigi* mdcDigi = hitIt->getDigi();
717//MdcHit *mdcHit = new MdcHit(mdcDigi,mdcDetector);
718MdcHit *mdcHit;
719vector<MdcHit*>::iterator imdcHit = mdcHitCol.begin();
720for(;imdcHit != mdcHitCol.end();imdcHit++){
721if((*imdcHit)->mdcId() == hitIt->getDigi()->identify()){
722mdcHit = *imdcHit;
723break;
724}
725}
726mdcHit->setCalibSvc(hitIt->getMdcCalibFunSvc());
727mdcHit->setCountPropTime(true);
728//if(mdcHit->driftDist(m_bunchT0,0)>ddCut)continue;//FIXME
729int newAmbig = 0;
730MdcRecoHitOnTrack mdcRecoHitOnTrack(*mdcHit, newAmbig, bunchT0*1.e9);
731MdcHitOnTrack* mdcHitOnTrack = &mdcRecoHitOnTrack;
732HepPoint3D midPoint = hitIt->getHitPosition();
733double fltLength = flightLength(midPoint);
734//double fltLength = flightLength(hitIt->getHitPosition());
735mdcHitOnTrack->setFltLen(fltLength);
736trkHitList->appendHot(mdcHitOnTrack);
737}else if(hitIt->getHitType()==HoughHit::CGEMHIT){
738CgemHitOnTrack* cgemHitOnTrack = new CgemHitOnTrack(*(hitIt->getCgemCluster()),bunchT0*1.e9);
739cgemHitOnTrack->setCgemGeomSvc(hitIt->getCgemGeomSvc());
740cgemHitOnTrack->setCgemCalibFunSvc(hitIt->getCgemCalibFunSvc());
741trkHitList->appendHot(cgemHitOnTrack);
742}
743}
744//---fitting
745TrkHitList* qhits = trkRecoTrk->hits();
746TrkErrCode err=qhits->fit();
747
748if(err.success()){
749const TrkFit* fitResult = trkRecoTrk->fitResult();
750TrkExchangePar par = fitResult->helix(0.);
751m_dr = -par.d0();
752m_phi0 = par.phi0()-M_PI/2;
753m_kappa = par.omega()*fabs(alpha());
754m_dz = par.z0();
755m_tanl = par.tanDip();
756while(m_phi0>2*M_PI)m_phi0-=2*M_PI;
757while(m_phi0<0)m_phi0+=2*M_PI;
758//HepVector hepVector(5,0);
759//hepVector(1) = m_dr;
760//hepVector(2) = m_phi0;
761//hepVector(3) = m_kappa;
762//hepVector(4) = m_dz;
763//hepVector(5) = m_tanl;
764//a(hepVector);
765updateHelix();
766
767m_trkRecoTrk = trkRecoTrk;
768m_chi2 = fitResult->chisq();
769}
770return err;
771}
772*/
773
774int HoughTrack::calculateZ_S(HoughHit* hit)
775{
776 int nTangency(0);
777 if(hit->getFlag()==0)return nTangency;
778 //if(hit->getPairHit()==NULL)return nTangency;
779 //if(hit->getPairHit()->getHalfCircle()!=1)return nTangency;
780 double s(0),z(0);
781 double xc = center().x();
782 double yc = center().y();
783 double rTrack = radius();// signed FIXME
784 if(hit->getHitType()==HoughHit::CGEMHIT){
785 if(hit->getFlag()!=2)return nTangency;
786 //for(HitVector_Iterator hitIt = m_hot.begin(); hitIt != m_hot.end();hitIt++)
787 for(vector<HoughHit*>::iterator hotIt = m_vecHitPnt.begin(); hotIt != m_vecHitPnt.end(); hotIt++){
788 if((*hotIt)->getFlag()!=0)continue;
789 if((*hotIt)->getHitType()==HoughHit::MDCHIT)continue;
790 if((*hotIt)->getCgemCluster()->getclusterid() == hit->getCgemCluster()->getclusterflagb()){
791 HepPoint3D point = hit->getHitPosition();
792 s = flightArc(point);
793 z = point.z();
794 // double rHit = point.perp();
795 // s = flightArc(rHit);
796 // z = getZfrom ... // FIXME
797 pair<double,double> sz(s,z);
798 hit->addSZ(sz);
799 nTangency++;
800 //double s = flightArc(hit->getHitPosition());
801 }
802 }
803 }else{
804 hit->resetSZ();
805 double drift = hit->getDriftDist();
806 const MdcGeoWire* wire = hit->getMdcGeomSvc()->Wire(hit->getLayer(),hit->getWire());
807 HepPoint3D westPoint = wire->Forward();
808 HepPoint3D eastPoint = wire->Backward();
809 double xEast = eastPoint.x()/10.0;
810 double xWest = westPoint.x()/10.0;
811 double yEast = eastPoint.y()/10.0;
812 double yWest = westPoint.y()/10.0;
813 double zEast = eastPoint.z()/10.0;
814 double zWest = westPoint.z()/10.0;
815 //cout<<"wast:x,y,z: "<<xWest<<", "<<yWest<<", "<<zWest<<endl;
816 //cout<<"east:x,y,z: "<<xEast<<", "<<yEast<<", "<<zEast<<endl;
817 //cout<<"Xc, Yc, Rc: "<<xc<<", "<<yc<<", "<<rTrack<<endl;
818 double west2east = sqrt((xEast-xWest)*(xEast-xWest)+(yEast-yWest)*(yEast-yWest));
819
820 double slope = (yEast-yWest)/(xEast-xWest);
821 double intercept = (yWest-slope*xWest+yEast-slope*xEast)/2;
822 double a = 1+slope*slope;
823 double b = -2*(xc+slope*yc-slope*intercept);
824 double c1 = xc*xc+(yc-intercept)*(yc-intercept)-(rTrack+drift)*(rTrack+drift);
825 double c2 = xc*xc+(yc-intercept)*(yc-intercept)-(rTrack-drift)*(rTrack-drift);
826 //double c1 = intercept*(2*yc-intercept)+drift*(drift+rTrack);
827 //double c2 = intercept*(2*yc-intercept)+drift*(drift-rTrack);
828 double delta1 = (b*b-4*a*c1);
829 double delta2 = (b*b-4*a*c2);
830 //cout<<"a,b,c1,c2: "<<a<<", "<<b<<", "<<c1<<", "<<c2<<endl;
831 //cout<<"delta: "<<delta1<<", "<<delta2<<endl;
832 if(delta1>=0){
833 double x1 = (-b+sqrt(delta1))/(2*a);
834 double x2 = (-b-sqrt(delta1))/(2*a);
835 double y1 = slope*x1+intercept;
836 double y2 = slope*x2+intercept;
837 if((x1-xWest)*(x1-xEast)<0){
838 HepPoint3D point(x1,y1,0);
839 s = flightArc(point);
840 //s = flightArc(hit->getHitPosition());
841 double l = sqrt((x1-xWest)*(x1-xWest)+(y1-yWest)*(y1-yWest));
842 z = zWest + l/west2east*fabs((zEast-zWest));
843 pair<double,double> sz(s,z);
844 hit->addSZ(sz);
845 nTangency++;
846 HepPoint3D position(x1,y1,z);
847 //hit->addPosition(position);
848 //cout<<"hit outside track, bigger root, x,y,z,s: "<<x1<<", "<<y1<<", "<<z<<", "<<s<<endl;
849 }
850 if((x2-xWest)*(x2-xEast)<0){
851 HepPoint3D point(x2,y2,0);
852 s = flightArc(point);
853 //s = flightArc(hit->getHitPosition());
854 double l = sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
855 z = zWest + l/west2east*fabs((zEast-zWest));
856 pair<double,double> sz(s,z);
857 hit->addSZ(sz);
858 nTangency++;
859 HepPoint3D position(x2,y2,z);
860 //cout<<"hit outside track, small root, x,y,z,s: "<<x2<<", "<<y2<<", "<<z<<", "<<s<<endl;
861 //hit->addPosition(position);
862 }
863 }
864
865 if( delta2>=0){
866 double x1 = (-b+sqrt(delta2))/(2*a);
867 double x2 = (-b-sqrt(delta2))/(2*a);
868 double y1 = slope*x1+intercept;
869 double y2 = slope*x2+intercept;
870 if((x1-xWest)*(x1-xEast)<0){
871 HepPoint3D point(x1,y1,0);
872 s = flightArc(point);
873 //s = flightArc(hit->getHitPosition());
874 double l = sqrt((x1-xWest)*(x1-xWest)+(y1-yWest)*(y1-yWest));
875 z = zWest + l/west2east*fabs((zEast-zWest));
876 pair<double,double> sz(s,z);
877 hit->addSZ(sz);
878 nTangency++;
879 HepPoint3D position(x1,y1,z);
880 //cout<<"hit outside track, bigger root, x,y,z,s: "<<x1<<", "<<y1<<", "<<z<<", "<<s<<endl;
881 //hit->addPosition(position);
882 }
883 if((x2-xWest)*(x2-xEast)<0){
884 HepPoint3D point(x2,y2,0);
885 s = flightArc(point);
886 //s = flightArc(hit->getHitPosition());
887 //double l = zWest + sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
888 //z = l/west2east*fabs((zEast-zWest));
889 double l = sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
890 z = zWest + l/west2east*fabs((zEast-zWest));
891 pair<double,double> sz(s,z);
892 hit->addSZ(sz);
893 nTangency++;
894 HepPoint3D position(x2,y2,z);
895 //cout<<"hit outside track, small root, x,y,z,s: "<<x2<<", "<<y2<<", "<<z<<", "<<s<<endl;
896 //hit->addPosition(position);
897 }
898 }
899 }
900 return nTangency;
901}
902
903bool innerLayer(HoughHit hit1, HoughHit hit2)
904{
905 return hit1.getLayer()<hit2.getLayer();
906}
907
908//void HoughTrack::sortHot()
909//{
910//std::sort(m_hot.begin(),m_hot.end(),innerLayer);
911//}
912
913bool sortByLayer(HoughHit* hit1, HoughHit* hit2)
914{
915 return hit1->getLayer()<hit2->getLayer();
916}
917
918void HoughTrack::sortHot(vector<HoughHit*>& hotList)
919{
920 std::sort(hotList.begin(),hotList.end(),sortByLayer);
921}
922
923void HoughTrack::updateHelix()
924{
925 pivot(HepPoint3D(0,0,0));
926 HepVector hepVector(5,0);
927 hepVector(1) = m_dr;
928 hepVector(2) = m_phi0;
929 hepVector(3) = m_kappa;
930 hepVector(4) = m_dz;
931 hepVector(5) = m_tanl;
932 a(hepVector);
933}
934
935void HoughTrack::update(double angle, double rho)
936{
937 m_angle = angle;
938 m_rho = rho;
939 m_phi0 = (rho>0)? angle:angle+M_PI;
940 m_phi0 = (m_charge<0)? m_phi0:m_phi0+M_PI;
941 if(m_phi0>2*M_PI)m_phi0-=2*M_PI;
942 if(m_phi0<0) m_phi0+=2*M_PI;
943 m_kappa = fabs(m_alpha*rho)*m_charge;
944 HepPoint3D pivot(0,0,0);
945 HepVector a(5,0);
946 HepSymMatrix Ea(5,0);
947 a(2) = m_phi0;
948 a(3) = m_kappa;
949 set(pivot,a,Ea);
950}
951
952
953void HoughTrack::updateCirclePar(double dr, double phi0, double kappa)
954{
955 m_dr = dr;
956 m_phi0=phi0;
957 m_kappa=kappa;
958 updateHelix();
959
960}
961
962void HoughTrack::print()
963{
964 cout
965 <<setw(12)<<"trkID:" <<setw(15)<<m_trkID
966 <<setw(12)<<"flag:" <<setw(15)<<m_flag
967 //<<setw(12)<<"q:" <<setw(15)<<m_charge
968 <<setw(12)<<"pt:" <<setw(15)<<pt()
969 <<setw(12)<<"angle:" <<setw(15)<<m_angle
970 <<setw(12)<<"rho:" <<setw(15)<<m_rho
971 <<endl
972 <<setw(12)<<"dr:" <<setw(15)<<dr()
973 <<setw(12)<<"phi0:" <<setw(15)<<phi0()
974 <<setw(12)<<"kappa:" <<setw(15)<<kappa()
975 <<setw(12)<<"dz:" <<setw(15)<<dz()
976 <<setw(12)<<"tanl:" <<setw(15)<<tanl()
977 <<setw(12)<<"chi2:" <<setw(15)<<getChi2()
978 <<endl
979 //<<setw(12)<<"dr:" <<setw(15)<<m_dr
980 //<<setw(12)<<"phi0:" <<setw(15)<<m_phi0
981 //<<setw(12)<<"kappa:" <<setw(15)<<m_kappa
982 //<<setw(12)<<"dz:" <<setw(15)<<m_dz
983 //<<setw(12)<<"tanl:" <<setw(15)<<m_tanl
984 //<<endl
985 //<<setw(12)<<"dAngle:" <<setw(15)<<m_dAngle
986 //<<setw(12)<<"dRho:" <<setw(15)<<m_dRho
987 //<<setw(12)<<"dTanl:" <<setw(15)<<m_dTanl
988 //<<setw(12)<<"dDz:" <<setw(15)<<m_dDz
989 //<<setw(12)<<"chi2:" <<setw(15)<<m_chi2
990 //<<setw(12)<<":"<<setw(15)<<
991 //<<setw(12)<<":"<<setw(15)<<
992 //<<setw(12)<<":"<<setw(15)<<
993 ;
994 int nHot = getHotList().size();
995 int nAxialHot = 0;
996 int nStereoHot = 0;
997 int nXCluster = 0;
998 int nXVCluster = 0;
999 vector<HoughHit*> hotList;
1000 vector<HoughHit*> axialHot = getVecHitPnt();
1001 vector<HoughHit*> stereoHot = getVecStereoHitPnt();
1002 for(vector<HoughHit*>::iterator hotIt = axialHot.begin(); hotIt != axialHot.end(); hotIt++){
1003 if((**hotIt).getHitType()==HoughHit::CGEMHIT)nXCluster++;
1004 if((**hotIt).getHitType()==HoughHit::MDCHIT)nAxialHot++;
1005 }
1006 for(vector<HoughHit*>::iterator hotIt = stereoHot.begin(); hotIt != stereoHot.end(); hotIt++){
1007 if((**hotIt).getHitType()==HoughHit::CGEMHIT)nXVCluster++;
1008 if((**hotIt).getHitType()==HoughHit::MDCHIT)nStereoHot++;
1009 }
1010 cout
1011 <<setw(12)<<"nHot:" <<setw(15)<<nHot
1012 <<setw(12)<<"nAxialHot0:" <<setw(15)<<nAxialHot
1013 <<setw(12)<<"nStereoHot0:" <<setw(15)<<nStereoHot
1014 <<setw(12)<<"nXCluster0:" <<setw(15)<<nXCluster
1015 <<setw(12)<<"nXVCluster0:" <<setw(15)<<nXVCluster
1016 <<endl;
1017 //cout<<endl;
1018}
1019
1020void HoughTrack::printHot()
1021{
1022 vector<HoughHit*> hitList = getHotList();
1023 //vector<HoughHit*> hitList = getVecStereoHitPnt();
1024 for(vector<HoughHit*>::iterator hotIt = hitList.begin(); hotIt != hitList.end();hotIt++){
1025 HoughHit* hitIt = *hotIt;
1026 hitIt->print();
1027 //if(hitIt->getLayer()>3)continue;
1028 //cout<<hitIt->getHitID();
1029 //cout<<" ("<<hitIt->getLayer()<<", "<<hitIt->getWire()<<") "<<hitIt->getFlag();
1030 //if(hitIt->getHitType()==HoughHit::CGEMMCHIT||hitIt->getHitType()==HoughHit::MDCMCHIT)cout<<" halfCircle:"<<hitIt->getHalfCircle();
1031 //else if(hitIt->getPairHit()!=NULL)cout<<" halfCircle:"<<hitIt->getPairHit()->getHalfCircle();
1032 //HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1033 //double res = hitIt->residual(this, positionOntrack, positionOnDetector);
1034 //cout<<" res:"<<setw(10)<<res<<" ";
1035 //cout<<endl;
1036 }
1037 cout<<endl;
1038 /*
1039 for(vector<HoughHit*>::iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
1040 cout<<setw(4)<<(*hitIt)->getHitID();
1041 cout<<"("<<setw(2)<<(*hitIt)->getLayer()<<","<<setw(3)<<(*hitIt)->getWire()<<") ";
1042 cout<<"("<<setw( 9)<<(*hitIt)->getHitPosition().x()<<","<<setw( 9)<<(*hitIt)->getHitPosition().y()<<","<<setw( 9)<<(*hitIt)->getHitPosition().z()<<") ";
1043 cout<<"flag:"<<setw(3)<<(*hitIt)->getFlag();
1044 cout<<"use:"<<setw(3)<<(*hitIt)->getUse();
1045 if((*hitIt)->getHitType()==HoughHit::CGEMMCHIT||(*hitIt)->getHitType()==HoughHit::MDCMCHIT){
1046 vector<int> trkID = (*hitIt)->getTrkID();
1047 cout<<"TrkID:"<<setw(4)<<trkID[0];
1048 }
1049 if((*hitIt)->getHitType()==HoughHit::CGEMHIT){
1050 if((*hitIt)->getPairHit()!=NULL)cout<<"TrkID:"<<setw(3)<<((*hitIt)->getPairHit()->getTrkID())[0];
1051 cout<<"res:"<<setw(10)<<(*hitIt)->residual(this)<<" ";
1052 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1053 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
1054 cout<<"res:"<<setw(10)<<res<<" ";
1055 //cout<<"["<<setw(4)<<(*hitIt)->getCgemCluster()->getclusterflagb()<<", "<<setw(4)<<(*hitIt)->getCgemCluster()->getclusterflage()<<"] ";
1056 }
1057 if((*hitIt)->getHitType()==HoughHit::MDCHIT){
1058 cout<<"TrkID:"<<setw(3)<<(*hitIt)->getDigi()->getTrackIndex();
1059 cout<<"res:"<<setw(10)<<(*hitIt)->residual(this)<<" ";
1060 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1061 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
1062 cout<<"res:"<<setw(10)<<res<<" ";
1063 //cout<<(*hitIt)->getDigi()->identify()<<" ";
1064 //cout<<"("<<(*hitIt)->getPairHit()->getLayer()<<", "<<(*hitIt)->getPairHit()->getWire()<<") ";
1065 //if((*hitIt)->getPairHit()!=NULL)cout<<(*hitIt)->getDriftDist()<<" "<<(*hitIt)->getPairHit()->getDriftDist()<<" "<<(*hitIt)->getDriftDist()-(*hitIt)->getPairHit()->getDriftDist()<<endl;
1066 }
1067 if((*hitIt)->getPairHit()!=NULL){
1068 cout<<setw(4)<<(*hitIt)->getPairHit()->getHitID();
1069 //if((*hitIt)->driftTime()>1500.)cout<<(*hitIt)->getDriftDist()<<" "<<(*hitIt)->getPairHit()->getDriftDist();
1070 }
1071 //HepPoint3D hitPoint = (*hitIt)->getHitPosition();
1072 //cout<<flightLength(hitPoint);
1073 //else cout<<"NULL";
1074 cout<<endl;
1075 }
1076
1077 //vector<HoughHit*> hitList = getVecHitPnt();
1078 hitList = getVecStereoHitPnt();
1079 for(vector<HoughHit*>::iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
1080 cout<<setw(4)<<(*hitIt)->getHitID();
1081 cout<<"("<<setw(2)<<(*hitIt)->getLayer()<<","<<setw(3)<<(*hitIt)->getWire()<<") ";
1082 cout<<"("<<setw( 9)<<(*hitIt)->getHitPosition().x()<<","<<setw( 9)<<(*hitIt)->getHitPosition().y()<<","<<setw( 9)<<(*hitIt)->getHitPosition().z()<<") ";
1083 cout<<"flag:"<<setw(3)<<(*hitIt)->getFlag();
1084 cout<<"use:"<<setw(3)<<(*hitIt)->getUse();
1085 if((*hitIt)->getHitType()==HoughHit::CGEMMCHIT||(*hitIt)->getHitType()==HoughHit::MDCMCHIT){
1086 vector<int> trkID = (*hitIt)->getTrkID();
1087 cout<<"TrkID:"<<setw(4)<<trkID[0];
1088 }
1089 if((*hitIt)->getHitType()==HoughHit::CGEMHIT){
1090 if((*hitIt)->getPairHit()!=NULL)cout<<"TrkID:"<<setw(3)<<((*hitIt)->getPairHit()->getTrkID())[0];
1091 cout<<"res:"<<setw(10)<<(*hitIt)->residual(this)<<" ";
1092 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1093 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
1094 cout<<"res:"<<setw(10)<<res<<" ";
1095 //cout<<"["<<setw(4)<<(*hitIt)->getCgemCluster()->getclusterflagb()<<", "<<setw(4)<<(*hitIt)->getCgemCluster()->getclusterflage()<<"] ";
1096 }
1097 if((*hitIt)->getHitType()==HoughHit::MDCHIT){
1098 cout<<"TrkID:"<<setw(3)<<(*hitIt)->getDigi()->getTrackIndex();
1099 cout<<"res:"<<setw(10)<<(*hitIt)->residual(this)<<" ";
1100 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1101 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
1102 cout<<"res:"<<setw(10)<<res<<" ";
1103 //cout<<(*hitIt)->getDigi()->identify()<<" ";
1104 //cout<<"("<<(*hitIt)->getPairHit()->getLayer()<<", "<<(*hitIt)->getPairHit()->getWire()<<") ";
1105 //if((*hitIt)->getPairHit()!=NULL)cout<<(*hitIt)->getDriftDist()<<" "<<(*hitIt)->getPairHit()->getDriftDist()<<" "<<(*hitIt)->getDriftDist()-(*hitIt)->getPairHit()->getDriftDist()<<endl;
1106 }
1107 if((*hitIt)->getPairHit()!=NULL){
1108 cout<<setw(4)<<(*hitIt)->getPairHit()->getHitID();
1109 //if((*hitIt)->driftTime()>1500.)cout<<(*hitIt)->getDriftDist()<<" "<<(*hitIt)->getPairHit()->getDriftDist();
1110 }
1111//HepPoint3D hitPoint = (*hitIt)->getHitPosition();
1112//cout<<flightLength(hitPoint);
1113//else cout<<"NULL";
1114cout<<endl;
1115}
1116//*/
1117}
1118
1119void HoughTrack::printRecMdcHitVec()
1120{
1121 cout<<"HoughTrack::printRecMdcHitVec(): "<<endl;
1122 vector<RecMdcHit>::iterator iter_recMdcHit = m_vecRecMdcHit.begin();
1123 for(; iter_recMdcHit!=m_vecRecMdcHit.end(); iter_recMdcHit++)
1124 {
1125 Identifier mdcid = iter_recMdcHit->getMdcId();
1126 int layer = MdcID::layer(mdcid);
1127 int wire = MdcID::wire(mdcid);
1128 cout<<" hit at layer "<<layer<<" wire "<<wire<<endl;
1129 }
1130}
1131
1132void HoughTrack::XhitCutWindow(double rho, int ilayer, double charge, double& cut1, double &cut2)
1133{
1134 static bool initialized = false;
1135 if(!initialized)
1136 {
1137 const int nPt=12;
1138 double rho[nPt] = {0.07, 0.056,0.045,0.038,0.0335,0.03,0.0198,0.0148,0.0074,0.00376,0.00303,0.00157};
1139 double cut1_Cgem_99[12]={-2.14,-1.36, -0.6 , -0.46, -0.35, -0.59, -0.32, -0.26, -0.22, -0.21, -0.21, -0.211};
1140 double cut2_Cgem_99[12]={ 0.5 , 1.77, 1.99, 1.94, 1.99, 1.9 , 0.31, 0.27, 0.24, 0.23, 0.23, 0.222};
1141 double cut1_ODC1_99[12]={-1.28,-0.71, -0.58, -0.62, -0.64, -0.68, -0.18, -0.14, -0.11, -0.1 , -0.1 , -0.11 };
1142 double cut2_ODC1_99[12]={ 0.9 , 0.74, 0.42, 0.37, 0.32, 0.37, 0.28, 0.25, 0.24, 0.24, 0.24, 0.23};
1143 double cut1_ODC2_99[12]={ -2.14, -1.25, -1.28, -0.86, -0.47, -0.78, -0.36, -0.44, -0.61, -0.67, -0.64, -0.82 };
1144 double cut2_ODC2_99[12]={ -1.35, 0.55 , 0.53 , 0.83 , 0.85 , 0.83 , 0.38 , 0.55 , 0.49 , 0.46 , 0.49 , 0.4};
1145 m_cut1_cgem = new TGraph(nPt, rho, cut1_Cgem_99);
1146 m_cut2_cgem = new TGraph(nPt, rho, cut2_Cgem_99);
1147 m_cut1_ODC1 = new TGraph(nPt, rho, cut1_ODC1_99);
1148 m_cut2_ODC1 = new TGraph(nPt, rho, cut2_ODC1_99);
1149 m_cut1_ODC2 = new TGraph(nPt, rho, cut1_ODC2_99);
1150 m_cut2_ODC2 = new TGraph(nPt, rho, cut2_ODC2_99);
1151 initialized = true;
1152 }
1153
1154 rho=fabs(rho);
1155 if(rho>0.07) rho=0.07;
1156 TGraph* g_cut1, *g_cut2;
1157 if(ilayer<=3) {
1158 g_cut1=m_cut1_cgem;
1159 g_cut2=m_cut2_cgem;
1160 }
1161 else if(ilayer<=19) {
1162 g_cut1=m_cut1_ODC1;
1163 g_cut2=m_cut2_ODC1;
1164 }
1165 else if(ilayer<=42) {
1166 g_cut1=m_cut1_ODC2;
1167 g_cut2=m_cut2_ODC2;
1168 }
1169 // charge < 0
1170 cut1=g_cut1->Eval(rho);
1171 cut2=g_cut2->Eval(rho);
1172 if(charge>0) {// charge > 0
1173 double cut=cut1;
1174 cut1=-1*cut2;
1175 cut2=-1*cut;
1176 }
1177 if(charge==0) {// charge = 0
1178 double cut=max(fabs(cut1),fabs(cut2));
1179 cut1=-1*cut;
1180 cut2=cut;
1181 }
1182 cut1=1.2*cut1;
1183 cut2=1.2*cut2;
1184
1185 //delete m_cut1_cgem;
1186 //delete m_cut2_cgem;
1187 //delete m_cut1_ODC1;
1188 //delete m_cut2_ODC1;
1189 //delete m_cut1_ODC2;
1190 //delete m_cut2_ODC2;
1191}
1192
1193void HoughTrack::clearHits()
1194{
1195 //m_hot.clear();
1196 m_vecHitPnt.clear();
1197 m_vecHitResidual.clear();
1198 m_vecHitChi2.clear();
1199 m_map_lay_d.clear();
1200 m_map_lay_hit.clear();
1201 m_setLayer.clear();
1202 m_unused=0;
1203 m_untried=0;
1204 m_nHitFirstHalf=0;
1205 m_nHitSecondHalf=0;
1206 m_nHitUnused_FirstHalf=0;
1207 m_nHitUnused_SecondHalf=0;
1208
1209
1210}
1211
1212// gap statistics, remove the last few hits after a big gap
1213int HoughTrack::XGapSize(std::set<int> aLaySet, int& gapMax)
1214{
1215 int nGap=0;
1216 gapMax=0;
1217 int iGapmax = -1;
1218 int nBigGap = 0;
1219 vector<int> vec_nGap;
1220 vector<int> vec_layerAfterGap;
1221 //cout<<"nGap = "<<nGap<<endl;
1222 int nLayer=aLaySet.size();
1223 if(nLayer<2) return 0;
1224 std::set<int>::iterator it = aLaySet.begin();
1225 int lastLayer = *it;
1226 // cout<<"lastLayer = "<<lastLayer<<endl;
1227 //it++;
1228 for(;it!=aLaySet.end(); it++)
1229 {
1230 int iLayer = *it;
1231 //cout<<"iLayer = "<<iLayer<<endl;
1232 if(iLayer>=8&&iLayer<=19) iLayer=iLayer-5;
1233 if(iLayer>=36&&iLayer<=42) iLayer=iLayer-21;
1234
1235 // switch(lastLayer)
1236 // {
1237 // case 2:
1238 // lastLayer=7;
1239 // break;
1240 // case 19:
1241 // lastLayer=35;
1242 // break;
1243 // default:
1244 // break;
1245 // }
1246 if(it!=aLaySet.begin()) {
1247 if(iLayer-1-lastLayer>0) {
1248 //if(nGap==0) cout<<endl;
1249 int aGap = iLayer-1-lastLayer;
1250 //cout<<aGap<<" gap(s) before layer "<<*it<<endl;
1251 nGap+=aGap;
1252 vec_nGap.push_back(aGap);
1253 vec_layerAfterGap.push_back(*it);
1254 if(aGap>2) nBigGap++;
1255 if(gapMax<aGap) {
1256 gapMax=aGap;
1257 iGapmax=vec_nGap.size()-1;
1258 }
1259 }
1260 }
1261 lastLayer=iLayer;
1262 //cout<<"nGap = "<<nGap<<endl;
1263 }
1264 // remove the last few hits if only one big gap
1265 if(nBigGap==1&&nGap/(nGap+nLayer)>0.3) {
1266 vector<double>::iterator it0_res = m_vecHitResidual.begin();
1267 vector<HoughHit*>::iterator iter0 = m_vecHitPnt.begin();
1268 vector<HoughHit*>::iterator iter = m_vecHitPnt.begin();
1269 for(; iter!=m_vecHitPnt.end(); iter++)
1270 {
1271 int iLayer = (*iter)->getLayer();
1272 if(iLayer>=vec_layerAfterGap[iGapmax])
1273 {
1274 dropHitPnt(*iter);
1275
1276 }
1277 }
1278 }
1279
1280 return nGap;
1281
1282}
1283/*
1284 HepVector HoughTrack::calCircle(vector<HoughHit*>& aVecHoughHit, double& chi2_new)
1285 {
1286//bool loc_debug = true;
1287bool loc_debug =false;
1288// trkpar: dr, phi0, kappa
1289HepVector vec_a = a();
1290HepVector trkpar(3);
1291trkpar(1)=vec_a(1);
1292trkpar(2)=vec_a(2);
1293trkpar(3)=vec_a(3);
1294double dr=trkpar(1);
1295double phi0=trkpar(2);
1296double kappa=trkpar(3);
1297double rad = m_alpha/trkpar(3);
1298double xc=(dr+rad)*cos(phi0);
1299double yc=(dr+rad)*sin(phi0);
1300if(loc_debug) cout<<" ~~~ CalculateCirclePar begin "<<endl;
1301if(loc_debug) cout<<"dr, phi0, kappa = "<<dr<<", "<<phi0<<", "<<kappa<<endl;
1302
1303// new parameters, diff
1304HepVector a_new(3), da(3);
1305
1306// initial charge
1307int charge=int(trkpar(3)/fabs(trkpar(3)));
1308
1309// U matrix = (A^T*W*A)^-1
1310HepSymMatrix U(3,0);
1311// A^T*W*(M-F(a))
1312HepVector ATWdM(3,0);
1313
1314double chi2_ini(0.0);
1315chi2_new=0.0;
1316
1317// loop the hits
1318int nHits = 0;
1319std::vector<HoughHit*>::iterator iter0=aVecHoughHit.begin();
1320std::vector<HoughHit*>::iterator iter=iter0;
1321for(; iter!=aVecHoughHit.end(); iter++)
1322{
1323//cout<<"start a hit"<<endl;
1324int sLayerType = (*iter)->getFlag();
1325if(sLayerType!=0) continue;
1326nHits++;
1327int iLayer=(*iter)->getLayer();
1328if(loc_debug) cout<<"iLayer = "<<iLayer;
1329double hit_x = (*iter)->getHitPosition().x();
1330double hit_y = (*iter)->getHitPosition().y();
1331// measurement and error
1332double m, err2, m_trkPar, dm;
1333// J_dmda for a hit
1334HepVector J(3,0);
1335if((*iter)->getHitType()==HoughHit::CGEMHIT) {// CGEM X-clusters
1336//cout<<"CGEM hit a"<<endl;
1337double r_cgem = sqrt(hit_x*hit_x+hit_y*hit_y);
1338// -- measurement estimated with the determine function
1339//cout<<"CGEM hit b"<<endl;
1340double dphi_circleInterCgem = IntersectCylinder(r_cgem);
1341double phi_circleInterCgem = x(dphi_circleInterCgem).phi();
1342//cout<<"CGEM hit c"<<endl;
1343double phi_hit = atan2(hit_y, hit_x);
1344double dphi = phi_hit-phi_circleInterCgem;
1345//cout<<"phi_hit, phi_circleInterCgem, dphi="<<phi_hit<<", "<< phi_circleInterCgem<<", "<< dphi<<endl;
1346while(dphi>pi) dphi-=2*pi;
1347while(dphi<-pi) dphi+=2*pi;
1348//cout<<"phi_hit, phi_circleInterCgem, dphi="<<phi_hit<<", "<< phi_circleInterCgem<<", "<< dphi<<endl;
1349dm=dphi*r_cgem;
1350//cout<<setprecision(5)<<"iLayer, dm = "<<iLayer<<", "<<dm<<endl;
1351if(loc_debug) {
1352cout<<", phi_hit = "<<phi_hit;
1353cout<<", phi_circleInterCgem = "<<phi_circleInterCgem;
1354cout<<", dphi_circleInterCgem = "<<dphi_circleInterCgem;
1355cout<<setprecision(5)<<", dm = "<<dm;
1356cout<<", r_cgem"<<r_cgem;
1357}
1358// -- turning angle
1359double Phi = dphi_circleInterCgem;
1360double cosPhi=cos(Phi);
1361double sinPhi=sin(Phi);
1362// -- Jacobi dm/dx
1363HepMatrix J_dmdx(1,2,0);
1364//cout<<"to fill J_dmdx(1)"<<endl;
1365double rXStrip = (*iter)->getCgemGeomSvc()->getReadoutPlane(iLayer,0)->getRX()/10;
1366J_dmdx(1,1)=-rXStrip*hit_y/(r_cgem*r_cgem);
1367//cout<<"to fill J_dmdx(2)"<<endl;
1368J_dmdx(1,2)= rXStrip*hit_x/(r_cgem*r_cgem);
1369//cout<<"filled J_dmdx(2)"<<endl;
1370// -- Jacobi dx/da
1371double dPhiDdr = -(kappa/m_alpha-cosPhi/(dr+m_alpha/kappa))/sinPhi;
1372double dPhiDkappa = -kappa*(dr*dr-r_cgem*r_cgem)*(2*m_alpha+dr*kappa)/(2*m_alpha*(m_alpha+dr*kappa)*(m_alpha+dr*kappa)*sinPhi);
1373double dxDdr = cos(phi0)+m_alpha/kappa*sin(phi0+Phi)*dPhiDdr;
1374double dyDdr = sin(phi0)-m_alpha/kappa*cos(phi0+Phi)*dPhiDdr;
1375double dxDphi0 = -dr*sin(phi0)+m_alpha/kappa*(-sin(phi0)+sin(phi0+Phi));
1376double dyDphi0 = dr*cos(phi0)+m_alpha/kappa*( cos(phi0)-cos(phi0+Phi));
1377double dxDkappa = -m_alpha/(kappa*kappa)*(cos(phi0)-cos(phi0+Phi))+m_alpha/kappa*sin(phi0+Phi)*dPhiDkappa;
1378double dyDkappa = -m_alpha/(kappa*kappa)*(sin(phi0)-sin(phi0+Phi))-m_alpha/kappa*cos(phi0+Phi)*dPhiDkappa;
1379HepMatrix J_dxda(2,3,0);
1380J_dxda(1,1) = dxDdr;
1381J_dxda(1,2) = dxDphi0;
1382J_dxda(1,3) = dxDkappa;
1383J_dxda(2,1) = dyDdr;
1384J_dxda(2,2) = dyDphi0;
1385J_dxda(2,3) = dyDkappa;
1386// -- J_dmda=dm/dx*dx/da
1387J=(J_dmdx*J_dxda).T();
1388// -- error
1389err2=0.0130*0.0130;// in cm
1390//cout<<"finish one Cgem hit"<<endl;
1391}
1392else continue;
1393
1394// else {// ODC
1395////cout<<"ODC hit a"<<endl;
1396//// -- measurement estimated with the determine function
1397//double rhit=sqrt(hit_x*hit_x+hit_y*hit_y);
1398//double phi_hit = atan2(hit_y, hit_x);
1399//double Rho = rad+trkpar(1);
1400//double l = fabs(rad+trkpar(1));
1401//double r_hitCent=sqrt((hit_x-xc)*(hit_x-xc)+(hit_y-yc)*(hit_y-yc));
1402//double d_trk = fabs(rad)-r_hitCent;
1403//double sign=0.0;
1404//if(d_trk!=0.0) sign=d_trk/fabs(d_trk);
1405//double d_measured = (*iter)->getDriftDist();
1406//dm = sign*d_measured-d_trk;
1407//// -- error
1408////err2=(*iter)->getDriftDistErr();
1409//err2=(*iter)->getMdcCalibFunSvc()->getSigma(iLayer,2,d_measured*10);// in mm
1410//err2=err2*err2/100.;// in cm
1411////cout<<"iLayer, d_measured, d_trk, dm = "<<iLayer<<", "<<sign*d_measured<<", "<<d_trk<<", "<<dm<<endl;
1412//if(loc_debug) cout<<" d_measured, d_trk, dm = "<<sign*d_measured<<", "<<d_trk<<", "<<dm;
1413//// -- turning angle
1414//HepPoint3D aPivot = pivot();
1415//Helix aHelix(aPivot, vec_a);
1416//aHelix.pivot((*iter)->getHitPosition());
1417//double phi0_new = aHelix.phi0();
1418////double Phi = dPhi((*iter)->getHitPosition());
1419//double Phi = phi0_new-phi0;
1420//while(Phi>M_PI) Phi-=2.0*M_PI;
1421//while(Phi<-M_PI) Phi+=2.0*M_PI;
1422//double cosPhi=cos(Phi);
1423//double sinPhi=sin(Phi);
1424//if(loc_debug) cout<<" tuning phi = "<<Phi;
1425//// --- Jacobi dd/da
1426//double dc=sqrt(rhit*rhit+Rho*Rho-2.*rhit*Rho*cos(phi_hit-phi0));
1427////cout<<"dc, r_hitCent, del = "<<dc<<", "<< r_hitCent<<", "<< dc-r_hitCent<<endl;
1428//double dDddr = -(Rho-rhit*cos(phi_hit-phi0))/dc;
1429//double dddphi0 = rhit*Rho*sin(phi_hit-phi0)/dc;
1430//double dddkappa= -rad/kappa;
1431//if(rad<0) dddkappa=rad/kappa;
1432//dddkappa+=rad*(Rho-rhit*cos(phi_hit-phi0))/(kappa*dc);
1433//J(1)=dDddr;
1434//J(2)=dddphi0;
1435//J(3)=dddkappa;
1436////cout<<"finish one ODC hit"<<endl;
1437//}
1438
1439// --- chi2 for initial trk parameters
1440double chi2_hit = dm*dm/err2;
1441if(loc_debug) cout<<", err = "<<sqrt(err2);
1442if(loc_debug) cout<<", chi2 = "<<chi2_hit<<endl;
1443chi2_ini+=dm*dm/err2;
1444// --- for U and ATWdM accumulating
1445for(int i=1; i<4; i++) {
1446 ATWdM(i)=ATWdM(i)+dm/err2*J(i);
1447 for(int j=1; j<4; j++) {
1448 U(i,j) = U(i,j)+J(i)*J(j)/err2;
1449 }
1450}
1451// cout<<"finish a hit at layer "<<iLayer<<endl;
1452}// loop hits
1453int ifail=99;
1454HepSymMatrix Uinv = U.inverse(ifail);
1455da=Uinv*ATWdM;
1456a_new=trkpar+da;
1457if(loc_debug) {
1458 cout<<"chi2_ini = "<<chi2_ini<<",chi2/ndof = "<<chi2_ini/(nHits-3)<<endl;
1459 cout<<"dr, phi0, kappa = "<<dr<<", "<<phi0<<", "<<kappa<<endl;
1460 cout<<"a_new = "<<a_new(1)<<", "<<a_new(2)<<", "<<a_new(3)<<endl;
1461 cout<<" CalculateCirclePar end ~~~ "<<endl;
1462}
1463return a_new;
1464}
1465//*/
1466
1467void HoughTrack::markUsedHot(vector<HoughHit*>& hitPntList, int use)
1468{
1469 vector<HoughHit*>::iterator iter = hitPntList.begin();
1470 for(; iter!=hitPntList.end(); iter++)
1471 {
1472 int useOld = (*iter)->getUse();
1473 if(use==-1&&useOld>0) continue;
1474 (*iter)->setUse(use);
1475
1476 if(use==1) {
1477 //cout<<"HoughTrack::markUsedHot add trk ID "<<getTrkID()<<endl;
1478 //(*iter)->addTrkPnt(this);
1479 (*iter)->addTrkID(getTrkID());
1480 }
1481 }
1482}
1483
1484void HoughTrack::markUsedHot(int use)
1485{
1486 vector<double>::iterator it0_res = m_vecHitResidual.begin();
1487 vector<HoughHit*>::iterator iter0 = m_vecHitPnt.begin();
1488 vector<HoughHit*>::iterator iter = m_vecHitPnt.begin();
1489 for(; iter!=m_vecHitPnt.end(); iter++)
1490 {
1491 int useOld = (*iter)->getUse();
1492 if(use==-1&&useOld>0) continue;
1493 (*iter)->setUse(use);
1494 if(use==1) {
1495 //cout<<"HoughTrack::markUsedHot(use) add trk ID "<<getTrkID()<<endl;
1496 //(*iter)->addTrkPnt(this);
1497 (*iter)->addTrkID(getTrkID());
1498 vector<double>::iterator it_res = it0_res+(iter-iter0);
1499 //cout<<"add residual "<<*it_res<<endl;
1500 (*iter)->addResid(*it_res);
1501 //vector<double> vecRes = (*iter)->getVecResid();
1502 //cout<<"VecResid.size = "<<vecRes.size()<<endl;
1503 }
1504 }
1505}
1506
1507bool HoughTrack::isAGoodCircle()
1508{
1509 bool debug = false; //debug=true;
1510 int NtotLayer = m_setLayer.size();
1511
1512 std::set<int>::iterator it = m_setLayer.begin();
1513 int lastLayer = *it;
1514 // check N_CGEM_layer, N_ODC1_nearStereo, N_ODC2_nearStereo
1515 int nCgemLayer(0), N_ODC1_nearStereo(0), N_ODC2_nearStereo(0);
1516 for(; it!=m_setLayer.end(); it++)
1517 {
1518 if((*it)<3) nCgemLayer++;
1519 else if((*it)>=16&&(*it)<=19) N_ODC1_nearStereo++;
1520 else if((*it)>=36&&(*it)<=39) N_ODC2_nearStereo++;
1521 }
1522 bool enoughVHits;
1523 //enoughVHits = nCgemLayer==3 ;
1524 enoughVHits = nCgemLayer>=2
1525 || (nCgemLayer>0 && (N_ODC1_nearStereo>=2 || N_ODC2_nearStereo>=2))
1526 //|| (N_ODC1_nearStereo>=2 && N_ODC2_nearStereo>=2 && NtotLayer>=8) ;
1527 || ( (N_ODC1_nearStereo>=2 || N_ODC2_nearStereo>=2) && NtotLayer>=8) ; // modified in 2020-08
1528
1529 // gap fraction
1530 double r_gap = 1.0*m_nGap/(m_nGap+NtotLayer);
1531 bool smallGap = r_gap<0.5;
1532 if(debug)
1533 {
1534 if(!enoughVHits) cout<<"not enough V hits, ";
1535 if(!smallGap) cout<<"gap too large 50%, "<<endl;
1536 if(m_unused<3) cout<<"N_unused = "<<m_unused<<"<3"<<endl;
1537 if(m_untried<3) cout<<"N_untried = "<<m_untried<<"<3"<<endl;
1538 }
1539
1540 // --- check conditions
1541 bool good = true;
1542 good=good && m_unused>=2; // 3 -> 2, modified 2020-10
1543 good=good && m_untried>=2;// 3 -> 2, modified 2020-10
1544 good=good && (nCgemLayer==3||NtotLayer>3);
1545 good=good && enoughVHits;
1546 good=good && smallGap;
1547 //good=good && m_maxGap<4;
1548 return good;
1549
1550}
1551
1552/*
1553 vector<double> HoughTrack::CircleFitByTaubin(vector< pair<double, double> > pos_hits)
1554 {
1555
1556 int nHits = pos_hits.size();
1557// for Taubin fit
1558int i, iter, IterMAX=99;
1559double Xi,Yi,Zi;
1560double Mz,Mxy,Mxx,Myy,Mxz,Myz,Mzz,Cov_xy,Var_z;
1561double A0,A1,A2,A22,A3,A33;
1562double Dy,xnew,x,ynew,y;
1563double DET,Xcenter,Ycenter;
1564
1565
1566// Compute x- and y- sample means (via a function in the class "data")
1567double meanX(0.), meanY(0.);
1568for (i=0; i<nHits; i++)
1569{
1570meanX+=pos_hits[i].first;
1571meanY+=pos_hits[i].second;
1572}
1573meanX/=nHits;
1574meanY/=nHits;
1575
1576
1577// computing moments
1578Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
1579for (i=0; i<nHits; i++)
1580{
1581Xi = pos_hits[i].first - meanX; // centered x-coordinates
1582Yi = pos_hits[i].second - meanY; // centered y-coordinates
1583Zi = Xi*Xi + Yi*Yi;
1584
1585Mxy += Xi*Yi;
1586Mxx += Xi*Xi;
1587Myy += Yi*Yi;
1588Mxz += Xi*Zi;
1589Myz += Yi*Zi;
1590Mzz += Zi*Zi;
1591}
1592Mxx /= nHits;
1593Myy /= nHits;
1594Mxy /= nHits;
1595Mxz /= nHits;
1596Myz /= nHits;
1597Mzz /= nHits;
1598// computing coefficients of the characteristic polynomial
1599
1600Mz = Mxx + Myy;
1601Cov_xy = Mxx*Myy - Mxy*Mxy;
1602Var_z = Mzz - Mz*Mz;
1603A3 = 4*Mz;
1604A2 = -3*Mz*Mz - Mzz;
1605A1 = Var_z*Mz + 4*Cov_xy*Mz - Mxz*Mxz - Myz*Myz;
1606A0 = Mxz*(Mxz*Myy - Myz*Mxy) + Myz*(Myz*Mxx - Mxz*Mxy) - Var_z*Cov_xy;
1607A22 = A2 + A2;
1608A33 = A3 + A3 + A3;
1609
1610// finding the root of the characteristic polynomial
1611// using Newton's method starting at x=0
1612// (it is guaranteed to converge to the right root)
1613
1614for (x=0.,y=A0,iter=0; iter<IterMAX; iter++) // usually, 4-6 iterations are enough
1615{
1616Dy = A1 + x*(A22 + A33*x);
1617xnew = x - y/Dy;
1618if ((xnew == x)||(!isfinite(xnew))) break;
1619ynew = A0 + xnew*(A1 + xnew*(A2 + xnew*A3));
1620if (abs(ynew)>=abs(y)) break;
1621x = xnew; y = ynew;
1622}
1623//if(iter==IterMAX) cout<<"IterMAX reached!"<<endl;
1624
1625// computing paramters of the fitting circle
1626
1627DET = x*x - x*Mz + Cov_xy;
1628Xcenter = (Mxz*(Myy - x) - Myz*Mxy)/DET/2;
1629Ycenter = (Myz*(Mxx - x) - Mxz*Mxy)/DET/2;
1630
1631// assembling the output
1632
1633double xc_fit = Xcenter + meanX;
1634double yc_fit = Ycenter + meanY;
1635double r_fit = sqrt(Xcenter*Xcenter + Ycenter*Ycenter + Mz);
1636//cout<<"fitted trk par: xc, yc, r = "<<xc_fit<<" ,"<<yc_fit<<", "<<r_fit<<endl;
1637//cout<<"after "<<iter<<" iterations! "<<endl;
1638//cout<<"CircleFitByTaubin "<<iter<<" iterations! "<<endl;
1639
1640vector<double> output;
1641output.push_back(xc_fit);
1642output.push_back(yc_fit);
1643output.push_back(r_fit);
1644output.push_back(iter);
1645
1646return output;
1647}
1648//*/
1649/*
1650 vector<double> HoughTrack::fitCircle_Taubin(vector<HoughHit*>& aVecHoughHit)
1651 {
1652 m_circleFitStat=0;
1653 bool loc_debug = true; loc_debug=false;
1654// trkpar: dr, phi0, kappa
1655HepVector vec_a = a();
1656HepVector trkpar(3);
1657trkpar(1)=vec_a(1);
1658trkpar(2)=vec_a(2);
1659trkpar(3)=vec_a(3);
1660double dr=trkpar(1);
1661double phi0=trkpar(2);
1662double kappa=trkpar(3);
1663double rad = m_alpha/trkpar(3);
1664double xc=(dr+rad)*cos(phi0);
1665double yc=(dr+rad)*sin(phi0);
1666if(loc_debug) {
1667cout<<" ~~~ CalculateCirclePar begin "<<endl;
1668cout<<"dr, phi0, kappa = "<<dr<<", "<<phi0<<", "<<kappa<<endl;
1669cout<<"xc, yc, rad = "<<xc<<", "<<yc<<", "<<rad<<endl;
1670}
1671
1672// initial charge
1673int charge=int(trkpar(3)/fabs(trkpar(3)));
1674
1675// loop the hits to give input vector for CircleFitByTaubin
1676vector< pair<double, double> > pos_hits;
1677int nHits = 0;
1678std::vector<HoughHit*>::iterator iter0=aVecHoughHit.begin();
1679std::vector<HoughHit*>::iterator iter=iter0;
1680for(; iter!=aVecHoughHit.end(); iter++)
1681{
1682//cout<<"start a hit"<<endl;
1683int sLayerType = (*iter)->getFlag();
1684if(sLayerType!=0) continue;
1685nHits++;
1686int iLayer=(*iter)->getLayer();
1687if(loc_debug) cout<<"iLayer = "<<iLayer<<" ";
1688double hit_x = (*iter)->getHitPosition().x();
1689double hit_y = (*iter)->getHitPosition().y();
1690// measurement and error
1691double m, err2, m_trkPar, dm;
1692if((*iter)->getHitType()==HoughHit::CGEMHIT) {// CGEM X-clusters
1693pair<double, double> aHitPos(hit_x,hit_y);
1694pos_hits.push_back(aHitPos);
1695//cout<<"CGEM hit a"<<endl;
1696if(loc_debug) {
1697cout<<" hit_x, hit_y = "<<hit_x<<", "<<hit_y<<endl;
1698}
1699}
1700else {// ODC
1701//cout<<"ODC hit a"<<endl;
1702// -- turning angle
1703HepPoint3D aPivot = pivot();
1704Helix aHelix(aPivot, vec_a);
1705aHelix.bFieldZ(-(aHelix.bFieldZ()));
1706aHelix.pivot((*iter)->getHitPosition());
1707double phi0_new = aHelix.phi0();
1708double dr_new = aHelix.dr();
1709double d_measured = (*iter)->getDriftDist();
1710d_measured=dr_new/fabs(dr_new)*d_measured;
1711double x_new = hit_x+d_measured*cos(phi0_new);
1712double y_new = hit_y+d_measured*sin(phi0_new);
1713if(loc_debug) {
1714cout<<" hit_x, hit_y, d_measured, x_new, y_new = "<<hit_x<<", "<<hit_y<<", "<<d_measured<<", "
1715<<x_new<<", "<<y_new<<endl;
1716cout<<"new pos: "<<aHelix.x(0.)<<endl;
1717}
1718pair<double, double> aHitPos(x_new,y_new);
1719//pair<double, double> aHitPos(hit_x,hit_y);
1720pos_hits.push_back(aHitPos);
1721//cout<<"finish one ODC hit"<<endl;
1722}
1723// cout<<"finish a hit at layer "<<iLayer<<endl;
1724}// loop hits
1725//if(loc_debug) {
1726//cout<<"chi2_ini = "<<chi2_ini<<",chi2/ndof = "<<chi2_ini/(nHits-3)<<endl;
1727//cout<<"dr, phi0, kappa = "<<dr<<", "<<phi0<<", "<<kappa<<endl;
1728//cout<<"a_new = "<<a_new(1)<<", "<<a_new(2)<<", "<<a_new(3)<<endl;
1729//cout<<" CalculateCirclePar end ~~~ "<<endl;
1730//}
1731
1732vector<double> par_new = CircleFitByTaubin(pos_hits);
1733double xc_fit = par_new[0];
1734double yc_fit = par_new[1];
1735double r_fit = par_new[2];
1736double dr_fit = sqrt(xc_fit*xc_fit+yc_fit*yc_fit)-r_fit;
1737double phi0_fit = atan2(yc_fit,xc_fit);
1738double kappa_fit= m_alpha/r_fit;
1739if(rad<0) {
1740 dr_fit = -dr_fit;
1741 phi0_fit+=M_PI;
1742 kappa_fit = -kappa_fit;
1743}
1744while(phi0_fit< 0 ) phi0_fit+=2*M_PI;
1745while(phi0_fit> 2*M_PI) phi0_fit-=2*M_PI;
1746vector<double> par_fit;
1747par_fit.push_back(dr_fit);
1748par_fit.push_back(phi0_fit);
1749par_fit.push_back(kappa_fit);
1750if(loc_debug) cout<<"after fit: dr, phi0, kappa = "<<dr_fit<<", "<<phi0_fit<<", "<<kappa_fit<<endl;
1751return par_fit;
1752}
1753//*/
1754/*
1755 double HoughTrack::chi2_circle(vector<HoughHit*>& aVecHoughHit)
1756 {
1757//bool loc_debug=true;
1758bool loc_debug=false;
1759double chi2=0.0;
1760
1761// trkpar: dr, phi0, kappa
1762HepVector vec_a = a();
1763HepVector trkpar(3);
1764trkpar(1)=vec_a(1);
1765trkpar(2)=vec_a(2);
1766trkpar(3)=vec_a(3);
1767double dr=trkpar(1);
1768double phi0=trkpar(2);
1769double kappa=trkpar(3);
1770double rad = m_alpha/trkpar(3);
1771double xc=(dr+rad)*cos(phi0);
1772double yc=(dr+rad)*sin(phi0);
1773
1774// loop the hits
1775int nHits = 0;
1776std::vector<HoughHit*>::iterator iter0=aVecHoughHit.begin();
1777std::vector<HoughHit*>::iterator iter=iter0;
1778for(; iter!=aVecHoughHit.end(); iter++)
1779{
1780//cout<<"start a hit"<<endl;
1781int sLayerType = (*iter)->getFlag();
1782if(sLayerType!=0) continue;
1783nHits++;
1784int iLayer=(*iter)->getLayer();
1785if(loc_debug) cout<<"iLayer = "<<iLayer;
1786double hit_x = (*iter)->getHitPosition().x();
1787double hit_y = (*iter)->getHitPosition().y();
1788// measurement and error
1789double m, err2, m_trkPar, dm;
1790if((*iter)->getHitType()==HoughHit::CGEMHIT) {// CGEM X-clusters
1791//cout<<"CGEM hit a"<<endl;
1792double r_cgem = sqrt(hit_x*hit_x+hit_y*hit_y);
1793// -- measurement estimated with the determine function
1794//cout<<"CGEM hit b"<<endl;
1795double dphi_circleInterCgem = IntersectCylinder(r_cgem);
1796double phi_circleInterCgem = x(dphi_circleInterCgem).phi();
1797//cout<<"CGEM hit c"<<endl;
1798double phi_hit = atan2(hit_y, hit_x);
1799double dphi = phi_hit-phi_circleInterCgem;
1800//cout<<"phi_hit, phi_circleInterCgem, dphi="<<phi_hit<<", "<< phi_circleInterCgem<<", "<< dphi<<endl;
1801while(dphi>pi) dphi-=2*pi;
1802while(dphi<-pi) dphi+=2*pi;
1803//cout<<"phi_hit, phi_circleInterCgem, dphi="<<phi_hit<<", "<< phi_circleInterCgem<<", "<< dphi<<endl;
1804dm=dphi*r_cgem;
1805//cout<<setprecision(5)<<"iLayer, dm = "<<iLayer<<", "<<dm<<endl;
1806if(loc_debug) {
1807cout<<", phi_hit = "<<phi_hit;
1808cout<<", phi_circleInterCgem = "<<phi_circleInterCgem;
1809cout<<", dphi_circleInterCgem = "<<dphi_circleInterCgem;
1810cout<<setprecision(5)<<", dm = "<<dm;
1811cout<<", r_cgem"<<r_cgem;
1812}
1813// -- error
1814err2=0.0130*0.0130;// in cm
1815//cout<<"finish one Cgem hit"<<endl;
1816}
1817else {// ODC
1818//cout<<"ODC hit a"<<endl;
1819// -- measurement estimated with the determine function
1820double rhit=sqrt(hit_x*hit_x+hit_y*hit_y);
1821double phi_hit = atan2(hit_y, hit_x);
1822double Rho = rad+trkpar(1);
1823double l = fabs(rad+trkpar(1));
1824double r_hitCent=sqrt((hit_x-xc)*(hit_x-xc)+(hit_y-yc)*(hit_y-yc));
1825double d_trk = fabs(rad)-r_hitCent;
1826double sign=0.0;
1827if(d_trk!=0.0) sign=d_trk/fabs(d_trk);
1828double d_measured = (*iter)->getDriftDist();
1829dm = sign*d_measured-d_trk;
1830// -- error
1831//err2=(*iter)->getDriftDistErr();
1832err2=(*iter)->getMdcCalibFunSvc()->getSigma(iLayer,2,d_measured*10);// in mm
1833err2=err2*err2/100.;// in cm
1834//cout<<"iLayer, d_measured, d_trk, dm = "<<iLayer<<", "<<sign*d_measured<<", "<<d_trk<<", "<<dm<<endl;
1835if(loc_debug) cout<<" d_measured, d_trk, dm = "<<sign*d_measured<<", "<<d_trk<<", "<<dm;
1836//cout<<"finish one ODC hit"<<endl;
1837}
1838// --- chi2 for initial trk parameters
1839double chi2_hit = dm*dm/err2;
1840if(loc_debug) cout<<", err = "<<sqrt(err2);
1841if(loc_debug) cout<<", chi2 = "<<chi2_hit<<endl;
1842//if((*iter)->getHitType()==HoughHit::CGEMHIT) // CGEM X-clusters
1843chi2+=dm*dm/err2;
1844
1845// cout<<"finish a hit at layer "<<iLayer<<endl;
1846}// loop hits
1847
1848return chi2;
1849
1850}
1851//*/
1852void HoughTrack::releaseSelHits()
1853{
1854 vector<HoughHit*>::iterator it = m_vecHitPnt.begin();
1855 for(; it!=m_vecHitPnt.end(); it++)
1856 {
1857 (*it)->dropTrkID(m_trkID);// erase trk id and residual from the Hough hit
1858 }
1859}
1860
1861void HoughTrack::dropHitPnt(HoughHit* aHitPnt)
1862{
1863 vector<HoughHit*>::iterator it0 = m_vecHitPnt.begin();
1864 vector<HoughHit*>::iterator result = find(m_vecHitPnt.begin(), m_vecHitPnt.end(), aHitPnt);
1865 if(result!=m_vecHitPnt.end()) {
1866 //cout<<"remove hit "<<(*result)->getHitID();
1867 m_vecHitPnt.erase(result);
1868 vector<double>::iterator itRes = m_vecHitResidual.begin()+(result-it0);
1869 if(itRes!=m_vecHitResidual.end()) {
1870 //cout<<" with residual "<<*itRes<<" from trk "<<getTrkID()<<endl;
1871 m_vecHitResidual.erase(itRes);
1872 }
1873 }
1874}
1875
1876void HoughTrack::dropVHitPnt(HoughHit* aHitPnt)
1877{
1878 vector<HoughHit*>::iterator it0 = m_vecStereoHitPnt.begin();
1879 vector<HoughHit*>::iterator result = find(m_vecStereoHitPnt.begin(), m_vecStereoHitPnt.end(), aHitPnt);
1880 if(result!=m_vecStereoHitPnt.end()) {
1881 //cout<<"remove hit "<<(*result)->getHitID();
1882 m_vecStereoHitPnt.erase(result);
1883 }
1884}
1885
1886void HoughTrack::dropRedundentCgemXHits()
1887{
1888 //cout<<"begin HoughTrack::dropRedundentCgemXHits"<<endl;
1889 vector<HoughHit*> hitPnt_cgem[3];
1890 HoughHit* hitPnt_cgem_keep[3]={NULL, NULL, NULL};
1891 double res_cgem_keep[3]={999., 999., 999.};
1892 HoughHit* hitPnt_cgem_resMin[3]={NULL, NULL, NULL};
1893 double res_cgem_min[3]={999., 999., 999.};
1894 vector<HoughHit*>::iterator iter0 = m_vecHitPnt.begin();
1895 vector<HoughHit*>::iterator iter = m_vecHitPnt.begin();
1896 vector<double>::iterator itRes = m_vecHitResidual.begin();
1897 for(; iter!=m_vecHitPnt.end(); iter++)
1898 {
1899 int iLayer = (*iter)->getLayer();
1900 if(iLayer<3) {
1901 hitPnt_cgem[iLayer].push_back(*iter);
1902 int idx = iter-iter0;
1903 double res = m_vecHitResidual[idx];
1904 vector<double> vecRes = (*iter)->getVecResid();
1905 int nTrk = vecRes.size();
1906 if(nTrk==1&&fabs(res)<fabs(res_cgem_keep[iLayer])) {
1907 res_cgem_keep[iLayer] = res;
1908 hitPnt_cgem_keep[iLayer] = *iter;
1909 }
1910 if(fabs(res)<fabs(res_cgem_min[iLayer])) {
1911 res_cgem_min[iLayer] = res;
1912 hitPnt_cgem_resMin[iLayer] = *iter;
1913 }
1914 }
1915 }
1916 for(int iLayer=0; iLayer<3; iLayer++)
1917 {
1918 int nHits = hitPnt_cgem[iLayer].size();
1919 if(nHits>1)
1920 {
1921 HoughHit* hitToKeep;
1922 if(hitPnt_cgem_keep[iLayer]!=NULL) hitToKeep = hitPnt_cgem_keep[iLayer];
1923 else hitToKeep = hitPnt_cgem_resMin[iLayer];
1924 //cout<<"HoughTrack::dropRedundentCgemXHits: "<<endl;
1925 //cout<<" keep hit "<<hitToKeep->getHitID()
1926 //<<", layer "<<hitToKeep->getLayer()
1927 //<<" at ("<<hitToKeep->getHitPosition().x()
1928 //<<", "<<hitToKeep->getHitPosition().y()<<")"<<endl;
1929 for(int iHit=0; iHit<nHits; iHit++)
1930 {
1931 HoughHit* aHoughHit = hitPnt_cgem[iLayer][iHit];
1932 if(aHoughHit==hitToKeep) continue;
1933 //cout<<" remove hit "<<aHoughHit->getHitID()
1934 //<<", layer "<<aHoughHit->getLayer()
1935 //<<" at ("<<aHoughHit->getHitPosition().x()
1936 //<<", "<<aHoughHit->getHitPosition().y()<<")"<<endl;
1937 dropHitPnt(aHoughHit);
1938 //aHoughHit->dropTrkID(m_trkID);
1939 }
1940 }
1941 }
1942 //cout<<"end HoughTrack::dropRedundentCgemXHits"<<endl;
1943}
1944
1945void HoughTrack::dropRedundentCgemVHits()
1946{
1947 //cout<<"begin HoughTrack::dropRedundentCgemXHits"<<endl;
1948 vector<HoughHit*> hitPnt_cgem[3];
1949 HoughHit* hitPnt_cgem_resMin[3]={NULL, NULL, NULL};
1950 double res_cgem_min[3]={999., 999., 999.};
1951 vector<HoughHit*>::iterator iter = m_vecStereoHitPnt.begin();
1952 for(; iter!=m_vecStereoHitPnt.end(); iter++)
1953 {
1954 int iLayer = (*iter)->getLayer();
1955 if(iLayer<3) {
1956 hitPnt_cgem[iLayer].push_back(*iter);
1957 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
1958 double res = (*iter)->residual(this, positionOntrack, positionOnDetector);
1959 if(fabs(res)<fabs(res_cgem_min[iLayer])) {
1960 res_cgem_min[iLayer] = res;
1961 hitPnt_cgem_resMin[iLayer] = *iter;
1962 }
1963 }
1964 }
1965 for(int iLayer=0; iLayer<3; iLayer++)
1966 {
1967 int nHits = hitPnt_cgem[iLayer].size();
1968 if(nHits>1)
1969 {
1970 HoughHit* hitToKeep = hitPnt_cgem_resMin[iLayer];
1971 for(int iHit=0; iHit<nHits; iHit++)
1972 {
1973 HoughHit* aHoughHit = hitPnt_cgem[iLayer][iHit];
1974 if(aHoughHit==hitToKeep) continue;
1975 dropVHitPnt(aHoughHit);
1976 }
1977 }
1978 }
1979 //cout<<"end HoughTrack::dropRedundentCgemXHits"<<endl;
1980}
1981
1982
1983int HoughTrack::getNHitsShared()
1984{
1985 int nShared = 0;
1986 vector<HoughHit*>::iterator iter = m_vecHitPnt.begin();
1987 for(; iter!=m_vecHitPnt.end(); iter++)
1988 {
1989 vector<double> vecRes = (*iter)->getVecResid();
1990 int nTrk = vecRes.size();
1991 if(nTrk>1) nShared++;
1992 }
1993 return nShared;
1994}
1995
1996int HoughTrack::findVHot(vector<HoughHit*>& hitList, int charge, int maxLayer, double cutFactor)
1997{
1998 bool printFlag(false); //printFlag=true;
1999 int nHot(0);
2000 m_vecStereoHitPnt.clear();
2001 for(vector<HoughHit*>::iterator hitIt = hitList.begin(); hitIt != hitList.end(); hitIt++){
2002 /*
2003 if(0 == flag){
2004 if((*hitIt)->getFlag() != 0) continue;
2005 if(judgeCharge(*hitIt)*charge<0) continue;// charge requirement
2006 //if((*hitIt)->getPairHit()==NULL)continue;//FIXME
2007 //if((*hitIt)->getPairHit()->getHalfCircle()!=1)continue;//FIXME
2008 double cut1, cut2;
2009 int iLayer = (*hitIt)->getLayer();
2010 int used = (*hitIt)->getUse();
2011 XhitCutWindow(m_rho, iLayer, charge, cut1, cut2);
2012 double res = driftDistRes(*hitIt);
2013 map<int,double>::iterator it_map;
2014 it_map=m_map_lay_d.find(iLayer);
2015 if(it_map==m_map_lay_d.end()||fabs(res)<fabs(it_map->second))
2016 {
2017 m_map_lay_d[iLayer]=res;
2018 m_map_lay_hit[iLayer]=(*hitIt);
2019 }
2020 if(res>cut1&&res<cut2)
2021 {
2022 //cout<<" selected!";
2023 //if(iLayer>3)
2024 m_vecHitPnt.push_back(*hitIt);
2025 m_vecHitResidual.push_back(res);
2026 //m_hot.push_back((*(*hitIt)));
2027 m_chi2 =+ res*res;
2028 m_setLayer.insert(iLayer);
2029 if(used==0) m_untried++;
2030 if(used==0||used==-1) m_unused++;
2031 if(judgeHalf(*hitIt)>0) {
2032 m_nHitFirstHalf++;
2033 if(used<=0) m_nHitUnused_FirstHalf++;
2034 }
2035 else {
2036 m_nHitSecondHalf++;
2037 if(used<=0) m_nHitUnused_SecondHalf++;
2038 }
2039 nHot++;
2040 }// within window
2041 //cout<<endl;
2042 }// X-view
2043 //*/
2044 //else {
2045 int layer = (*hitIt)->getLayer();
2046 int wire = (*hitIt)->getWire();
2047 if((*hitIt)->getFlag() == 0) continue;
2048 if((*hitIt)->getHitType()==HoughHit::CGEMHIT){
2049 if(calculateZ_S(*hitIt)==0){
2050 continue;
2051 }
2052 }
2053 //if((*hitIt)->getPairHit()==NULL)continue;//FIXME
2054 //if((*hitIt)->getPairHit()->getHalfCircle()>1)continue;//FIXME
2055 double Pt = fabs(pt());
2056 double cut[2] = {0};
2057 if(layer<3){
2058 if(m_cut[0][layer+3]!=NULL)cut[0]=m_cut[0][layer+3]->Eval(Pt);
2059 if(m_cut[1][layer+3]!=NULL)cut[1]=m_cut[1][layer+3]->Eval(Pt);
2060 }else{
2061 if(m_cut[0][layer]!=NULL)cut[0]=m_cut[0][layer]->Eval(Pt);
2062 if(m_cut[1][layer]!=NULL)cut[1]=m_cut[1][layer]->Eval(Pt);
2063 }
2064 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
2065 double res = (*hitIt)->residual(this, positionOntrack, positionOnDetector);
2066 if(judgeCharge(positionOntrack.x(),positionOntrack.y())*charge<0) continue;// charge requirement
2067 bool isNewHot(true);
2068 for(vector<HoughHit*>::iterator hotIt = m_vecStereoHitPnt.begin(); hotIt != m_vecStereoHitPnt.end(); hotIt++){
2069 if((*hitIt)->getHitID() == (*hotIt)->getHitID()){
2070 isNewHot = false;
2071 break;
2072 }
2073 }
2074 if(!isNewHot){
2075 //cout<<" OK!";
2076 //cout<<endl;
2077 continue;
2078 }
2079 if(printFlag)(*hitIt)->print();
2080 if(printFlag)cout<<"res:"<<setw(12)<<res;
2081 if(printFlag)cout<<" win: "<<setw(5)<<cut[0]<<" ~ "<<setw(5)<<cut[1];
2082 double ext=cutFactor;
2083 if(ext*cut[0]<res&&res<ext*cut[1]){
2084 if((*hitIt)->getLayer()>=maxLayer)
2085 {
2086 //if(printFlag) cout<<;
2087 continue;
2088 }
2089 m_vecStereoHitPnt.push_back(*hitIt);
2090 if(printFlag)cout<<" selected!";
2091 nHot++;
2092 /*
2093 if(maxLayer<20){
2094 if((*hitIt)->getLayer()>=maxLayer)continue;
2095 //HepPoint3D szz(minS,zTrk,zHit);
2096 //(*hitIt)->setHitPosition(szz);
2097 m_vecStereoHitPnt.push_back(*hitIt);
2098 //m_vecStereoHitRes.push_back(minRes);
2099 nHot++;
2100 //cout<<" OK!";
2101 //(*hitIt)->print();
2102 }
2103 else{
2104 //if((*hitIt)->getLayer()!=maxLayer)continue;
2105 if((*hitIt)->getLayer()>=maxLayer)continue;
2106 //HepPoint3D szz(minS,zTrk,zHit);
2107 //(*hitIt)->setHitPosition(szz);
2108 m_vecStereoHitPnt.push_back(*hitIt);
2109 //m_vecStereoHitRes.push_back(minRes);
2110 nHot++;
2111 //cout<<" OK!";
2112 //(*hitIt)->print();
2113 }
2114 //*/
2115 }
2116 if(printFlag)cout<<endl;
2117 //if(layer<3)cout<<endl;
2118 //}// V-view
2119 }// loop hits
2120 return nHot;
2121}
2122
2123//type:0->X cluster + Axial hit,
2124//type:1->XV cluster + stereo hits,
2125//type:2->XV cluster + Axial hit + stereo hits,
2126//type:3->X cluster + XV cluster + Axial hit + stereo hits
2127vector<HoughHit*> HoughTrack::getHotList(int type)
2128{
2129 vector<HoughHit*> hotList;
2130 vector<HoughHit*> axialHot = getVecHitPnt();
2131 vector<HoughHit*> stereoHot = getVecStereoHitPnt();
2132
2133 for(vector<HoughHit*>::iterator hotIt = axialHot.begin(); hotIt != axialHot.end(); hotIt++){
2134 if(type==1)continue;
2135 if(type==2&&(**hotIt).getHitType()==HoughHit::CGEMHIT)continue;
2136 hotList.push_back(*hotIt);
2137 }
2138 for(vector<HoughHit*>::iterator hotIt = stereoHot.begin(); hotIt != stereoHot.end(); hotIt++){
2139 if(type>0)hotList.push_back(*hotIt);
2140 }
2141
2142 /*
2143 // --- CGEM XV cluster
2144 for(vector<HoughHit*>::iterator hotIt = stereoHot.begin(); hotIt != stereoHot.end(); hotIt++){
2145 if((**hotIt).getLayer()<3)hotList.push_back(*hotIt);
2146 //cout<<(**hotIt).getHitID()<<" "<<(**hotIt).getLayer()<<endl;
2147 }
2148
2149 // --- MDC hits, layer 9-20
2150 for(vector<HoughHit*>::iterator hotIt = axialHot.begin(); hotIt != axialHot.end(); hotIt++){
2151 //if(m_fitFlag==1&&(**hotIt).getLayer()<3)hot.push_back(**hotIt);
2152 if((**hotIt).getLayer()>3&&(**hotIt).getLayer()<20)hotList.push_back(*hotIt);
2153 //cout<<(**hotIt).getHitID()<<" "<<(**hotIt).getLayer()<<endl;
2154 }
2155
2156 // --- MDC hits, layer 21-36
2157 for(vector<HoughHit*>::iterator hotIt = stereoHot.begin(); hotIt != stereoHot.end(); hotIt++){
2158 if((**hotIt).getLayer()>=20&&(**hotIt).getLayer()<36)hotList.push_back(*hotIt);
2159 //int layer = (*hotIt)->getLayer();
2160 //int wire = (*hotIt)->getWire();
2161 //cout<<"("<<setw(2)<<layer<<","<<setw(3)<<wire<<") ";
2162 //cout<<(**hotIt).residual(&(*trkIT));
2163 //cout<<endl;
2164 //cout<<(**hotIt).getHitID()<<" "<<(**hotIt).getLayer()<<endl;
2165 }
2166
2167 // --- MDC hits, layer 37-43
2168 for(vector<HoughHit*>::iterator hotIt = axialHot.begin(); hotIt != axialHot.end(); hotIt++){
2169 if((**hotIt).getLayer()>=36)hotList.push_back(*hotIt);
2170 //cout<<(**hotIt).getHitID()<<" "<<(**hotIt).getLayer()<<endl;
2171 }
2172 //cout<<"nHot: "<<hot.size()<<endl;
2173 // */
2174
2175 sortHot(hotList);
2176 return hotList;
2177}
2178
2179TrkErrCode HoughTrack::fitCircle(const MdcDetector* mdcDetector, TrkContextEv* trkContextEv, double bunchT0)
2180{
2181 bool debug=false; //debug=true;
2182 m_circleFitStat=0;
2183 double d0 = -m_dr;
2184 double phi0 = m_phi0+M_PI/2;
2185 while(phi0>M_PI)phi0-=2*M_PI;
2186 while(phi0<-M_PI)phi0+=2*M_PI;
2187 double omega = m_kappa/fabs(alpha());
2188 //double z0 = m_dz;
2189 //double tanl = m_tanl;
2190 double z0 = 0;
2191 double tanl = 0;
2192 //cout<<"hough params:"<<endl;
2193 //cout
2194 //<<setw(15)<<d0
2195 //<<setw(15)<<phi0
2196 //<<setw(15)<<omega
2197 //;
2198 //cout<<endl;
2199 TrkExchangePar circleTrk(d0,phi0,omega,z0,tanl);
2200 TrkCircleMaker circlefactory;
2201 double chisum =0.;
2202 TrkRecoTrk* trkRecoTrk = circlefactory.makeTrack(circleTrk, chisum, *trkContextEv, bunchT0);
2203 bool permitFlips = true;
2204 bool lPickHits = true;
2205 circlefactory.setFlipAndDrop(*trkRecoTrk, permitFlips, lPickHits);
2206 TrkHitList* trkHitList = trkRecoTrk->hits();
2207
2208 //---convert hits
2209 vector<MdcHit*> mdcHitVector;
2210 mdcHitVector.clear();
2211 vector<CgemHitOnTrack*> cgemHitVector;
2212 cgemHitVector.clear();
2213
2214 for(vector<HoughHit*>::iterator iter = m_vecHitPnt.begin(); iter != m_vecHitPnt.end(); iter++){
2215 HoughHit* hitIt = (*iter);
2216 if(hitIt->getFlag()!=0)continue;
2217 //hitIt->print();
2218 if(hitIt->getHitType()==HoughHit::MDCHIT || hitIt->getHitType()==HoughHit::MDCMCHIT){
2219 const MdcDigi* mdcDigi = hitIt->getDigi();
2220 MdcHit *mdcHit = new MdcHit(mdcDigi,mdcDetector);
2221 mdcHitVector.push_back(mdcHit);
2222 mdcHit->setCalibSvc(hitIt->getMdcCalibFunSvc());
2223 mdcHit->setCountPropTime(true);
2224 //if(mdcHit->driftDist(m_bunchT0,0)>ddCut)continue;//FIXME
2225 int newAmbig = 0;
2226 MdcRecoHitOnTrack mdcRecoHitOnTrack(*mdcHit, newAmbig, bunchT0*1.e9);
2227 MdcHitOnTrack* mdcHitOnTrack = &mdcRecoHitOnTrack;
2228 HepPoint3D midPoint = hitIt->getHitPosition();
2229 double fltLength = flightLength(midPoint);
2230 //double fltLength = flightLength(hitIt->getHitPosition());
2231 mdcHitOnTrack->setFltLen(fltLength);
2232 trkHitList->appendHot(mdcHitOnTrack);
2233 }else if(hitIt->getHitType()==HoughHit::CGEMHIT || hitIt->getHitType()==HoughHit::CGEMMCHIT){
2234 if(m_useCgemInGlobalFit<1)continue;
2235 CgemHitOnTrack* cgemHitOnTrack = new CgemHitOnTrack(*(hitIt->getCgemCluster()),bunchT0*1.e9);
2236 cgemHitVector.push_back(cgemHitOnTrack);
2237 cgemHitOnTrack->setCgemGeomSvc(hitIt->getCgemGeomSvc());
2238 cgemHitOnTrack->setCgemCalibFunSvc(hitIt->getCgemCalibFunSvc());
2239 trkHitList->appendHot(cgemHitOnTrack);
2240 }
2241 }
2242 //---fitting
2243 TrkHitList* qhits = trkRecoTrk->hits();
2244 //cout<<"num of qhits to fit = "<<qhits->nHit()<<endl;
2245 //cout<<"before TrkHitList->fit()"<<endl;
2246 TrkErrCode err=qhits->fit();
2247 //cout<<"after TrkHitList->fit()"<<endl;
2248
2249 if(err.success()){
2250 //cout<<" circle fits well "<<endl;
2251 m_circleFitStat=1;
2252 const TrkFit* fitResult = trkRecoTrk->fitResult();
2253 TrkExchangePar par = fitResult->helix(0.);
2254 //if(fabs(m_kappa-par.omega()*fabs(alpha()))/fabs(m_kappa)>0.3)cout<<"fit divergent"<<endl;
2255 m_dr = -par.d0();
2256 m_phi0 = par.phi0()-M_PI/2;
2257 m_kappa = par.omega()*fabs(alpha());
2258 m_dz = par.z0();
2259 m_tanl = par.tanDip();
2260
2261 //cout<<"phi0 = "<<m_phi0<<endl;
2262 while(m_phi0>2.*M_PI) m_phi0-=2.*M_PI;
2263 while(m_phi0<0.) m_phi0+=2.*M_PI;
2264 //cout<<"phi0 = "<<m_phi0<<endl;
2265 //HepVector hepVector(5,0);
2266 //hepVector(1) = m_dr;
2267 //hepVector(2) = m_phi0;
2268 //hepVector(3) = m_kappa;
2269 //hepVector(4) = m_dz;
2270 //hepVector(5) = m_tanl;
2271 //a(hepVector);
2272 updateHelix();
2273 //cout<<"update Helix"<<endl;
2274
2275 //m_trkRecoTrk = trkRecoTrk;
2276 m_chi2 = fitResult->chisq();
2277 if(debug) cout<<"chi2 obtained "<<m_chi2<<endl;
2278 //cout<<"fitting params:"<<endl;
2279 //cout
2280 //<<setw(15)<<par.d0()
2281 //<<setw(15)<<par.phi0()
2282 //<<setw(15)<<par.omega()
2283 //;
2284 //cout<<endl;
2285
2286 //m_vecMdcDigiPnt.clear();
2287 //m_vecCgemClusterPnt.clear();
2288
2289 int nHotKept = 0;
2290 for(TrkHitList::hot_iterator hot = qhits->begin();hot!=qhits->end();hot++)
2291 {
2292 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ )
2293 {
2294 const MdcRecoHitOnTrack* recoHot;
2295 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2296 if(!(recoHot->isActive())) continue;
2297 nHotKept++;
2298 //const MdcDigi* aMdcDigi = recoHot->mdcHit()->digi();
2299 //m_vecMdcDigiPnt.push_back(aMdcDigi);
2300 }
2301 else if(typeid(*hot)==typeid(CgemHitOnTrack))
2302 {
2303 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2304 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2305 int clusterid = recCgemCluster->getclusterid();
2306 int stat = recoHot->isActive();
2307 if(stat==0) continue;//
2308 nHotKept++;
2309 //m_vecCgemClusterPnt.push_back(recCgemCluster);
2310 }
2311 }
2312 if(debug) cout<<nHotKept<<" hits kept after TrkFitter"<<endl;
2313 }
2314 for(vector<MdcHit*>::iterator iter = mdcHitVector.begin(); iter != mdcHitVector.end(); iter++){
2315 delete *iter;
2316 }
2317 for(vector<CgemHitOnTrack*>::iterator iter = cgemHitVector.begin(); iter != cgemHitVector.end(); iter++){
2318 delete *iter;
2319 }
2320 if(m_clearTrack)delete trkRecoTrk;
2321 return err;
2322}
2323
2324int HoughTrack::fitCircle(DotsHelixFitter* fitter, double bunchT0)
2325{
2326 bool debug=false; //debug=true;
2327 // --- fill vector<MdcDigi*> and vector<RecCgemCluster*> from m_vecHitPnt
2328 ///*
2329 m_vecMdcDigiPnt.clear();
2330 m_vecCgemClusterPnt.clear();
2331 for(vector<HoughHit*>::iterator hotIt = m_vecHitPnt.begin(); hotIt != m_vecHitPnt.end(); hotIt++)
2332 {
2333 if((*hotIt)->getFlag()!=0)continue;
2334 if((*hotIt)->getHitType()==HoughHit::MDCHIT)
2335 {
2336 const MdcDigi* aMdcDigi = (*hotIt)->getDigi();
2337 m_vecMdcDigiPnt.push_back(aMdcDigi);
2338 }
2339 else if((*hotIt)->getHitType()==HoughHit::CGEMHIT)
2340 {
2341 const RecCgemCluster* aRecCgemCluster = (*hotIt)->getCgemCluster();
2342 m_vecCgemClusterPnt.push_back(aRecCgemCluster);
2343 }
2344 }
2345 //*/
2346
2347 // --- set
2348 HepPoint3D pivot(0,0,0);
2349 HepVector a(5,0);
2350 HepSymMatrix Ea(5,0);
2351
2352 // --- initial circle parameters from Hough transform
2353 double phi0_Hough = (m_rho>0)? m_angle:m_angle+M_PI;
2354 phi0_Hough = (m_charge<0)? phi0_Hough:phi0_Hough+M_PI;
2355 if(phi0_Hough>2*M_PI)phi0_Hough-=2*M_PI;
2356 if(phi0_Hough<0) phi0_Hough+=2*M_PI;
2357 double kappa_Hough = fabs(m_alpha*m_rho)*m_charge;
2358 a(2) = phi0_Hough;
2359 a(3) = kappa_Hough;
2360
2361 // --- from fit?
2362 //a(1) = m_dr;
2363 //a(2) = m_phi0;
2364 //a(3) = m_kappa;
2365 //a(5) = 0.00001;
2366
2367 KalmanFit::Helix ini_helix(pivot,a,Ea);
2368 if(debug) cout<<"circle ini_helix: "<<ini_helix.a()<<endl;
2369 fitter->setInitialHelix(ini_helix);
2370 fitter->setDChits(m_vecMdcDigiPnt, bunchT0);
2371 fitter->setCgemClusters(m_vecCgemClusterPnt);
2372 if(debug) cout<<"finish DotsHelixFitter setting"<<endl;
2373
2374 // --- circle fit
2375 //fitter->fitCircleOnly();
2376 fitter->fitModeCircle();
2377 int fitFlag = 0;
2378 while(1)
2379 {
2380 fitFlag = fitter->calculateNewHelix();
2381 if(fitFlag==0)
2382 {
2383 HepVector a_new = fitter->getHelix();
2384 if(debug) cout<<"circle helix from DotsHelixFitter: "<<a_new<<", chi2="<<fitter->getChi2()<<" "<<fitter->getNActiveHits()<<" hits kept"<<endl;
2385 if( fitter->getChi2() > 25*(fitter->getNActiveHits()) && fitter->deactiveHits(10.0))
2386 {
2387 continue;
2388 if(debug) cout<<"too large chi2, drop one worst DC hit!"<<endl;
2389 }
2390 else break;
2391 }
2392 else
2393 {
2394 if(debug) cout<<"DotsHelixFitter failed, fitFlag = "<<fitFlag<<endl;
2395 break;
2396 }
2397 }
2398 return fitFlag;
2399
2400}
2401
2402TrkErrCode HoughTrack::fitHelix(const MdcDetector* mdcDetector, TrkContextEv* trkContextEv, double bunchT0, vector<MdcHit*>& mdcHitCol, vector<HoughHit*>& hot)
2403{
2404 int nHot(0);
2405 double d0 = -m_dr;
2406 double phi0 = m_phi0+M_PI/2;
2407 double omega = m_kappa/fabs(alpha());
2408 double z0 = m_dz;
2409 double tanl = m_tanl;
2410 while(phi0>M_PI)phi0-=2*M_PI;
2411 while(phi0<-M_PI)phi0+=2*M_PI;
2412
2413 TrkExchangePar helixTrk(d0,phi0,omega,z0,tanl);
2414 TrkHelixMaker helixfactory;
2415 double chisum =0.;
2416 m_trkRecoTrk = helixfactory.makeTrack(helixTrk, chisum, *trkContextEv, bunchT0);
2417 bool permitFlips = true;
2418 bool lPickHits = true;
2419 helixfactory.setFlipAndDrop(*m_trkRecoTrk, permitFlips, lPickHits);
2420 TrkHitList* trkHitList = m_trkRecoTrk->hits();
2421
2422 //--- convert hits
2423 //cout<<maxLayer<<endl;
2424 for(vector<HoughHit*>::iterator hitIt = hot.begin(); hitIt != hot.end();hitIt++){
2425 //(*hitIt)->print();
2426 //cout<<endl;
2427 if((*hitIt)->getHitType()==HoughHit::MDCHIT /*|| (*hitIt)->getHitType()==HoughHit::MDCMCHIT*/){
2428 bool sameHit(false);
2429 const TrkHitList* aList = m_trkRecoTrk->hits();
2430 for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2431 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack)){
2432 const MdcRecoHitOnTrack* recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2433 int layerId = recoHot->mdcHit()->layernumber();
2434 int wireId = recoHot->mdcHit()->wirenumber();
2435 if((*hitIt)->getLayer()==layerId && (*hitIt)->getWire()==wireId)sameHit = true;
2436 }
2437 }
2438 if(sameHit)continue;
2439
2440 //const MdcDigi* mdcDigi = (*hitIt)->getDigi();
2441 //MdcHit *mdcHit = new MdcHit(mdcDigi,mdcDetector);
2442 MdcHit *mdcHit;
2443 vector<MdcHit*>::iterator imdcHit = mdcHitCol.begin();
2444 for(;imdcHit != mdcHitCol.end();imdcHit++){
2445 if((*imdcHit)->mdcId() == (*hitIt)->getDigi()->identify()){
2446 mdcHit = *imdcHit;
2447 break;
2448 }
2449 }
2450
2451 mdcHit->setCalibSvc((*hitIt)->getMdcCalibFunSvc());
2452 mdcHit->setCountPropTime(true);
2453 int newAmbig = 0;
2454 MdcRecoHitOnTrack mdcRecoHitOnTrack(*mdcHit, newAmbig, bunchT0*1.e9);
2455 MdcHitOnTrack* mdcHitOnTrack = &mdcRecoHitOnTrack;
2456 HepPoint3D midPoint = (*hitIt)->getHitPosition();
2457 //double fltLength = flightLength(midPoint);
2458 double fltLength = flightArc(midPoint);
2459 //double fltLength = flightLength((*hitIt)->getHitPosition());
2460 mdcHitOnTrack->setFltLen(fltLength);
2461 trkHitList->appendHot(mdcHitOnTrack);
2462 nHot++;
2463 }else if((*hitIt)->getHitType()==HoughHit::CGEMHIT/* || (*hitIt)->getHitType()==HoughHit::CGEMMCHIT*/){
2464 if(m_useCgemInGlobalFit<2)continue;
2465 bool sameHit(false);
2466 const TrkHitList* aList = m_trkRecoTrk->hits();
2467 for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2468 if(typeid(*hot)==typeid(CgemHitOnTrack)){
2469 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2470 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2471 int clusterid = recCgemCluster->getclusterid();
2472 if((*hitIt)->getCgemCluster()->getclusterid()==clusterid)sameHit = true;
2473 }
2474 }
2475 if(sameHit)continue;
2476
2477 CgemHitOnTrack* cgemHitOnTrack = new CgemHitOnTrack(*((*hitIt)->getCgemCluster()),bunchT0*1.e9);
2478 m_cgemHitVector.push_back(cgemHitOnTrack);
2479 if((*hitIt)->getFlag()==0)continue;
2480 cgemHitOnTrack->setCgemGeomSvc((*hitIt)->getCgemGeomSvc());
2481 cgemHitOnTrack->setCgemCalibFunSvc((*hitIt)->getCgemCalibFunSvc());
2482 trkHitList->appendHot(cgemHitOnTrack);
2483 nHot++;
2484 }
2485 }
2486 //cout<<trkHitList->nHit()<<endl;
2487 //---fitting
2488 TrkHitList* qhits = m_trkRecoTrk->hits();
2489 TrkErrCode err=qhits->fit();
2490 //err.print(std::cout);
2491 //cout<<endl;
2492 //cout<<endl;
2493
2494 //const TrkHitList* aList = m_trkRecoTrk->hits();
2495 //for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2496 //if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack)){
2497 //const MdcRecoHitOnTrack* recoHot;
2498 //recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2499 //int layerId = recoHot->mdcHit()->layernumber();
2500 //int wireId = recoHot->mdcHit()->wirenumber();
2501 //double res=999.,rese=999.;
2502 //if (recoHot->resid(res,rese,false)){
2503 //}else{}
2504 //std::cout<<"("<<layerId<<","<<wireId<<") res:"<<res<<std::endl;
2505 //}
2506 //else if(typeid(*hot)==typeid(CgemHitOnTrack)){
2507 //const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2508 //const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2509 //int clusterid = recCgemCluster->getclusterid();
2510 //cout<<"clusterid:"<< recCgemCluster->getclusterid()<<" layer:"<<recCgemCluster->getlayerid()<<endl;
2511 //}
2512 //}
2513
2514 if(err.success()){
2515 const TrkFit* fitResult = m_trkRecoTrk->fitResult();
2516 TrkExchangePar par = fitResult->helix(0.);
2517 m_dr = -par.d0();
2518 m_phi0 = par.phi0()-M_PI/2;
2519 m_kappa = par.omega()*fabs(alpha());
2520 m_dz = par.z0();
2521 m_tanl = par.tanDip();
2522 while(m_phi0>2*M_PI)m_phi0-=2*M_PI;
2523 while(m_phi0<0)m_phi0+=2*M_PI;
2524 //HepVector hepVector(5,0);
2525 //a(hepVector);
2526 updateHelix();
2527
2528 m_chi2 = fitResult->chisq();
2529 //cout<<"fitting params:"<<endl;
2530 //cout
2531 //<<setw(15)<<par.d0()
2532 //<<setw(15)<<par.phi0()
2533 //<<setw(15)<<par.omega()
2534 //<<setw(15)<<par.z0()
2535 //<<setw(15)<<par.tanDip()
2536 //;
2537 //cout<<endl;
2538 }
2539 //cout<<"getId()="<<m_trkRecoTrk->id()<<endl;
2540 return err;
2541}
2542
2543TrkErrCode HoughTrack::fitHelix(const MdcDetector* mdcDetector, const BField* bfield, double bunchT0, vector<HoughHit*>& hot,int maxLayer)
2544{
2545 // --- make TrkRecoTrk and trkHitList
2546 int nHot(0);
2547 double d0 = -m_dr;
2548 double phi0 = m_phi0+M_PI/2;
2549 double omega = m_kappa/fabs(alpha());
2550 double z0 = m_dz;
2551 double tanl = m_tanl;
2552 while(phi0>M_PI)phi0-=2*M_PI;
2553 while(phi0<-M_PI)phi0+=2*M_PI;
2554
2555 TrkExchangePar helixTrk(d0,phi0,omega,z0,tanl);
2556 TrkHelixMaker helixfactory;
2557 double chisum =0.;
2558 long idnum(0);
2559 TrkRecoTrk* trkRecoTrk = helixfactory.makeTrack(helixTrk, chisum, bfield, bunchT0);
2560 bool permitFlips = true;
2561 bool lPickHits = true;
2562 helixfactory.setFlipAndDrop(*trkRecoTrk, permitFlips, lPickHits);
2563 TrkHitList* trkHitList = trkRecoTrk->hits();
2564
2565 //--- convert hits and put them into trkHitList
2566 //cout<<maxLayer<<endl;
2567 vector<MdcHit*> mdcHitVector;
2568 mdcHitVector.clear();
2569 vector<CgemHitOnTrack*> cgemHitVector;
2570 cgemHitVector.clear();
2571
2572 for(vector<HoughHit*>::iterator hitIt = hot.begin(); hitIt != hot.end();hitIt++){
2573 //(*hitIt)->print();
2574 //cout<<endl;
2575 if((*hitIt)->getHitType()==HoughHit::MDCHIT /*|| (*hitIt)->getHitType()==HoughHit::MDCMCHIT*/){
2576 if((*hitIt)->getFlag()!=0&&(*hitIt)->getLayer()>maxLayer)continue;
2577 bool sameHit(false);
2578 const TrkHitList* aList = trkRecoTrk->hits();
2579 for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2580 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack)){
2581 const MdcRecoHitOnTrack* recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2582 int layerId = recoHot->mdcHit()->layernumber();
2583 int wireId = recoHot->mdcHit()->wirenumber();
2584 if((*hitIt)->getLayer()==layerId && (*hitIt)->getWire()==wireId)sameHit = true;
2585 }
2586 }
2587 if(sameHit)continue;
2588
2589 const MdcDigi* mdcDigi = (*hitIt)->getDigi();
2590 MdcHit *mdcHit = new MdcHit(mdcDigi,mdcDetector);
2591 mdcHitVector.push_back(mdcHit);
2592
2593 mdcHit->setCalibSvc((*hitIt)->getMdcCalibFunSvc());
2594 mdcHit->setCountPropTime(true);
2595 int newAmbig = 0;
2596 MdcRecoHitOnTrack mdcRecoHitOnTrack(*mdcHit, newAmbig, bunchT0*1.e9);
2597 MdcHitOnTrack* mdcHitOnTrack = &mdcRecoHitOnTrack;
2598 HepPoint3D midPoint = (*hitIt)->getHitPosition();
2599 //double fltLength = flightLength(midPoint);
2600 double fltLength = flightArc(midPoint);
2601 //double fltLength = flightLength((*hitIt)->getHitPosition());
2602 mdcHitOnTrack->setFltLen(fltLength);
2603 trkHitList->appendHot(mdcHitOnTrack);
2604 nHot++;
2605 }else if((*hitIt)->getHitType()==HoughHit::CGEMHIT/* || (*hitIt)->getHitType()==HoughHit::CGEMMCHIT*/){
2606 if(m_useCgemInGlobalFit<3)continue;
2607 bool sameHit(false);
2608 const TrkHitList* aList = trkRecoTrk->hits();
2609 for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2610 if(typeid(*hot)==typeid(CgemHitOnTrack)){
2611 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2612 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2613 int clusterid = recCgemCluster->getclusterid();
2614 if((*hitIt)->getCgemCluster()->getclusterid()==clusterid)sameHit = true;
2615 }
2616 }
2617 if(sameHit)continue;
2618
2619 CgemHitOnTrack* cgemHitOnTrack = new CgemHitOnTrack(*((*hitIt)->getCgemCluster()),bunchT0*1.e9);
2620 cgemHitVector.push_back(cgemHitOnTrack);
2621 if((*hitIt)->getFlag()==0)continue;
2622 cgemHitOnTrack->setCgemGeomSvc((*hitIt)->getCgemGeomSvc());
2623 cgemHitOnTrack->setCgemCalibFunSvc((*hitIt)->getCgemCalibFunSvc());
2624 trkHitList->appendHot(cgemHitOnTrack);
2625 nHot++;
2626 }
2627 }
2628 //cout<<trkHitList->nHit()<<endl;
2629 //---fitting
2630 TrkHitList* qhits = trkRecoTrk->hits();
2631 TrkErrCode err=qhits->fit();
2632 //err.print(std::cout);cout<<endl;cout<<endl;
2633
2634 //const TrkHitList* aList = trkRecoTrk->hits();
2635 //for (TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2636 //if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack)){
2637 //const MdcRecoHitOnTrack* recoHot;
2638 //recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2639 //int layerId = recoHot->mdcHit()->layernumber();
2640 //int wireId = recoHot->mdcHit()->wirenumber();
2641 //double res=999.,rese=999.;
2642 //if (recoHot->resid(res,rese,false)){
2643 //}else{}
2644 //std::cout<<"("<<layerId<<","<<wireId<<") res:"<<res<<std::endl;
2645 //}
2646 //else if(typeid(*hot)==typeid(CgemHitOnTrack)){
2647 //const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2648 //const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2649 //int clusterid = recCgemCluster->getclusterid();
2650 //cout<<"clusterid:"<< recCgemCluster->getclusterid()<<" layer:"<<recCgemCluster->getlayerid()<<endl;
2651 //}
2652 //}
2653
2654 // --- if fit successes, update the helix parameters
2655 if(err.success()){
2656 const TrkFit* fitResult = trkRecoTrk->fitResult();
2657 TrkExchangePar par = fitResult->helix(0.);
2658 m_dr = -par.d0();
2659 m_phi0 = par.phi0()-M_PI/2;
2660 m_kappa = par.omega()*fabs(alpha());
2661 m_dz = par.z0();
2662 m_tanl = par.tanDip();
2663 while(m_phi0>2*M_PI)m_phi0-=2*M_PI;
2664 while(m_phi0<0)m_phi0+=2*M_PI;
2665 updateHelix();
2666 m_chi2 = fitResult->chisq();
2667
2668 int nHotKept = 0;
2669 for(TrkHitList::hot_iterator hot = qhits->begin();hot!=qhits->end();hot++)
2670 {
2671 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ )
2672 {
2673 const MdcRecoHitOnTrack* recoHot;
2674 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2675 if(!(recoHot->isActive())) continue;
2676 nHotKept++;
2677 //const MdcDigi* aMdcDigi = recoHot->mdcHit()->digi();
2678 //m_vecMdcDigiPnt.push_back(aMdcDigi);
2679 }
2680 else if(typeid(*hot)==typeid(CgemHitOnTrack))
2681 {
2682 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2683 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2684 int clusterid = recCgemCluster->getclusterid();
2685 int stat = recoHot->isActive();
2686 if(stat==0) continue;//
2687 nHotKept++;
2688 //m_vecCgemClusterPnt.push_back(recCgemCluster);
2689 }
2690 }
2691
2692 cout<<"chi2= "<<m_chi2
2693 <<", "<<nHotKept<<" hits kept after TrkFitter"
2694 <<", helix from TrkFitter = "<<a()
2695 <<endl;
2696 }
2697 for(vector<MdcHit*>::iterator iter = mdcHitVector.begin(); iter != mdcHitVector.end(); iter++){
2698 delete *iter;
2699 }
2700 for(vector<CgemHitOnTrack*>::iterator iter = cgemHitVector.begin(); iter != cgemHitVector.end(); iter++){
2701 delete *iter;
2702 }
2703 if(m_clearTrack)delete trkRecoTrk;
2704 //cout<<"getId():"<<trkRecoTrk->id()<<endl;
2705 return err;
2706}
2707
2708int HoughTrack::fitHelix(DotsHelixFitter* fitter, double bunchT0, RecCgemClusterCol::iterator recCgemClusterColBegin, double averageChi2cut)
2709{
2710 bool debug=false; //debug=true;
2711 // --- fill vector<MdcDigi*> and vector<RecCgemCluster*> from m_vecHitPnt
2712 //cout<<"before fill"<<endl;
2713 m_vecMdcDigiPnt.clear();
2714 m_vecCgemClusterPnt.clear();
2715 for(vector<HoughHit*>::iterator hotIt = m_vecHitPnt.begin(); hotIt != m_vecHitPnt.end(); hotIt++)
2716 {
2717 //if((*hotIt)->getFlag()!=0)continue;
2718 if((*hotIt)->getHitType()==HoughHit::MDCHIT)
2719 {
2720 const MdcDigi* aMdcDigi = (*hotIt)->getDigi();
2721 m_vecMdcDigiPnt.push_back(aMdcDigi);
2722 }
2723 else if((*hotIt)->getHitType()==HoughHit::CGEMHIT)
2724 {
2725 const RecCgemCluster* aRecCgemCluster = (*hotIt)->getCgemCluster();
2726 m_vecCgemClusterPnt.push_back(aRecCgemCluster);
2727 }
2728 }
2729 for(vector<HoughHit*>::iterator hotIt = m_vecStereoHitPnt.begin(); hotIt != m_vecStereoHitPnt.end(); hotIt++)
2730 {
2731 if((*hotIt)->getHitType()==HoughHit::MDCHIT)
2732 {
2733 const MdcDigi* aMdcDigi = (*hotIt)->getDigi();
2734 m_vecMdcDigiPnt.push_back(aMdcDigi);
2735 }
2736 else if((*hotIt)->getHitType()==HoughHit::CGEMHIT)
2737 {
2738 const RecCgemCluster* aRecCgemCluster = (*hotIt)->getCgemCluster();
2739 if(aRecCgemCluster->getflag()==2)
2740 {
2741 int clusterId = aRecCgemCluster->getclusterflage();
2742 RecCgemClusterCol::iterator cgemClusterIter = recCgemClusterColBegin+clusterId;
2743 //cout<<"V cluster ID = "<<setw(10)<<clusterId<<setw(10)<<(*cgemClusterIter)->getclusterid();
2744 if(clusterId==(*cgemClusterIter)->getclusterid())
2745 m_vecCgemClusterPnt.push_back(*cgemClusterIter);
2746 }
2747 }
2748 }
2749 //cout<<"end filling"<<endl;
2750
2751
2752 // --- set
2753 HepPoint3D pivot(0,0,0);
2754 HepVector a(5,0);
2755 HepSymMatrix aEa(5,0);
2756 a(1) = m_dr;
2757 a(2) = m_phi0;
2758 a(3) = m_kappa;
2759 a(4) = m_dz;
2760 a(5) = m_tanl;
2761 KalmanFit::Helix ini_helix(pivot,a,aEa);
2762 if(debug) cout<<"fitHelix ini_helix: "<<ini_helix.a()<<endl;
2763 fitter->setInitialHelix(ini_helix);
2764 fitter->setDChits(m_vecMdcDigiPnt, bunchT0);
2765 fitter->setCgemClusters(m_vecCgemClusterPnt);
2766 int nMdcDigi = m_vecMdcDigiPnt.size();
2767
2768 // --- helix fit
2769 fitter->fitModeHelix();
2770 int fitFlag = 0;
2771 while(1)
2772 {
2773 fitFlag = fitter->calculateNewHelix();
2774 if(fitFlag==0)
2775 {
2776 HepVector a_new = fitter->getHelix();
2777 if(debug) {
2778 cout<<"helix from DotsHelixFitter: "<<a_new<<", chi2="<<fitter->getChi2()<<" "<<fitter->getNActiveHits()<<" hits kept"<<endl;
2779 cout<<"nMdcDigi = "<<nMdcDigi<<endl;
2780 }
2781 //if( fitter->getChi2() > averageChi2cut*(fitter->getNActiveHits()) && fitter->deactiveHits())
2782 if( fitter->getChi2() > averageChi2cut*(fitter->getNActiveHits()) && fitter->deactiveHits(43,--nMdcDigi))
2783 {
2784 //cout<<"too large chi2, drop the worst DC hit!"<<endl;
2785 if(debug) cout<<"too large chi2, drop the outmost DC hit!"<<endl;
2786 continue;
2787 }
2788 else
2789 {
2790 if(debug) cout<<"stop droping hit! "<<endl;
2791 break;
2792 }
2793 }
2794 else
2795 {
2796 if(debug) cout<<"DotsHelixFitter failed, fitFlag = "<<fitFlag<<endl;
2797 break;
2798 }
2799 }
2800 if(fitFlag==0) //update dr ... and Ea
2801 {
2802 HepVector a = fitter->getHelix();
2803 m_dr = a[0];
2804 m_phi0 = a[1];
2805 m_kappa = a[2];
2806 m_dz = a[3];
2807 m_tanl = a[4];
2808 m_chi2 = fitter->getChi2();
2809 Ea(fitter->getEa());
2810 }
2811 return fitFlag;
2812}
2813
2814void HoughTrack::clearMemory()
2815{
2816 for(vector<CgemHitOnTrack*>::iterator iter = m_cgemHitVector.begin(); iter != m_cgemHitVector.end(); iter++){
2817 delete *iter;
2818 }
2819 if(m_clearTrack)delete m_trkRecoTrk;
2820}
2821//cout<<__FILE__<<" "<<__LINE__<<endl;
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
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:101
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: Gam4pikp.cxx:37
s
XmlRpcServer s
Definition: HelloServer.cpp:11
sortByLayer
bool sortByLayer(HoughHit *hit1, HoughHit *hit2)
Definition: HoughTrack.cxx:913
innerLayer
bool innerLayer(HoughHit hit1, HoughHit hit2)
Definition: HoughTrack.cxx:903
cut
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared cut
Definition: KarFin.h:27
m_dz
NTuple::Array< double > m_dz
Definition: MdcHistItem.h:117
m_phi0
NTuple::Item< double > m_phi0
Definition: MdcHistItem.h:80
m_chi2
NTuple::Item< double > m_chi2
Definition: MdcHistItem.h:91
m_tanl
NTuple::Item< double > m_tanl
Definition: MdcHistItem.h:83
M_PI
#define M_PI
Definition: TConstant.h:4
CgemHitOnTrack::setCgemCalibFunSvc
void setCgemCalibFunSvc(const CgemCalibFunSvc *svc)
Definition: CgemHitOnTrack.h:115
CgemHitOnTrack::baseHit
const RecCgemCluster * baseHit() const
Definition: CgemHitOnTrack.h:119
CgemHitOnTrack::setCgemGeomSvc
void setCgemGeomSvc(const CgemGeomSvc *svc)
Definition: CgemHitOnTrack.h:114
DotsHelixFitter::setInitialHelix
void setInitialHelix(KalmanFit::Helix aHelix)
Definition: DotsHelixFitter.cxx:702
DotsHelixFitter::deactiveHits
int deactiveHits(double chi_cut=10, int nMax=1)
Definition: DotsHelixFitter.cxx:600
DotsHelixFitter::setDChits
void setDChits(vector< const MdcDigi * > aVecMdcDigi, double T0)
Definition: DotsHelixFitter.cxx:139
DotsHelixFitter::setCgemClusters
void setCgemClusters(vector< const RecCgemCluster * > aVecCgemCluster)
Definition: DotsHelixFitter.cxx:159
DotsHelixFitter::getEa
const HepSymMatrix & getEa()
Definition: DotsHelixFitter.h:44
DotsHelixFitter::getHelix
HepVector getHelix()
Definition: DotsHelixFitter.h:43
Helix::flightLength
double flightLength(HepPoint3D &hit) const
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:812
Helix::center
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:212
Helix::set
void set(const HepPoint3D &pivot, const HepVector &a, const HepSymMatrix &Ea)
sets helix pivot position, parameters, and error matrix.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:392
Helix::operator=
Helix & operator=(const Helix &)
Copy operator.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:403
Helix::Ea
const HepSymMatrix & Ea(void) const
returns error matrix.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:278
Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:152
Helix::radius
double radius(void) const
returns radious of helix.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:224
Helix::dr
double dr(void) const
returns an element of parameters.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:236
Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:272
Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:218
HoughHit::getMdcGeomSvc
MdcGeomSvc * getMdcGeomSvc() const
Definition: HoughHit.h:65
HoughHit::getDigi
const MdcDigi * getDigi() const
Definition: HoughHit.h:42
HoughHit::getCgemCluster
const RecCgemCluster * getCgemCluster() const
Definition: HoughHit.h:41
HoughHit::getDriftDist
double getDriftDist() const
Definition: HoughHit.h:54
HoughHit::getMdcCalibFunSvc
MdcCalibFunSvc * getMdcCalibFunSvc() const
Definition: HoughHit.h:66
HoughHit::addSZ
void addSZ(S_Z sz)
Definition: HoughHit.h:106
HoughHit::getHitPosition
HepPoint3D getHitPosition() const
Definition: HoughHit.h:51
HoughHit::getHitType
HitType getHitType() const
Definition: HoughHit.h:40
HoughHit::getCgemGeomSvc
CgemGeomSvc * getCgemGeomSvc() const
Definition: HoughHit.h:67
HoughHit::getFlag
int getFlag() const
Definition: HoughHit.h:47
HoughHit::print
void print()
Definition: HoughHit.cxx:327
HoughHit::getCgemCalibFunSvc
CgemCalibFunSvc * getCgemCalibFunSvc() const
Definition: HoughHit.h:68
HoughHit::getLayer
int getLayer() const
Definition: HoughHit.h:45
HoughHit::getWire
int getWire() const
Definition: HoughHit.h:46
HoughHit::resetSZ
void resetSZ()
Definition: HoughHit.h:97
HoughHit::CGEMMCHIT
@ CGEMMCHIT
Definition: HoughHit.h:27
HoughHit::CGEMHIT
@ CGEMHIT
Definition: HoughHit.h:27
HoughHit::MDCHIT
@ MDCHIT
Definition: HoughHit.h:27
HoughHit::MDCMCHIT
@ MDCMCHIT
Definition: HoughHit.h:27
HoughHit::setResidual
void setResidual(double residual)
Definition: HoughHit.h:90
HoughTrack::clearHits
void clearHits()
Definition: HoughTrack.cxx:1193
HoughTrack::print
void print()
Definition: HoughTrack.cxx:962
HoughTrack::updateCirclePar
void updateCirclePar(double dr, double phi0, double kappa)
Definition: HoughTrack.cxx:953
HoughTrack::findXHot
int findXHot(vector< HoughHit * > &hitList, int charge)
Definition: HoughTrack.cxx:481
HoughTrack::clearMemory
void clearMemory()
Definition: HoughTrack.cxx:2814
HoughTrack::getChi2
double getChi2() const
Definition: HoughTrack.h:58
HoughTrack::m_useCgemInGlobalFit
static int m_useCgemInGlobalFit
Definition: HoughTrack.h:146
HoughTrack::isAGoodCircle
bool isAGoodCircle()
Definition: HoughTrack.cxx:1507
HoughTrack::driftDistRes
double driftDistRes(HoughHit *hit)
Definition: HoughTrack.cxx:636
HoughTrack::sortHot
void sortHot(vector< HoughHit * > &hotList)
Definition: HoughTrack.cxx:918
HoughTrack::m_clearTrack
static int m_clearTrack
Definition: HoughTrack.h:145
HoughTrack::findVHot
int findVHot(vector< HoughHit * > &hitList, int charge, int maxLayer, double cutFactor=1.0)
Definition: HoughTrack.cxx:1996
HoughTrack::markUsedHot
void markUsedHot(vector< HoughHit * > &hitPntList, int use=1)
Definition: HoughTrack.cxx:1467
HoughTrack::dropHitPnt
void dropHitPnt(HoughHit *aHitPnt)
Definition: HoughTrack.cxx:1861
HoughTrack::getNHitsShared
int getNHitsShared()
Definition: HoughTrack.cxx:1983
HoughTrack::fitCircle
TrkErrCode fitCircle(const MdcDetector *mdcDetector, TrkContextEv *trkContextEv, double bunchT0)
Definition: HoughTrack.cxx:2179
HoughTrack::getTrkID
int getTrkID() const
Definition: HoughTrack.h:49
HoughTrack::printRecMdcHitVec
void printRecMdcHitVec()
Definition: HoughTrack.cxx:1119
HoughTrack::judgeHalf
int judgeHalf(HoughHit *hit)
Definition: HoughTrack.cxx:596
HoughTrack::dropRedundentCgemXHits
void dropRedundentCgemXHits()
Definition: HoughTrack.cxx:1886
HoughTrack::releaseSelHits
void releaseSelHits()
Definition: HoughTrack.cxx:1852
HoughTrack::getVecHitPnt
vector< HoughHit * > getVecHitPnt()
Definition: HoughTrack.h:83
HoughTrack::updateHelix
void updateHelix()
Definition: HoughTrack.cxx:923
HoughTrack::update
void update(double angle, double rho)
Definition: HoughTrack.cxx:935
HoughTrack::operator=
HoughTrack & operator=(const HoughTrack &other)
Definition: HoughTrack.cxx:255
HoughTrack::getVecStereoHitPnt
vector< HoughHit * > getVecStereoHitPnt()
Definition: HoughTrack.h:134
HoughTrack::judgeCharge
int judgeCharge(HoughHit *hit)
Definition: HoughTrack.cxx:609
HoughTrack::getHotList
vector< HoughHit * > getHotList(int type=2)
Definition: HoughTrack.cxx:2127
HoughTrack::dropVHitPnt
void dropVHitPnt(HoughHit *aHitPnt)
Definition: HoughTrack.cxx:1876
HoughTrack::m_cut
static TGraph * m_cut[2][43]
Definition: HoughTrack.h:147
HoughTrack::calculateZ_S
int calculateZ_S(HoughHit *hit)
Definition: HoughTrack.cxx:774
HoughTrack::dropRedundentCgemVHits
void dropRedundentCgemVHits()
Definition: HoughTrack.cxx:1945
HoughTrack::printHot
void printHot()
Definition: HoughTrack.cxx:1020
HoughTrack::fitHelix
int fitHelix(DotsHelixFitter *fitter, double bunchT0, RecCgemClusterCol::iterator recCgemClusterColBegin, double averageChi2cut=25)
Definition: HoughTrack.cxx:2708
KalmanFit::Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-15-14/KalFitAlg/helix/Helix.h:265
MdcGeoWire::Forward
HepPoint3D Forward(void) const
Definition: MdcGeoWire.h:129
MdcGeoWire::Backward
HepPoint3D Backward(void) const
Definition: MdcGeoWire.h:128
MdcGeomSvc::Wire
const MdcGeoWire *const Wire(unsigned id)
Definition: MdcGeomSvc.cxx:770
MdcHit
Definition: MdcHit.h:44
MdcHit::setCalibSvc
void setCalibSvc(const MdcCalibFunSvc *calibSvc)
Definition: MdcHit.cxx:136
MdcHit::layernumber
unsigned layernumber() const
Definition: MdcHit.h:61
MdcHit::wirenumber
unsigned wirenumber() const
Definition: MdcHit.h:62
MdcHit::setCountPropTime
void setCountPropTime(const bool count)
Definition: MdcHit.h:86
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition: MdcID.cxx:49
MdcID::wire
static int wire(const Identifier &id)
Definition: MdcID.cxx:54
MdcRecoHitOnTrack::mdcHit
const MdcHit * mdcHit() const
Definition: MdcRecoHitOnTrack.cxx:38
RecCgemCluster::getclusterid
int getclusterid(void) const
Definition: RecCgemCluster.h:53
RecCgemCluster::getflag
int getflag(void) const
Definition: RecCgemCluster.h:57
RecCgemCluster::getclusterflagb
int getclusterflagb(void) const
Definition: RecCgemCluster.h:71
RecCgemCluster::getclusterflage
int getclusterflage(void) const
Definition: RecCgemCluster.h:72
TrkAbsFit::chisq
virtual double chisq() const =0
TrkErrCode::success
int success() const
Definition: TrkErrCode.h:62
TrkExchangePar::phi0
double phi0() const
Definition: TrkExchangePar.h:43
TrkExchangePar::omega
double omega() const
Definition: TrkExchangePar.h:44
TrkExchangePar::z0
double z0() const
Definition: TrkExchangePar.h:45
TrkExchangePar::d0
double d0() const
Definition: TrkExchangePar.h:42
TrkExchangePar::tanDip
double tanDip() const
Definition: TrkExchangePar.h:46
TrkFit
Definition: TrkFit.h:23
TrkFit::helix
virtual TrkExchangePar helix(double fltL) const =0
TrkFundHit::begin
TrkFundHit::hot_iterator begin() const
Definition: TrkFundHit.h:113
TrkHitList::fit
TrkErrCode fit()
Definition: TrkHitList.cxx:59
TrkHitList::begin
hot_iterator begin() const
Definition: TrkHitList.h:45
TrkHitList::appendHot
TrkHitOnTrk * appendHot(const TrkHitOnTrk *theHot)
Definition: TrkHitList.cxx:91
TrkHitList::end
hot_iterator end() const
Definition: TrkHitList.h:46
TrkHitOnTrk::setFltLen
void setFltLen(double f)
Definition: TrkHitOnTrk.h:147
TrkHitOnTrk::isActive
bool isActive() const
Definition: TrkHitOnTrk.h:200
TrkRecoTrk::hits
TrkHitList * hits()
Definition: TrkRecoTrk.h:107
TrkRecoTrk::fitResult
const TrkFit * fitResult() const
Definition: TrkRecoTrk.cxx:387
TrkSimpleMaker::setFlipAndDrop
bool setFlipAndDrop(TrkRecoTrk &, bool allowFlips, bool allowDrops) const
TrkSimpleMaker::makeTrack
TrkRecoTrk * makeTrack(const TrkExchangePar &helix, const double chi2, const TrkContext &, double trackT0) const
c1
TCanvas * c1
Definition: tau_mode.c:75