CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcHoughFinder/MdcHoughFinder-00-00-12/src/HoughTrack.cxx
Go to the documentation of this file.
1#include "MdcHoughFinder/HoughTrack.h"
2#include "MdcHoughFinder/HoughMap.h"
3#include "MdcHoughFinder/HoughStereo.h"
4//#include "MdcData/CgemHitOnTrack.h"
5#include "CgemData/CgemHitOnTrack.h"
6#include <vector>
7#include <map>
8#include "TF1.h"
9#include "TGraph.h"
10using namespace std;
11
12bool HoughTrack::m_debug=0;
13bool HoughTrack::m_globalfit=0;
14
15HoughTrack::HoughTrack(){
16 p_trk=NULL;
17 p_trk2D=NULL;
18}
19HoughTrack::~HoughTrack(){
20}
21
22HoughTrack& HoughTrack::operator=(const HoughTrack& other){
23 _dist=(other._dist);
24 _charge=(other._charge);
25 _ptLeast=(other._ptLeast);
26 _pt2D=(other._pt2D);
27 _pt3D=(other._pt3D);
28 _pz=(other._pz);
29 _p=(other._p);
30
31 _d0=other._d0;
32 _omega=other._omega;
33 _phi0=other._phi0;
34 _tanl=(other._tanl);
35 _z0=(other._z0);
36 _tanl_zs=(other._tanl_zs);
37 _z0_zs=(other._z0_zs);
38
39 _centerPeak=(other._centerPeak);
40 _Hough2D=(other._Hough2D);
41 _Hough3D=(other._Hough3D);
42 _bunchTime=(other._bunchTime);
43 _centerX=(other._centerX);
44 _centerY=(other._centerY);
45 _centerR=(other._centerR);
46 _chi2_aver=(other._chi2_aver);
47 _nfit=(other._nfit);
48 _chi2_aver2D=(other._chi2_aver2D);
49 _nfit2D=(other._nfit2D);
50 _recHitVec=(other._recHitVec);
51 _stat2D=(other._stat2D);
52 _stat3D=(other._stat3D);
53 p_trk=other.p_trk;
54 p_trk2D=other.p_trk2D;
55 _maprho = other._maprho;
56 _maptheta= other._maptheta;
57 t_pro_correct= other.t_pro_correct;
58 _houghList = other._houghList;
59 vec_mdcHit = other.vec_mdcHit;
60 _trkid = other._trkid;
61 // if(_recHitVec.size() != 0 ){
62 // for(int i =0;i<_recHitVec.size();i++){
63 // delete _recHitVec[i];
64 // _recHitVec[i]=NULL;
65 // }
66 // _recHitVec.clear();
67 // }
68 // for(int i =0;i<other._recHitVec.size();i++){
69 // //const HoughHit* p_hit= new HoughHit( *(other._recHitVec[i]));
70 // HoughRecHit p_hit ( (trackHitList[i])->digi(),0.,0. );
71 // _recHitVec.push_back(p_hit);
72 // }
73}
74
75HoughTrack::HoughTrack(const HoughTrack& other):
76 _dist(other._dist),
77 _charge(other._charge),
78 _ptLeast(other._ptLeast),
79 _pt2D(other._pt2D),
80 _pt3D(other._pt3D),
81 _pz(other._pz),
82 _p(other._p),
83
84 _d0(other._d0),
85 _omega(other._omega),
86 _phi0(other._phi0),
87 _z0(other._z0),
88 _tanl(other._tanl),
89 _z0_zs(other._z0_zs),
90 _tanl_zs(other._tanl_zs),
91
92 _centerPeak(other._centerPeak),
93 _Hough2D(other._Hough2D),
94 _Hough3D(other._Hough3D),
95 _bunchTime(other._bunchTime),
96 _centerX(other._centerX),
97 _centerY(other._centerY),
98 _centerR(other._centerR),
99 _chi2_aver(other._chi2_aver),
100 _nfit(other._nfit),
101 _chi2_aver2D(other._chi2_aver2D),
102 _nfit2D(other._nfit2D),
103 _stat2D(other._stat2D),
104 _stat3D(other._stat3D),
105 _recHitVec(other._recHitVec),
106 p_trk(other.p_trk),
107 p_trk2D(other.p_trk2D),
108 _maprho(other._maprho),
109 _maptheta(other._maptheta),
110 t_pro_correct(other.t_pro_correct),
111 _houghList(other._houghList),
112 vec_mdcHit ( other.vec_mdcHit),
113 _trkid(other._trkid)
114
115{
116 // _recHitVec.clear();
117 //if(_recHitVec.size() != 0 ){
118 // //cout<<"clear first "<<endl;
119 // for(int i =0;i<_recHitVec.size();i++){
120 // delete _recHitVec[i];
121 // _recHitVec[i]=NULL;
122 // }
123 // _recHitVec.clear();
124 //}
125 //for(int i =0;i<other._recHitVec.size();i++){
126 // //const HoughHit* p_hit= new HoughHit( *(other._recHitVec[i]));
127 // HoughRecHit p_hit ( (trackHitList[i])->digi(),0.,0. );
128 // _recHitVec.push_back(p_hit);
129 //}
130}
131/*
132HoughTrack& HoughTrack::add(const HoughTrack& other)
133{
134 for(int i =0;i<other._recHitVec.size();i++){
135 int same=0;
136 for(int j =0;j<_recHitVec.size();j++){
137 if( _recHitVec[j].getDetectorType() == MDC && other._recHitVec[i].getDetectorType() == MDC && _recHitVec[j].digi() == other._recHitVec[i].digi() ) { same=1;}
138 if( _recHitVec[j].getDetectorType() == CGEM && other._recHitVec[i].getDetectorType() == CGEM && _recHitVec[j].getCluster() == other._recHitVec[i].getCluster() ) { same=1;}
139 }
140 if( same==0 ) {
141 //const HoughHit* p_hit= new HoughHit( *(other._recHitVec[i]));
142 //HoughRecHit p_hit ( (other._recHitVec[i]).digi(),0.,1);
143 HoughRecHit p_hit ( (other._recHitVec[i]));
144 _recHitVec.push_back(p_hit);
145 }
146 //cout<<" after copy: "<<i<<" new: "<<p_hit<<" old: "<<other._recHitVec[i]<<endl;
147 }
148 return *this;
149}
150*/
151HoughTrack::HoughTrack(const HoughPeak& centerPeak , std::vector<const HoughHit*> trackHitList,double rho,double theta, int trackId){
152 _centerPeak=centerPeak;
153 _recHitVec.clear();
154 _stat2D=0;
155 _stat3D=0;
156 _tanl=-999;
157 _z0=-999;
158 _ptLeast=-999;
159 _pt2D=-999;
160 _pt3D=-999;
161 _p=-999;
162 _d0=-999;
163 _omega=-999;
164 _phi0=-999;
165 _z0=-999;
166 _tanl=-999;
167 _chi2_aver=-999;
168 _nfit=-999;
169 p_trk=NULL;
170 p_trk2D=NULL;
171 _charge=0; //undeter mine
172 _maprho=rho;
173 _maptheta=theta;
174 _trkid = trackId;
175 //if(_recHitVec.size()!= 0 ){
176 // cout<<" Mind trackhitlist size !=0 "<<endl;
177 // for(int i =0;i<_recHitVec.size();i++){
178 // delete _recHitVec[i];
179 // _recHitVec[i]=NULL;
180 // }
181 // _recHitVec.clear();
182 //}
183 int t_size = trackHitList.size();
184 for(int i =0;i<t_size;i++){
185 HoughRecHit p_hit ( *(trackHitList[i]),0.,0.,1);
186 p_hit.setPtr2D(&_Hough2D);
187 p_hit.setflag(0);
188 _recHitVec.push_back(p_hit);
189 }
190}
191
192bool digi_in_track(const HoughRecHit& hita,const HoughRecHit& hitb){
193 return hita.getLayerId() < hitb.getLayerId();
194}
195void HoughTrack::sortHit(){
196 std::sort (_recHitVec.begin(),_recHitVec.end(),digi_in_track);
197}
198int HoughTrack::fit2D(double bunchtime){
199 sortHit();
200 // cout<<"in fit2d"<<endl;
201 _bunchTime=bunchtime;
202 _stat2D = fitLeast();
203 cald_layer();
204 // outerHit();
205 if(m_globalfit)_stat2D = fit_global2D(_recHitVec);
206 else _stat2D =1;
207 //_stat2D =1;
208 //printRecHit();
209 //collectAxialHit();
210 //_Hough2D.print();
211 //cald_layer();
212 return _stat2D;
213}
214
215bool less_layer(const int& a,const int& b){
216 return a < b;
217}
218void HoughTrack::outerHit(){
219 vector<int> vec_layer_num[43];
220 for(int ilay=0;ilay<43;ilay++){
221 for(int ihit=0;ihit<_recHitVec.size();ihit++){
222 if(_recHitVec[ihit].getDetectorType()==CGEM)continue;
223 if( _recHitVec[ihit].getLayerId()==ilay && _recHitVec[ihit].getflag()==0 ) {
224 vec_layer_num[ilay].push_back( _recHitVec[ihit].getWireId() );
225 }
226 }
227 std::sort( vec_layer_num[ilay].begin(),vec_layer_num[ilay].end(),less_layer);
228 //for(int j=0;j<vec_layer_num[ilay].size();j++) cout<<"("<<ilay<<","<<vec_layer_num[ilay][j]<<")"<<endl;
229 if( vec_layer_num[ilay].size() < 4 ) continue;
230 for(int j=0;j<vec_layer_num[ilay].size();j++) {
231 // if( (vec_layer_num[ilay][j]+1 == vec_layer_num[ilay][j+1]) && (vec_layer_num[ilay][j+1]+1 == vec_layer_num[ilay][j+2]) && (vec_layer_num[ilay][j+2]+1 == vec_layer_num[ilay][j+3]) ) {
232 for(int jhit=0;jhit<_recHitVec.size();jhit++) {
233 if(_recHitVec[jhit].getDetectorType()==CGEM)continue;
234 if( (ilay==_recHitVec[jhit].getLayerId()) && (vec_layer_num[ilay][j]==_recHitVec[jhit].getWireId() ) ) {_recHitVec[jhit].setflag(-1);}
235 }
236 //}
237 }
238 }
239}
240
241
242//int HoughTrack::collectAxialHit(){
243// for( int ilay=8;ilay<20;ilay++){
244// for( int ihit=0;ihit<_recHitVec.size();ihit++){
245// if( _recHitVec[ihit].getLayerId()==ilay && _recHitVec[ihit].getflag()==0) {
246// _recHitVec[ihit].setflag(1);
247// }
248// }
249// }
250//}
251/*
252int HoughTrack::collectAxialHit(){
253 vector<HoughRecHit> hitCol;
254 int lay_contain[12]; //FIXME
255 for( int ilay=8;ilay<20;ilay++){
256 int lay_count=0;
257 for( int ihit=0;ihit<_recHitVec.size();ihit++){
258 if( _recHitVec[ihit].getLayerId()==ilay && _recHitVec[ihit].getflag()==0) {
259 hitCol.push_back( _recHitVec[ihit]);
260 lay_count++;
261 }
262 }
263 lay_contain[ilay-8]=lay_count;
264 }
265 //int **a;
266 int **a;
267 int *len;
268 int p[12];
269 a = (int**)malloc(sizeof(int*) * 12);
270 len = (int*)malloc(sizeof(int) * 12);
271 for (int i = 0; i < 12; i++) {
272 int add_lay_count=0;
273 for (int s = 0; s < i; s++) {
274 add_lay_count+=lay_contain[s];
275 }
276 int n = lay_contain[i];
277 len[i] = n;
278 a[i] = (int*)malloc(sizeof(int) * n);
279 for (int j = 0; j < n; j++) {
280 a[i][j]= add_lay_count+j;
281 cout<<"aij "<<a[i][j]<<endl;
282 cout<<"("<<hitCol[add_lay_count+j].getLayerId()<<","<<hitCol[add_lay_count+j].getWireId()<<")"<<endl;
283 }
284 }
285 int numall=0;
286 fun(0,a,len,p,hitCol,numall);
287 cout<<" num all "<<numall<<endl;
288}
289void HoughTrack::fun(int level,int** a,int* len,int* p,vector<HoughRecHit>& hitCol,int& numall) {
290 if (level < 12) {
291 cout<<"level "<<level<<" len "<<len[level]<<endl;
292 for (int i = 0; i < len[level]; i++) {
293 cout<<i<<" i "<<endl;
294 p[level] = i;
295 cout<<i<<" i "<<endl;
296 cout<<"fun "<<level+1<<endl;
297 fun(level+1,a,len,p,hitCol,numall);
298 }
299 } else {
300 vector<HoughRecHit> recHitCol;
301 for (int i = 0; i < 12; i++) {
302 int num=a[i][p[i]];
303 cout<<"("<<hitCol[num].getLayerId()<<","<<hitCol[num].getWireId()<<")"<<endl;
304 recHitCol.push_back( hitCol[num] );
305 }
306 fit_global2D(recHitCol);
307 numall++;
308 }
309}
310*/
311int HoughTrack::fit3D(){
312 sortHit();
313 // cout<<"in fit3d"<<endl;
314 int nhit_zs = calzs();
315 if( nhit_zs >=3 ) fitzs();
316 else {
317 cout<<" not enough ster hit in zs fit"<<endl;
318 _tanl=0.;_z0=0.;
319 }
320 outerHit();
321 if(m_debug>0) { cout<< " before 3d fit "<<endl; this->printRecHit();}
322 if(m_globalfit)_stat3D = fit_global3D(0);
323 else _stat3D =1;
324 return _stat3D;
325}
326
327int HoughTrack::fit3D_inner(){ //for multi turn track
328 sortHit();
329 int nhit_zs = calzs();
330 cutMultiCirHit();
331 if( nhit_zs >=3 ) fitzs();
332 else {
333 // cout<<" not enough ster hit in zs fit"<<endl;
334 _tanl=0.;_z0=0.;
335 }
336 cutMultiCirHit_after_zs();
337 outerHit();
338 if(m_debug>0) { cout<< " before 3d fit "<<endl; this->printRecHit();}
339 if(m_globalfit)_stat3D = fit_global3D(0);
340 else _stat3D =1;
341 return _stat3D;
342}
343
344int HoughTrack::cutMultiCirHit(){
345 for(int i=0;i<_recHitVec.size();i++){
346 if( _recHitVec[i].getLayerId()>8 ) continue;
347 if( (fabs(_recHitVec[i].getzAmb(0))>10) && (fabs(_recHitVec[i].getzAmb(1))>10) ) {
348 _recHitVec[i].setflag(1);
349 // cout<<"reject ("<<_recHitVec[i].getLayerId()<<","<<_recHitVec[i].getWireId()<<")"<<_recHitVec[i].getCirList()<<endl;
350 }
351 }
352 _Hough3D.setRecHit(_recHitVec);
353}
354
355int HoughTrack::cutMultiCirHit_after_zs(){
356 for(int i=0;i<_recHitVec.size();i++){
357 int layer = _recHitVec[i].getLayerId();
358 double zl = _recHitVec[i].getzAmb(0);
359 double zr = _recHitVec[i].getzAmb(1);
360 double sl = _recHitVec[i].getsAmb(0);
361 double sr = _recHitVec[i].getsAmb(1);
362 double dl = zl - (sl*_tanl+_z0);
363 double dr = zr - (sr*_tanl+_z0);
364 if (layer<=8 && fabs(dl)>10 && fabs(dr)>10 ) _recHitVec[i].setflag(1);
365 if (layer>8 && fabs(dl)>20 && fabs(dr)>20 ) _recHitVec[i].setflag(1);
366 // if (layer>30) _recHitVec[i].setflag(1);
367 }
368 _Hough3D.setRecHit(_recHitVec);
369}
370
371int HoughTrack::fitLeast(){
372 _Hough2D=Hough2D(_recHitVec,_bunchTime);
373 double circleR = fabs(1./(_maprho));
374 double circleX = (1./_maprho)*cos(_maptheta);
375 double circleY = (1./_maprho)*sin(_maptheta);
376 double d0 = sqrt( circleX*circleX + circleY*circleY )-circleR;
377 double phi0 = atan2(circleY,circleX)+M_PI/2.;
378 double omega = 1/circleR;
379 if(_charge==-1){
380 //d0=d0;
381 d0=-d0;
382 omega=-1.*fabs(omega);
383 }
384 if(_charge==1){
385 //d0=-d0;
386 d0=d0;
387 omega=1.*fabs(omega);
388 phi0=phi0-M_PI;
389 }
390 _Hough2D.setCharge(_charge);
391 _Hough2D.setCirX( circleX);
392 _Hough2D.setCirY( circleY );
393 _Hough2D.setCirR( circleR );
394 _Hough2D.setD0( d0 );
395 _Hough2D.setPhi0( phi0 );
396 _Hough2D.setOmega( omega );
397 _Hough2D.setPt( circleR/333.567 );
398
399 //fill HoughTrack
400 _centerX=_Hough2D.getCirX();
401 _centerY=_Hough2D.getCirY();
402 _centerR=_Hough2D.getCirR();
403 //_pt=_Hough2D.getPt();
404 _pt2D=_Hough2D.getPt();
405 _d0=_Hough2D.getD0();
406 _omega=_Hough2D.getOmega();
407 _phi0=_Hough2D.getPhi0();
408 hitOnTrack();
409 if(m_debug>0) {cout<<" after least 2d "<<endl; printRecHit();}
410 return 1;
411
412}
413
414int HoughTrack::fit_global2D(vector<HoughRecHit>& recHitVec){
415 //_Hough2D=Hough2D(_Hough2D,_recHitVec,bunchtime);
416 //_Hough2D.setRecHit(_recHitVec);
417 _Hough2D.setRecHit(recHitVec);
418 //map<int, const HepVector >::iterator map_it = g_tkParTruth.begin();
419 //cout<<"size:"<<g_tkParTruth.size()<<endl;
420 //cout<<map_it->first<<" "<<map_it->second<<endl;
421 //double d0 = -(map_it->second)[0];
422 //double phi0 = (map_it->second)[1]+M_PI/2;
423 //double omega = (map_it->second)[2]/333.567;
424 //double z0 = (map_it->second)[3];
425 //double tanl = (map_it->second)[4];
426 //while(phi0> M_PI)phi0-=2*M_PI;
427 //while(phi0<-M_PI)phi0+=2*M_PI;
428 //cout<<" fit3D:"<<setw(12)<<d0<<" "<<setw(12)<<phi0<<" "<<setw(12)<<omega<<" "<<setw(12)<<z0<<" "<<setw(12)<<tanl<<endl;
429 //_Hough2D.setD0(d0);
430 //_Hough2D.setPhi0(phi0);
431 //_Hough2D.setOmega(omega);
432 //_Hough3D.setZ0(z0);
433 //_Hough3D.setTanl(tanl);
434 int Stat_2d=_Hough2D.fit();
435 p_trk2D = _Hough2D.getTrk();
436 if(Stat_2d==1){
437 _centerX=_Hough2D.getCirX();
438 _centerY=_Hough2D.getCirY();
439 _centerR=_Hough2D.getCirR();
440 //_pt=_Hough2D.getPt();
441 _pt2D=_Hough2D.getPt();
442 _d0=_Hough2D.getD0();
443 _omega=_Hough2D.getOmega();
444 _phi0=_Hough2D.getPhi0();
445 _chi2_aver2D = _Hough2D.getChi2_2D();
446 _nfit2D= _Hough2D.getNfit();
447 if(m_debug>0) cout<<"pt after global 2d: "<<_pt2D<<endl;
448 if(m_debug>0) cout<<"after global 2d "<<endl;
449 hitOnTrack();
450 if(m_debug>0) { cout<<" after global 2d "<<endl; printRecHit();}
451 }
452 else {
453 if(m_debug>0) cout<<" 2d global fail"<<endl;
454 //cout<<" 2d global fail"<<endl;
455 }
456 //cout<<_d0<<" "<<_phi0<<" "<<_omega<<" "<<endl;
457 return Stat_2d;
458}
459
460int HoughTrack::fit_global3D(int time){
461 if( time==0 ) _Hough3D=Hough3D(_Hough2D,_recHitVec,_bunchTime,_tanl,_z0,vec_mdcHit);
462 //map<int, const HepVector >::iterator map_it = g_tkParTruth.begin();
463 //cout<<"size:"<<g_tkParTruth.size()<<endl;
464 //cout<<map_it->first<<" "<<map_it->second<<endl;
465 //double d0 = -(map_it->second)[0];
466 //double phi0 = (map_it->second)[1]+M_PI/2;
467 //double omega = (map_it->second)[2]/333.567;
468 //double z0 = (map_it->second)[3];
469 //double tanl = (map_it->second)[4];
470 //while(phi0> M_PI)phi0-=2*M_PI;
471 //while(phi0<-M_PI)phi0+=2*M_PI;
472 //cout<<" fit3D:"<<setw(12)<<d0<<" "<<setw(12)<<phi0<<" "<<setw(12)<<omega<<" "<<setw(12)<<z0<<" "<<setw(12)<<tanl<<endl;
473 //_Hough3D.setD0(d0);
474 //_Hough3D.setPhi0(phi0);
475 //_Hough3D.setOmega(omega);
476 //_Hough3D.setZ0(z0);
477 //_Hough3D.setTanl(tanl);
478 if(m_debug>0) _Hough3D.print();
479 // _Hough3D.setCharge(_charge);
480 int Stat_3d=_Hough3D.fit();
481 p_trk = _Hough3D.getTrk();
482 if(Stat_3d==1){
483 _centerX=_Hough3D.getCirX();
484 _centerY=_Hough3D.getCirY();
485 _centerR=_Hough3D.getCirR();
486 _pt3D=_Hough3D.getPt();
487 _p=_Hough3D.getP();
488 _pz=_Hough3D.getPz();
489 _d0=_Hough3D.getD0();
490 _omega=_Hough3D.getOmega();
491 _phi0=_Hough3D.getPhi0();
492 _tanl=_Hough3D.getTanl();
493 _z0=_Hough3D.getZ0();
494 _chi2_aver = _Hough3D.getChi2();
495 _nfit= _Hough3D.getNfit();
496 if(m_debug>0) cout<<"pt after global 3d: "<<_pt3D<<endl;
497 // hitOnTrack();
498 }
499 else {
500 if(m_debug>0) cout<<" 3d global fail"<<endl;
501 //cout<<" 3d global fail"<<endl;
502 }
503 return Stat_3d;
504}
505
506void HoughTrack::hitOnTrack(){
507 if(m_debug>0) cout<<" calculate hit on track "<<endl;
508 for(int ihit=0;ihit<_recHitVec.size();ihit++){
509 std::pair<double,double> theta_l = calcuArcTrack( (_recHitVec[ihit]) );
510 int flag = judge_half(_recHitVec[ihit]);
511 double dist=calcuDistToTrack( (_recHitVec[ihit]));
512 double distToCir=calcuDistToCir( (_recHitVec[ihit]));
513 _recHitVec[ihit].setflag(flag); //FIXME
514 _recHitVec[ihit].setDisToTrack(dist);
515 _recHitVec[ihit].setDisToCir(distToCir);
516 _recHitVec[ihit].setRet(theta_l);
517 }
518 // cut_axial_inner();
519}
520
521int HoughTrack::judge_half(const HoughRecHit& hit){
522 int cir;
523 double xhit =hit.getMidX();
524 double yhit =hit.getMidY();
525 double x_cir=_Hough2D.getCirX();
526 double y_cir=_Hough2D.getCirY();
527 double r_cir=_Hough2D.getCirR();
528 if(_charge==-1){
529 if( (x_cir*yhit - y_cir*xhit >=0) ) cir=0;
530 else cir=1;
531 }
532 else if(_charge==1){
533 if( (x_cir*yhit - y_cir*xhit <=0) ) cir=0;
534 else cir=1;
535 }
536 else cout<<" charge is not set !!!!!!!!!!"<<endl;
537
538 return cir;
539}
540
541double HoughTrack::calcuDistToTrack(const HoughRecHit& hit){
542 double xhit =hit.getMidX();
543 double yhit =hit.getMidY();
544 double x_cir=_Hough2D.getCirX();
545 double y_cir=_Hough2D.getCirY();
546 double r_cir=_Hough2D.getCirR();
547 double dist=sqrt( pow( (xhit-x_cir),2)+pow( (yhit-y_cir),2) )-r_cir;
548 return dist;
549}
550
551double HoughTrack::calcuDistToCir(const HoughRecHit& hit){
552 double xhit =hit.getMidX();
553 double yhit =hit.getMidY();
554 double x_cir=_Hough2D.getCirX();
555 double y_cir=_Hough2D.getCirY();
556 double r_cir=_Hough2D.getCirR();
557 double dist=sqrt( pow( (xhit-x_cir),2)+pow( (yhit-y_cir),2) );
558 return dist;
559}
560//double HoughTrack::InterSect(const HoughRecHit& hit){
561// double xeast = _rechit->getEastPoint().x();
562// double xwest = _rechit->getWestPoint().x();
563// double yeast = _rechit->getEastPoint().y();
564// double ywest = _rechit->getWestPoint().y();
565// double k = (ywest-yeast)/(xwest-xeast);
566// double b = -k*xeast+yeast;
567// double drift = _rechit->getDriftDist();
568// double x_cir=_Hough2D.getCirX();
569// double y_cir=_Hough2D.getCirY();
570// double r_cir=_Hough2D.getCirR();
571//}
572
573std::pair<double,double> HoughTrack::calcuArcTrack(const HoughRecHit& hit){
574 double xhit =hit.getMidX();
575 double yhit =hit.getMidY();
576 double x_cir=_centerX;
577 double y_cir=_centerY;
578 double r_cir=_centerR;
579 std::pair<double,double> theta_l;
580
581 double theta_temp;
582 double l_temp;
583 //if(x_cir==0||xhit-x_cir==0){
584 //theta_temp=0;
585 //}
586 //else{
587 theta_temp=M_PI-atan2(yhit-y_cir,xhit-x_cir)+atan2(y_cir,x_cir);
588 if(theta_temp>2*M_PI){
589 theta_temp=theta_temp-2*M_PI;
590 }
591 if(theta_temp<0){
592 theta_temp=theta_temp+2*M_PI;
593 }
594 //}
595 if(_charge==-1) {
596 l_temp=r_cir*theta_temp;
597 }
598 if(_charge==1) {
599 theta_temp=2*M_PI-theta_temp;
600 l_temp=r_cir*(theta_temp);
601 }
602 theta_l.first=theta_temp;
603 theta_l.second=l_temp;
604 return theta_l;
605}
606
607int HoughTrack::calzs(){
608 int n_zs=0;
609 for( int i =0; i<_recHitVec.size(); i++){
610 if ( _recHitVec[i].getSlayerType() ==0 || _recHitVec[i].getflag()!=0 ) continue;
611 if ( _recHitVec[i].getDetectorType()==MDC){
612 // if ( _recHitVec[i].getLayerId()>=8) continue;
613 HoughStereo zs( _bunchTime, &_Hough2D, &(_recHitVec[i]) );
614 // cout<<"("<<_recHitVec[i].getLayerId()<<","<<_recHitVec[i].getWireId()<<") stat "<<stat<<" style "<<_recHitVec[i].getStyle() <<endl;
615 // cout<<endl;
616
617 //zs.setAmb(-1); //right
618 //int stat_left = zs.cald();
619
620 //zs.setAmb(1); //left
621 //int stat_right= zs.cald();
622 //if( stat_left==-1 && stat_right==-1 ) _recHitVec[i].setflag(-999);
623
624 int stat= zs.cald();
625 _recHitVec[i].setnsol(stat);
626 if( stat==0 ) _recHitVec[i].setflag(-999);
627
628 zs.setRecHit();
629 _recHitVec[i].setAmb(-999); // check
630 //cout<<__FILE__<<" ("<<_recHitVec[i].getLayerId()<<", "<<_recHitVec[i].getWireId()<<") "<<_recHitVec[i].getflag()<<" "<<stat<<endl;
631 //if( _recHitVec[i].getLayerId()<8 ) n_zs++;
632 n_zs++;
633 if(m_debug>0) zs.print();
634 //zs.print();
635 }else if( _recHitVec[i].getDetectorType()==CGEM){
636 double z=(_recHitVec[i].getCluster())->getRecZ()/10.;
637 std::pair<double,double> theta_l = calcuArcTrack(_recHitVec[i]);
638 double l=theta_l.second;
639 _recHitVec[i].setzAmb( 0 , z);
640 _recHitVec[i].setsAmb( 0 , l);
641 _recHitVec[i].setzAmb( 1 , z);
642 _recHitVec[i].setsAmb( 1 , l);
643 _recHitVec[i]. setzPos(z);
644 _recHitVec[i]. setsPos(l);
645 HepPoint3D leftPoint(_recHitVec[i].getMidX(),_recHitVec[i].getMidY(),z);
646 HepPoint3D rightPoint(_recHitVec[i].getMidX(),_recHitVec[i].getMidY(),z);
647 _recHitVec[i].setLeftPoint(leftPoint);
648 _recHitVec[i].setRightPoint(rightPoint);
649 n_zs++;
650 }
651 }
652 // int naddStero = cutNoise_inner();
653 //return naddStero;
654 return n_zs;
655}
656void HoughTrack::fitzs(){
657 //for( int i =0; i<_recHitVec.size(); i++)cout<<__FILE__<<" ("<<_recHitVec[i].getLayerId()<<", "<<_recHitVec[i].getWireId()<<") "<<_recHitVec[i].getflag()<<endl;
658 HoughZsFit zsfit(&_recHitVec);
659 //for( int i =0; i<_recHitVec.size(); i++)cout<<__FILE__<<" ("<<_recHitVec[i].getLayerId()<<", "<<_recHitVec[i].getWireId()<<") "<<_recHitVec[i].getflag()<<endl;
660 _tanl=zsfit.getTanl();
661 _z0=zsfit.getZ0();
662 _tanl_zs=zsfit.getTanl();
663 _z0_zs=zsfit.getZ0();
664 t_pro_correct = zsfit.getPro();
665 if(m_debug>0) printRecHit();
666}
667
668int HoughTrack::getHitNum(int select ) const{
669 int size=_recHitVec.size();
670 int n0,n1,n2,n3,n4,n5,n6;
671 n0=n1=n2=n3=n4=n5=n6=0;
672 for(int i=0;i<size;i++){
673 int cir= _recHitVec[i].getCirList();
674 int style= _recHitVec[i].getStyle();
675 n0++;
676 if( style == 1 ) n1++;
677 if( style == 2 ) n2++;
678 if(_recHitVec[i].getDetectorType()==MDC){
679 int index = _recHitVec[i].digi()->getTrackIndex();
680 if( index < 0 ) n3++;
681 if( index >=0 && cir<=1 && style==0) n4++;
682 if( index <0 || cir>1 ) n5++;
683 if( index >=0 && cir==0 && style==0) n6++;
684 }
685 }
686 if( select == 0 ) return n0;
687 if( select == 1 ) return n1;
688 if( select == 2 ) return n2;
689 if( select == 3 ) return n3;
690 if( select == 4 ) return n4;
691 if( select == 5 ) return n5;
692 if( select == 6 ) return n6;
693}
694int HoughTrack::getHitNumA(int select ) const{
695 int size=_recHitVec.size();
696 int n0,n1,n2,n3,n4,n5,n6;
697 n0=n1=n2=n3=n4=n5=n6=0;
698 for(int i=0;i<size;i++){
699 int cir= _recHitVec[i].getCirList();
700 int type = _recHitVec[i].getSlayerType();
701 int style= _recHitVec[i].getStyle();
702 if( type ==0 ) n0++;
703 if( type==0 && style == 1 ) n1++;
704 if( type==0 && style == 2 ) n2++;
705 if(_recHitVec[i].getDetectorType()==MDC){
706 int index = _recHitVec[i].digi()->getTrackIndex();
707 if( type==0 && index < 0 ) n3++;
708 if( type==0 && index >=0 && cir<=1 && style==0 ) n4++;
709 if( type==0 && (index <0 || cir>1) ) n5++;
710 if( type==0 && index >=0 && cir==0 && style==0 ) n6++;
711 }
712 }
713 if( select == 0 ) return n0;
714 if( select == 1 ) return n1;
715 if( select == 2 ) return n2;
716 if( select == 3 ) return n3;
717 if( select == 4 ) return n4;
718 if( select == 5 ) return n5;
719 if( select == 6 ) return n6;
720}
721int HoughTrack::getHitNumS(int select ) const{
722 int size=_recHitVec.size();
723 int n0,n1,n2,n3,n4,n5,n6;
724 n0=n1=n2=n3=n4=n5=n6=0;
725 for(int i=0;i<size;i++){
726 int cir= _recHitVec[i].getCirList();
727 int type = _recHitVec[i].getSlayerType();
728 int style= _recHitVec[i].getStyle();
729 if( type !=0 ) n0++;
730 if( type!=0 && style == 1 ) n1++;
731 if( type!=0 && style == 2 ) n2++;
732 if(_recHitVec[i].getDetectorType()==MDC){
733 int index = _recHitVec[i].digi()->getTrackIndex();
734 if( type!=0 && index < 0 ) n3++;
735 if( type!=0 && index >=0 && cir<=1 && style==0) n4++;
736 if( type!=0 && (index <0 || cir>1) ) n5++;
737 if( type!=0 && index >=0 && cir==0 && style==0) n6++;
738 }
739 }
740 if( select == 0 ) return n0;
741 if( select == 1 ) return n1;
742 if( select == 2 ) return n2;
743 if( select == 3 ) return n3;
744 if( select == 4 ) return n4;
745 if( select == 5 ) return n5;
746 if( select == 6 ) return n6;
747}
748
749int HoughTrack::judgeTrack(){
750 // hit number test
751 int axialHit=0;
752 // int stereoHit=0;
753 // int innerHit=0;
754 vector<int> vec_layer;
755 for( int i =0;i<_recHitVec.size();i++){
756 int layer = _recHitVec[i].getLayerId();
757 int slant= _recHitVec[i].getSlayerType();
758 if( 0==slant ) axialHit++;
759 // if( 0!=slant ) stereoHit++;
760 // if( layer<8 ) innerHit++;
761 }
762 // if ( axialHit <5 || stereoHit <3 || innerHit ==0 ) return 0;
763 if ( axialHit <5 ) return 0;
764 else return 1;
765 // // continues hit test
766 // bool continues=false;
767 // for(int i =0;i<vec_layer.size();i++){
768 // int layer = vec_layer[i];
769 // vector<int>::iterator iter = find( vec_layer.begin(),vec_layer.end(),layer+1 );
770 // if( iter != vec_layer.end() ){
771 // vector<int>::iterator iter2 = find( vec_layer.begin(),vec_layer.end(),layer-1 );
772 // if( iter2 != vec_layer.end() ) continues = true;
773 // }
774 // }
775}
776
777int HoughTrack::find_stereo_hit(){
778 vector<HoughRecHit> vec_rec; // store hit temp
779 int naddStero=0;
780 for(int ihit=0;ihit<_houghList.size();ihit++){
781 const HoughHit hit= _houghList[ihit];
782 if ( hit.getDetectorType() == CGEM) continue;
783 if ( hit.getSlayerType()==0 ) continue;
784 if ( hit.driftTime(_bunchTime,0)>800 ) continue;
785 HoughRecHit p_hit( hit,0,0,1);
786 p_hit.setPtr2D(&_Hough2D);
787 std::pair<double,double> theta_l = calcuArcTrack( p_hit );
788 double dist=calcuDistToTrack( p_hit );
789 double distToCir=calcuDistToCir( p_hit );
790 int flag = judge_half( p_hit );
791 int layer = p_hit.getLayerId();
792 int wire= p_hit.getWireId();
793 double disCut;
794 if(layer<8) disCut=4;
795 else disCut=6;
796 // double dist=sqrt( pow( (xhit-x_cir),2)+pow( (yhit-y_cir),2) )-r_cir;
797 // cout<<"find stereo hit dist?("<<layer<<" ,"<<wire<<" ) "<<dist<<endl;
798 if( m_debug >0 )cout<<"("<<layer<<","<<wire<<") "<<" rec hit dist theta "<<dist<<" "<< theta_l.first<<endl;
799 //cout<<__FILE__<<" ("<<layer<<","<<wire<<") "<<flag<<" rec hit dist theta "<<dist<<" "<< theta_l.first<<endl;
800 if( fabs(dist) < disCut ) {
801 p_hit.setDisToTrack(dist);
802 p_hit.setDisToCir(distToCir);
803 p_hit.setRet(theta_l);
804 p_hit.setflag(flag);
805 // if( p_hit.getflag()==0 ) vec_rec.push_back(p_hit);
806 // vec_rec.push_back(p_hit);
807 if( p_hit.getflag()!=0 ) continue;
808 _recHitVec.push_back(p_hit);
809 naddStero++;
810 }
811 }
812 return naddStero;
813}
814/*
815int HoughTrack::cut_axial_inner(){
816 for( int ilay=0;ilay<43;ilay++){
817 //if( (ilay<7&&ilay<20) || ilay>35 ) continue;
818 if ((ilay>=0&&ilay<=7)||(ilay>=20&&ilay<=35)) continue;
819 // for( int ilay=0;ilay<8;ilay++)
820 double disToCir=9999;
821 int max=-999;
822 int count=0;
823 int count_cut=0;
824 for( int ihit=0;ihit<_recHitVec.size();ihit++){
825 if( _recHitVec[ihit].getLayerId()==ilay && _recHitVec[ihit].getflag()==0 ) {
826 count++;
827 if( (fabs(_recHitVec[ihit].getDisToCir())<disToCir) ) {
828 disToCir= fabs(_recHitVec[ihit].getDisToCir()); max=ihit;
829 count_cut=1;
830 }
831 }
832 }
833 if( m_debug >0 ) cout<<"ilay count count_cut "<<ilay<<" "<<count<<" "<<count_cut<<endl;
834 if( count_cut !=0 && count>1 ) _recHitVec[max].setflag(-999);
835 if( m_debug >0 && count_cut!=0 && count>1) cout<<"delete ("<<_recHitVec[max].getLayerId()<<","<<_recHitVec[max].getWireId()<<")"<<endl;
836 }
837 int size_of_axial;
838 for( int ihit=0;ihit<_recHitVec.size();ihit++){
839 if(_recHitVec[ihit].getSlayerType()==0 && _recHitVec[ihit].getflag()==0 ) size_of_axial++;
840 }
841 return size_of_axial;
842}
843
844int HoughTrack::cutNoise_inner(){
845 // for( int ilay=0;ilay<43;ilay++)
846 // if( (ilay>7&&ilay<20) || ilay>35 ) continue;
847 for( int ilay=0;ilay<8;ilay++){
848 double z_0=0;
849 double z_1=0;
850 int min=-999;
851 int count=0;
852 int count_cut=0;
853 // cout<<" size "<<_recHitVec.size()<<endl;
854 for( int ihit=0;ihit<_recHitVec.size();ihit++){
855 if( _recHitVec[ihit].getLayerId()==ilay && _recHitVec[ihit].getflag()==0) {
856 count++;
857 if( ((fabs(_recHitVec[ihit].getzAmb(0))>z_0) && (fabs(_recHitVec[ihit].getzAmb(1)) > z_1)) ) {
858 if( z_0!=0 && z_1!=0 ) count_cut++;
859 z_0= fabs(_recHitVec[ihit].getzAmb(0));z_1= fabs(_recHitVec[ihit].getzAmb(1)) ; min=ihit;
860 }
861 else if( ((fabs(_recHitVec[ihit].getzAmb(0))<z_0) && (fabs(_recHitVec[ihit].getzAmb(1)) < z_1)) ) {
862 count_cut++;
863 }
864 }
865 }
866 if( m_debug >0 ) cout<<"ilay count count_cut "<<ilay<<" "<<count<<" "<<count_cut<<endl;
867 if( count_cut>0 && count>1 ) _recHitVec[min].setflag(-999);
868 if( m_debug >0 && count_cut>0 && count>1) cout<<"delete ("<<_recHitVec[min].getLayerId()<<","<<_recHitVec[min].getWireId()<<")"<<endl;
869 }
870 int size_of_stereo;
871 for( int ihit=0;ihit<_recHitVec.size();ihit++){
872 if(_recHitVec[ihit].getSlayerType()!=0 && _recHitVec[ihit].getflag()==0 ) size_of_stereo++;
873 }
874 return size_of_stereo;
875}
876*/
877int HoughTrack::trackCharge2D(){
878 int n_neg(0);
879 int n_pos(0);
880 for(int ihit=0;ihit<_recHitVec.size();ihit++){
881 if(_recHitVec[ihit].getSlayerType()!=0 ) continue;
882 double xhit =_recHitVec[ihit].getMidX();
883 double yhit =_recHitVec[ihit].getMidY();
884 double x_cir=_Hough2D.getCirX();
885 double y_cir=_Hough2D.getCirY();
886 double r_cir=_Hough2D.getCirR();
887 if( (x_cir*yhit - y_cir*xhit >=0) ) n_neg++;
888 if( (x_cir*yhit - y_cir*xhit <=0) ) n_pos++;
889 }
890 if( m_debug >0 ) cout<<" in track charge 2d "<<n_neg<<" "<<n_pos<<endl;
891 if( (_charge==-1&&n_neg<3) || (_charge==1&&n_pos<3) ) return 0;
892 else if( _charge ==0 ) cout<<" wrong ! in track charge 2D not set charge "<<endl;
893 else return 1;
894}
895
896int HoughTrack::trackCharge3D(){
897 int n_neg(0);
898 int n_pos(0);
899 for(int ihit=0;ihit<_recHitVec.size();ihit++){
900 if(_recHitVec[ihit].getSlayerType()==0 ) continue;
901 if(_recHitVec[ihit].getflag()!=0)continue;
902 //if(_recHitVec[ihit].getLayerId()>7 ) continue;
903 double xhit =_recHitVec[ihit].getMidX();
904 double yhit =_recHitVec[ihit].getMidY();
905 double x_cir=_Hough2D.getCirX();
906 double y_cir=_Hough2D.getCirY();
907 double r_cir=_Hough2D.getCirR();
908 if( (x_cir*yhit - y_cir*xhit >=0) ) n_neg++;
909 if( (x_cir*yhit - y_cir*xhit <=0) ) n_pos++;
910 }
911 if( m_debug >0 ) cout<<" in track charge 3d "<<n_neg<<" "<<n_pos<<endl;
912 if( (_charge==-1&&n_neg<2) || (_charge==1&&n_pos<2) ) return 0;
913 else if( _charge ==0 ) cout<<" wrong ! in track charge 3D not set charge "<<endl;
914 else return 1;
915}
916/*
917 void HoughTrack::Leastfit(vector<double> u,vector<double> v,double& k ,double& b){
918 int N = u.size();
919 if (N<=2) return;
920 double *x=new double[N];
921 double *y=new double[N];
922 for(int i=0;i<N;i++){
923 x[i]=u[i];
924 y[i]=v[i];
925 }
926
927 TF1 *fpol1=new TF1("fpol1","pol1",-50,50);
928 TGraph *tg=new TGraph(N,x,y);
929 tg->Fit("fpol1","QN");
930 double ktemp =fpol1->GetParameter(1);
931 double btemp =fpol1->GetParameter(0);
932 k = ktemp;
933 b = btemp;
934 delete []x;
935 delete []y;
936 delete fpol1;
937 delete tg;
938 }
939 */
940void HoughTrack::cald_layer(){
941 /*
942 //in truth
943 double k,b,theta,rho,x_cross,y_cross;
944 vector<double> vtemp,utemp;
945 std::vector<HoughRecHit>::iterator iter = _recHitVec.begin();
946 for(int iHit=0;iter!= _recHitVec.end(); iter++,iHit++){
947 const HoughRecHit h = (*iter);
948 // if( h.getCirList()!=0 ) continue;
949 if ( h.digi()->getTrackIndex()>=0 && h.getStyle()==0&& h.getSlayerType()==0 && h.getCirList()==0 && utemp.size()<10) // ??use 2nd half
950 {
951 utemp.push_back(h.getUTruth());
952 vtemp.push_back(h.getVTruth());
953 }
954 }
955 Leastfit(utemp,vtemp,k,b);
956 //calcu truth
957 //k,b from truth
958 x_cross = -b/(k+1/k);
959 y_cross = b/(1+k*k);
960 rho=sqrt(x_cross*x_cross+y_cross*y_cross);
961 theta=atan2(y_cross,x_cross);
962 //
963 */
964 //cout<<"track center position "<<_centerX<<" "<<_centerY<<endl;
965 std::vector<HoughRecHit>::iterator iter0= _recHitVec.begin();
966 for(; iter0!= _recHitVec.end(); iter0++){
967 HoughRecHit *hit = &(*iter0);
968 if( hit->getSlayerType()!=0 ) continue;
969 //if(hit->getDetectorType() == CGEM)continue;
970 // if( hit->getCirList()!=0 ) continue; // use in learn distribute
971 //double cirr_hit = hit->getDriftDistTruth();
972 int layer = hit->getLayerId();
973 int wire = hit->getWireId();
974 double cirx_hit = hit->getMidX();
975 double ciry_hit = hit->getMidY();
976 double cirr_hit = hit->getDriftDist();
977 //if( _charge == 1 && (cirx_track*ciry_hit - ciry_track*cirx_hit>=0) ) continue; // suppose charge -1 FIXME
978 //if( _charge ==-1 && (cirx_track*ciry_hit - ciry_track*cirx_hit<=0) ) continue; // suppose charge -1 FIXME
979 double deltaD ;
980 if(hit->getDetectorType() == MDC){
981 double l1l2 = sqrt( (cirx_hit-_centerX)*(cirx_hit-_centerX)+(ciry_hit-_centerY)*(ciry_hit-_centerY) );
982 if( l1l2>_centerR ) deltaD = l1l2-_centerR-cirr_hit;
983 if( l1l2<=_centerR ) deltaD = l1l2-_centerR+cirr_hit;
984 }else if(hit->getDetectorType() == CGEM){
985 double r_cgem = hit->getCgemGeomSvc()->getCgemLayer(layer)->getMiddleROfGapD()/10.;
986 double phi_cross = intersect_cylinder(_charge,_centerX,_centerY,r_cgem);
987 double phih_cluster = hit->getCluster()->getrecphi();
988 double dphi = phih_cluster - phi_cross;
989 if(dphi>M_PI)dphi-=2*M_PI;
990 if(dphi<-1*M_PI)dphi+=2*M_PI;
991 deltaD = r_cgem*dphi;
992 }
993 hit->setDeltaD(deltaD);
994 //cout<<layer<<" "<<deltaD<<endl;
995 //cout<<l1l2<<" "<<hit->getDisToCir()<<endl;
996
997
998 //cal flight length
999 /*
1000 double theta_temp;
1001 double l_temp;
1002 if(_centerX==0||cirx_hit-_centerX==0){
1003 theta_temp=0;
1004 }
1005 else{
1006 theta_temp=M_PI-atan2(ciry_hit-_centerY,cirx_hit-_centerX)+atan2(_centerY,_centerX);
1007 if(theta_temp>2*M_PI){
1008 theta_temp=theta_temp-2*M_PI;
1009 }
1010 if(theta_temp<0){
1011 theta_temp=theta_temp+2*M_PI;
1012 }
1013 }
1014 //cout<<" charge "<<_charge <<" "<<theta_temp<<endl;
1015 if(_charge==-1) {
1016 l_temp=_centerR*theta_temp;
1017 }
1018 if(_charge==1) {
1019 theta_temp=2*M_PI-theta_temp;
1020 l_temp=_centerR*(theta_temp);
1021 }
1022 // cout<<"("<<hit->getLayerId()<<","<<hit->getWireId()<<") "<<l_temp<<endl;
1023 hit->setFltLen(l_temp);
1024 */
1025 // cout<<"int map deltaD: ("<<hit->getLayerId()<<","<<hit->getWireId()<<")"<< hit->getDeltaD()<<endl;
1026 }
1027}
1028
1029//calculate phi value, only valid for the first half circle
1030double HoughTrack::intersect_cylinder(double r_cylinder, double r_center, double phi_center, int charge)
1031{
1032 double phi;
1033 if(charge==0)return 9999;
1034 while(phi_center<0) phi_center += 2*M_PI;
1035 while(phi_center>2*M_PI) phi_center -= 2*M_PI;
1036 if(r_center<r_cylinder/2) return -9999;
1037 double dphi = acos( r_cylinder/(2*r_center) );
1038 if(charge<0) phi = phi_center + dphi;
1039 else phi = phi_center - dphi;
1040 while(phi<0) phi += 2*M_PI;
1041 while(phi>2*M_PI) phi -= 2*M_PI;
1042 return phi;
1043}
1044double HoughTrack::intersect_cylinder(int charge, double x_center, double y_center, double r_cylinder)
1045{
1046 double phi;
1047 double phi_center = atan2(y_center,x_center);
1048 double r_center = sqrt(x_center*x_center+y_center*y_center);
1049 if(charge==0)return 9999;
1050 while(phi_center<0) phi_center += 2*M_PI;
1051 while(phi_center>2*M_PI) phi_center -= 2*M_PI;
1052 if(r_center<r_cylinder/2) return -9999;
1053 double dphi = acos( r_cylinder/(2*r_center) );
1054 if(charge<0) phi = phi_center + dphi;
1055 else phi = phi_center - dphi;
1056 while(phi<0) phi += 2*M_PI;
1057 while(phi>2*M_PI) phi -= 2*M_PI;
1058 return phi;
1059}
1060
1061int HoughTrack::find_pair_hit(){
1062 int nster1=0;
1063 int nster2=0;
1064 int nster3=0;
1065 int naxial1=0;
1066 int naxial2=0;
1067 int naxial3=0;
1068 for(int ihit=0;ihit<_houghList.size();ihit++){
1069 const HoughHit hit= _houghList[ihit];
1070 if ( hit.getDetectorType() == CGEM) continue;
1071 if ( hit.driftTime(_bunchTime,0)>1000 ) continue;
1072 if ( hit.getLayerId()>=24 ) continue;
1073 HoughRecHit p_hit( hit,0,0,1);
1074 p_hit.setPtr2D(&_Hough2D);
1075 std::pair<double,double> theta_l = calcuArcTrack( p_hit );
1076 double dist=calcuDistToTrack( p_hit );
1077 double distToCir=calcuDistToCir( p_hit );
1078 int flag = judge_half( p_hit );
1079 int layer = p_hit.getLayerId();
1080 int wire= p_hit.getWireId();
1081 int slayer= p_hit.getSlayerType();
1082 double disCut;
1083 if( slayer == 0 ) disCut = 6;
1084 else {
1085 if(layer<8) disCut=6;
1086 else disCut=6;
1087 }
1088 if( m_debug >0 )cout<<"("<<layer<<","<<wire<<") "<<" pair dist flag "<<dist<<" "<< flag<<endl;
1089 if( fabs(dist) < disCut ) {
1090 if( flag !=1 ) continue; //find hit in the 2nd half
1091 if ( layer<4) nster1++;
1092 else if ( layer<8) nster2++;
1093 else if ( layer<12) naxial1++;
1094 else if ( layer<16) naxial2++;
1095 else if ( layer<20) naxial3++;
1096 else nster3++;
1097 }
1098 }
1099 if(m_debug >0){
1100 cout <<"naxial_1 "<<naxial1<<endl;
1101 cout <<"naxial_2 "<<naxial2<<endl;
1102 cout <<"naxial_3 "<<naxial3<<endl;
1103 cout <<"stereo_1 "<<nster1<<endl;
1104 cout <<"stereo_2 "<<nster2<<endl;
1105 cout <<"stereo_3 "<<nster3<<endl;
1106 }
1107 //if (nster1>=2 && nster2>=2 && nster3>=2 && naxial1>=2 && naxial2>=2 && naxial3>=2 ) return 1;
1108 if ( nster1>=2 && nster3>=2 && naxial1>=2 && naxial2>=2 && naxial3>=2 ) return 1;
1109 else return 0;
1110}
1111void HoughTrack::print(){
1112 cout<<"print HoughTrack : "<<p_trk->id()<<endl;
1113 cout<<"peakNum:"<<_centerPeak.getPeakNum()<<" peak height:"<<_centerPeak.peakHeight()<<" rho,theta:"<<_centerPeak.getRho()<<","<<_centerPeak.getTheta()<<" circle center(x,y):("<<(1./_centerPeak.getRho())*cos(_centerPeak.getTheta())<<","<<(1./_centerPeak.getRho())*sin(_centerPeak.getTheta())<<")"<<endl;
1114 cout<<"par : "<<_d0<<","<<_omega<<","<<_phi0<<","<<_z0<<","<<_tanl<<", pt: "<<_pt3D<<endl;
1115 TrkHitList* qhits = p_trk->hits();
1116 TrkHotList::hot_iterator hotIter= qhits->hotList().begin();
1117 //if( typeid(*hotIter)!=typeid(MdcRecoHitOnTrack) && typeid(*hotIter)!=typeid(MdcHitOnTrack) && typeid(*hotIter)!=typeid(TrkHitOnTrk) )return;
1118 //int lay=((MdcHit*)(hotIter->hit()))->layernumber();
1119 //int wir=((MdcHit*)(hotIter->hit()))->wirenumber();
1120 while (hotIter!=qhits->hotList().end()) {
1121 if( typeid(*hotIter)==typeid(MdcRecoHitOnTrack) || typeid(*hotIter)==typeid(MdcHitOnTrack) /*|| typeid(*hotiter)!=typeid(trkhitontrk)*/ ){
1122 cout <<"(" <<((MdcHit*)(hotIter->hit()))->layernumber()
1123 <<","<<((MdcHit*)(hotIter->hit()))->wirenumber()
1124 <<":"<<hotIter->isActive()<<") ";
1125 cout << "nuse "<<hotIter->hit()->nUsedHits()<<endl;
1126 }else if(typeid(*hotIter)==typeid(CgemHitOnTrack)){
1127 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hotIter));
1128 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
1129 cout<<"(" <<recCgemCluster->getlayerid()
1130 <<","<<recCgemCluster->getsheetid()<<") "
1131 <<"clusterID:"<<recCgemCluster->getclusterid()
1132 <<" flag:"<<recCgemCluster->getflag()<<endl;
1133 }
1134 hotIter++;
1135 }
1136}
1137void HoughTrack::printHoughTrack()
1138{
1139 cout<<".........................................print HoughTrack "<<_trkid<<"............................................."<<endl;
1140 //cout<<"q d0 phi0 omega d0 tanl: "<<_charge<<" "<<_d0<<" "<<_phi0<<" "<<_omega<<" "<<_tanl<<" "<<_z0<<endl;
1141 //cout<<"q pt Xc Yc R:"<<_charge<<" "<<_ptLeast<<" "<<_centerX<<" "<<_centerY<<" "<<_centerR<<endl;
1142 cout<<"(Xc, Yc, R):"<<"("<<_centerX<<", "<<_centerY<<", "<<_centerR<<") (theta,rho):"<<"("<<_centerPeak.getTheta()<<", "<<_centerPeak.getRho()<<")"<<" ,pt_2D: "<<_pt2D*_charge<<" pt_3D: "<<_pt3D*_charge<<endl;
1143 cout<<endl;
1144}
1145void HoughTrack::printRecHit() const{
1146 cout<<"print rec hit"<<endl;
1147 double rho = _centerPeak.getRho();
1148 double theta= _centerPeak.getTheta();
1149 int size=_recHitVec.size();
1150 for(int i=0;i<size;i++){
1151 int layer= _recHitVec[i].getLayerId();
1152 int wire = _recHitVec[i].getWireId();
1153 int slant= _recHitVec[i].getSlayerType();
1154 int flag= _recHitVec[i].getflag();
1155 int style= _recHitVec[i].getStyle();
1156 int cirlist= _recHitVec[i].getCirList();
1157 if( slant ==0 ) std::cout<<"axial hit ("<<layer<<","<<wire<<") "<<_recHitVec[i].getDisToTrack()<<" "<<_recHitVec[i].getDisToCir()<<" "<<flag<<" "<<style<<" "<<cirlist<<std::endl;
1158 }
1159 for(int i=0;i<size;i++){
1160 int layer= _recHitVec[i].getLayerId();
1161 int wire = _recHitVec[i].getWireId();
1162 int slant= _recHitVec[i].getSlayerType();
1163 int flag= _recHitVec[i].getflag();
1164 int style= _recHitVec[i].getStyle();
1165 int cirlist= _recHitVec[i].getCirList();
1166 if( slant !=0 ) std::cout<<"stereo hit ("<<layer<<","<<wire<<") "<<_recHitVec[i].getDisToTrack()<<" "<<flag<<" "<<style<<" "<<cirlist<<std::endl;
1167 }
1168}
1169
time
Double_t time
Definition: DataBase/tau_mode.c:7
sin
double sin(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
digi_in_track
bool digi_in_track(const HoughRecHit &hita, const HoughRecHit &hitb)
Definition: MdcHoughFinder/MdcHoughFinder-00-00-12/src/HoughTrack.cxx:192
less_layer
bool less_layer(const int &a, const int &b)
Definition: MdcHoughFinder/MdcHoughFinder-00-00-12/src/HoughTrack.cxx:215
n2
int n2
Definition: SD0Tag.cxx:55
n1
int n1
Definition: SD0Tag.cxx:54
M_PI
#define M_PI
Definition: TConstant.h:4
Helix::dr
double dr(void) const
returns an element of parameters.
Definition: InstallArea/include/TrackUtil/TrackUtil/Helix.h:236
Hough2D::fit
virtual int fit()
Definition: Hough2D.cxx:140
Hough3D::print
void print()
Definition: Hough3D.cxx:332
Hough3D::fit
int fit()
Definition: Hough3D.cxx:80
HoughRecHit::setRet
void setRet(std::pair< double, double > theta_l)
Definition: InstallArea/include/MdcHoughFinder/MdcHoughFinder/HoughRecHit.h:29
HoughStereo::setRecHit
void setRecHit()
Definition: HoughStereo.cxx:315
HoughTrack::operator=
HoughTrack & operator=(const HoughTrack &other)
Definition: HoughTransAlg/HoughTransAlg-00-00-14/src/HoughTrack.cxx:254
HoughTrack::intersect_cylinder
double intersect_cylinder(double r_cylinder, double r_center, double phi_center, int charge)
Definition: MdcHoughFinder/MdcHoughFinder-00-00-12/src/HoughTrack.cxx:1030
TrkHitList::hotList
const TrkHotList & hotList() const
Definition: TrkHitList.cxx:51
TrkRecoTrk::id
const TrkId & id() const
Definition: TrkRecoTrk.cxx:134
EvtCyclic3::other
Index other(Index i, Index j)
Definition: EvtCyclic3.cc:118