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