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