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