CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
Classes | Public Member Functions | Static Public Member Functions | Static Public Attributes | List of all members
HoughZsFit Class Reference

#include <HoughZsFit.h>

Classes

struct  Line
 

Public Member Functions

 HoughZsFit (vector< HoughRecHit > *recHitCol)
 
 ~HoughZsFit ()
 
void leastFit (Line &linefit, int N)
 
int leastLine (int, double x[], double y[], double &, double &, double &)
 
void doit ()
 
void initPoint ()
 
double getTanl () const
 
double getZ0 () const
 
double getPro () const
 
void sortHit ()
 
void print ()
 
int cgemfit ()
 
void hough ()
 
 HoughZsFit (vector< HoughRecHit > *recHitCol)
 
 ~HoughZsFit ()
 
void leastFit (Line &linefit, int N)
 
int leastLine (int, double x[], double y[], double &, double &, double &)
 
void doit ()
 
void initPoint ()
 
double getTanl () const
 
double getZ0 () const
 
double getPro () const
 
void sortHit ()
 
void print ()
 
int cgemfit ()
 
void hough ()
 

Static Public Member Functions

static void setmethod (int m)
 
static void setmethod (int m)
 

Static Public Attributes

static int m_debug =0
 
static int m_recMethod =1
 
static double m_factor =5
 
static int m_nPoint3D =5
 
static vector< double > m_cut_deltaZ
 

Detailed Description

Definition at line 4 of file InstallArea/include/MdcHoughFinder/MdcHoughFinder/HoughZsFit.h.

Constructor & Destructor Documentation

◆ HoughZsFit() [1/2]

HoughZsFit::HoughZsFit ( vector< HoughRecHit > *  recHitCol)

Definition at line 35 of file HoughZsFit.cxx.

36{
37 vector<HoughRecHit>::iterator iter=recHitCol->begin();
38 for(;iter!=recHitCol->end();iter++){
39 //if( fabs((*iter).getzAmb(0))>10 && fabs((*iter).getzAmb(1))>10 ) continue;
40 if( (*iter).getDetectorType()==MDC && (*iter).getSlayerType() !=0 && (*iter).getflag()==0 ) _recStereoHit.push_back(&(*iter));
41 if( (*iter).getDetectorType()==CGEM && (*iter).getCluster()->getflag()==2 ) _recStereoHit.push_back(&(*iter));
42 //if( (*iter).getDetectorType()==CGEM && (*iter).getCluster()->getflag()==2 ) _cgemCluster[(*iter).getLayerId()].push_back(&(*iter));
43 }
44 _hitSize=_recStereoHit.size();
45 _pro_correct =0;
46 sortHit();
47 if(m_debug>0) print();
48 if(m_debug>0) cout<<"size of rechit in zs fit "<<_hitSize<<endl;
49 //cout<<"size of rechit in zs fit "<<_hitSize<<endl;
50
51// _vecPoint = new Point*[_hitSize];
52// for(int i=0;i<_hitSize;i++){
53// _vecPoint[i] = new Point[2];
54// }
55// initPoint();
56 if(m_recMethod==0){
57 if(_hitSize>=3) doit();
58 else cgemfit();
59 }else{
60 hough();
61 }
62 //{_tanl=0;_z0=0;}
63}
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:101
recHitCol
std::vector< HoughRecHit > recHitCol
Definition: InstallArea/include/MdcHoughFinder/MdcHoughFinder/Hough2D.h:24
HoughZsFit::doit
void doit()
Definition: HoughZsFit.cxx:68
HoughZsFit::hough
void hough()
Definition: HoughZsFit.cxx:588
HoughZsFit::print
void print()
Definition: HoughZsFit.cxx:454
HoughZsFit::cgemfit
int cgemfit()
Definition: HoughZsFit.cxx:460
HoughZsFit::sortHit
void sortHit()
Definition: HoughZsFit.cxx:451

◆ ~HoughZsFit() [1/2]

HoughZsFit::~HoughZsFit ( )

Definition at line 65 of file HoughZsFit.cxx.

65 {
66}

◆ HoughZsFit() [2/2]

HoughZsFit::HoughZsFit ( vector< HoughRecHit > *  recHitCol)

◆ ~HoughZsFit() [2/2]

HoughZsFit::~HoughZsFit ( )

Member Function Documentation

◆ cgemfit() [1/2]

int HoughZsFit::cgemfit ( )

Definition at line 460 of file HoughZsFit.cxx.

460 {
461 std::vector<int> fired_layer;
462 //---find the cgem layers that have at least one cluster
463 for(int ilayer=2;ilayer>=0;ilayer--){
464 int ncluster = _cgemCluster[ilayer].size();
465 if(ncluster!=0)fired_layer.push_back(ilayer);
466 }
467 if(fired_layer.size()<2){
468 _tanl=0;
469 _z0=0;
470 return -1;
471 }
472 //---find the correct cluster by s1/s2=z1/z2
473 double s1,s2,z1,z2;
474 double d1 =9999,d2=9999;
475 HoughRecHit *cluster1,*cluster2,*cluster3,*cluster4;
476 for(int i=0;i<_cgemCluster[fired_layer[0]].size();i++){
477 _cgemCluster[fired_layer[0]][i]->setflag(-999);
478 s1 = _cgemCluster[fired_layer[0]][i]->getsPos();
479 z1 = _cgemCluster[fired_layer[0]][i]->getzPos();
480 if(s1==0 || z1==0)continue;
481 for(int j=0;j<_cgemCluster[fired_layer[1]].size();j++){
482 _cgemCluster[fired_layer[1]][j]->setflag(-999);
483 s2 = _cgemCluster[fired_layer[1]][j]->getsPos();
484 z2 = _cgemCluster[fired_layer[1]][j]->getzPos();
485 if(s2==0 || z2==0)continue;
486 double delta=fabs(s1/s2-z1/z2);
487 if(delta<d1){
488 d1 = delta;
489 cluster1 = _cgemCluster[fired_layer[0]][i];
490 cluster2 = _cgemCluster[fired_layer[1]][j];
491 }
492 }
493 if(fired_layer.size()==3){
494 for(int j=0;j<_cgemCluster[fired_layer[2]].size();j++){
495 _cgemCluster[fired_layer[2]][j]->setflag(-999);
496 s2 = _cgemCluster[fired_layer[2]][j]->getsPos();
497 z2 = _cgemCluster[fired_layer[2]][j]->getzPos();
498 if(s2==0 || z2==0)continue;
499 double delta=fabs(s1/s2-z1/z2);
500 if(delta<d2){
501 d2 = delta;
502 cluster3 = _cgemCluster[fired_layer[0]][i];
503 cluster4 = _cgemCluster[fired_layer[2]][j];
504 }
505 }
506 }
507 }
508 Line linefit;
509 linefit._pointCol.push_back(*cluster1);
510 linefit._pointCol.push_back(*cluster2);
511 cluster1->setflag(0);
512 cluster2->setflag(0);
513 leastFit(linefit,2);
514 int line_size=2;
515 if(fired_layer.size()==3){
516 linefit._pointCol.push_back(*cluster4);
517 cluster4->setflag(0);
518 leastFit(linefit,3);
519 line_size++;
520 }
521
522 //---fit with ODC hits
523 for(int j=0;j<_hitSize;j++){
524 double chi_last=linefit._chi;
525 double k_last=linefit._k;
526 double b_last=linefit._b;
527 linefit._pointCol.push_back( *(_recStereoHit[j]) );
528 linefit._pointCol[line_size].setAmb(0);
529 if(m_debug>0) cout<<"Add point left "<<"("<<linefit._pointCol[line_size].getLayerId()<<","<<linefit._pointCol[line_size].getWireId()<<") "<<linefit._pointCol[line_size].getStyle()<<" "<<linefit._pointCol[line_size].getsPos()<<" "<<linefit._pointCol[line_size].getzPos()<<endl;
530 leastFit(linefit,line_size+1);
531 double chil=linefit._chi;
532 double kl=linefit._k;
533 double bl=linefit._b;
534 if (linefit._pointCol[line_size].getsPos()==-99) {
535 chil = 9999;
536 }
537 if(m_debug>0) cout<<"left "<<line_size<<" chi k b "<<chil<<" "<<kl<<" "<<bl<<endl;
538
539 linefit._pointCol[line_size].setAmb(1);
540 if(m_debug>0) cout<<"Add point right "<<"("<<linefit._pointCol[line_size].getLayerId()<<","<<linefit._pointCol[line_size].getWireId()<<") "<<linefit._pointCol[line_size].getStyle()<<" "<<linefit._pointCol[line_size].getsPos()<<" "<<linefit._pointCol[line_size].getzPos()<<endl;
541 leastFit(linefit,line_size+1);
542 double chir=linefit._chi;
543 double kr=linefit._k;
544 double br=linefit._b;
545 if (linefit._pointCol[line_size].getsPos()==-99) {
546 chir = 9999;
547 }
548 if(m_debug>0) cout<<"right "<<line_size<<" chi k b "<<chir<<" "<<kr<<" "<<br<<endl;
549
550 if(chil<chir) {
551 linefit._pointCol[line_size].setAmb(0);
552 linefit._chi=chil;
553 linefit._k=kl;
554 linefit._b=bl;
555 linefit._ambig.push_back(0);
556 _recStereoHit[j]->setAmb(0);
557 }
558 else{
559 linefit._ambig.push_back(1);
560 _recStereoHit[j]->setAmb(1);
561 }
562 line_size++;
563
564 double dChi = chi_last-linefit._chi;
565 double dChi_n = (linefit._chi)/(line_size-1)-(chi_last/line_size);
566 if(m_debug>0) cout<<"dChi: "<<dChi<<endl;
567 if(m_debug>0) cout<<"dChi/n: "<<dChi_n<<endl;
568 if(m_debug>0) cout<<endl;
569 // if( fabs(dChi) > 500.) //FIXME
570 if( fabs(dChi_n) > 25.) //FIXME
571 {
572 _recStereoHit[j]->setAmb(-999);
573 _recStereoHit[j]->setflag(-999);
574 linefit._pointCol.pop_back();
575 linefit._chi=chi_last;
576 linefit._k=k_last;
577 linefit._b=b_last;
578 linefit._ambig.at(j)=-999;
579 line_size--;
580 }
581 }
582 _tanl=linefit._k;
583 _z0=linefit._b;
584 if(m_debug>0) cout<<"z0 tanl : "<<_z0<<" "<<_tanl<<endl;
585 return 0;
586}
HoughZsFit::leastFit
void leastFit(Line &linefit, int N)
Definition: HoughZsFit.cxx:395

Referenced by HoughZsFit().

◆ cgemfit() [2/2]

int HoughZsFit::cgemfit ( )

◆ doit() [1/2]

void HoughZsFit::doit ( )

Definition at line 68 of file HoughZsFit.cxx.

68 {
69
70// Point point0;
71// point0.first=0.;
72// point0.second=0.;
73
74 Line linefit[8];
75
76 for(int iline=0;iline<8;iline++){
77 linefit[iline]._pointCol.push_back( *(_recStereoHit[0]) );
78 linefit[iline]._pointCol.push_back( *(_recStereoHit[1]) );
79 linefit[iline]._pointCol.push_back( *(_recStereoHit[2]) );
80 }
81 if(m_debug>0) {
82 for(int i =0;i<3;i++){
83 cout<<" the first 3 hits ("<<i<<" "<< _recStereoHit[i]->getLayerId()<<","<<_recStereoHit[i]->getWireId()<<")"<<endl;
84 cout<<" left "<<(_recStereoHit[i])->getsAmb(0)<<" "<<(_recStereoHit[i])->getzAmb(0)<<endl;
85 cout<<" right "<<(_recStereoHit[i])->getsAmb(1)<<" "<<(_recStereoHit[i])->getzAmb(1)<<endl;
86 }
87 }
88
89 linefit[0]._pointCol[0].setAmb(0);
90 linefit[0]._pointCol[1].setAmb(0);
91 linefit[0]._pointCol[2].setAmb(0);
92
93 linefit[1]._pointCol[0].setAmb(0);
94 linefit[1]._pointCol[1].setAmb(0);
95 linefit[1]._pointCol[2].setAmb(1);
96
97 linefit[2]._pointCol[0].setAmb(0);
98 linefit[2]._pointCol[1].setAmb(1);
99 linefit[2]._pointCol[2].setAmb(0);
100
101 linefit[3]._pointCol[0].setAmb(0);
102 linefit[3]._pointCol[1].setAmb(1);
103 linefit[3]._pointCol[2].setAmb(1);
104
105
106 linefit[4]._pointCol[0].setAmb(1);
107 linefit[4]._pointCol[1].setAmb(0);
108 linefit[4]._pointCol[2].setAmb(0);
109
110 linefit[5]._pointCol[0].setAmb(1);
111 linefit[5]._pointCol[1].setAmb(0);
112 linefit[5]._pointCol[2].setAmb(1);
113
114 linefit[6]._pointCol[0].setAmb(1);
115 linefit[6]._pointCol[1].setAmb(1);
116 linefit[6]._pointCol[2].setAmb(0);
117
118 linefit[7]._pointCol[0].setAmb(1);
119 linefit[7]._pointCol[1].setAmb(1);
120 linefit[7]._pointCol[2].setAmb(1);
121
122 for(int i=0;i<8;i++){
123 linefit[i]._ambig.push_back( linefit[i]._pointCol[0].getAmbig() );
124 linefit[i]._ambig.push_back( linefit[i]._pointCol[1].getAmbig() );
125 linefit[i]._ambig.push_back( linefit[i]._pointCol[2].getAmbig() );
126 if(m_debug>0) cout<<" ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^begin line "<<i<<endl;
127 leastFit(linefit[i],3);
128 //judst if there is noise in the 3 outer
129 int line_size=3;
130 for(int j=3;j<_hitSize;j++){
131 double chi_last=linefit[i]._chi;
132 double k_last=linefit[i]._k;
133 double b_last=linefit[i]._b;
134 if(m_debug>0) cout<<"last "<<j<<" chi k b "<<chi_last<<" "<<k_last<<" "<<b_last<<endl;
135 linefit[i]._pointCol.push_back( *(_recStereoHit[j]) );
136 //int layer =linefit[i]._pointCol.back().getLayerId();
137 //int wire=linefit[i]._pointCol.back().getWireId();
138 //cout<<"???ambig ("<<layer<<","<<wire<<") "<<linefit[i]._pointCol.back().getLrTruth()<<endl;
139 linefit[i]._pointCol[line_size].setAmb(0);
140 if(m_debug>0) cout<<"Add point left "<<"("<<linefit[i]._pointCol[line_size].getLayerId()<<","<<linefit[i]._pointCol[line_size].getWireId()<<") "<<linefit[i]._pointCol[line_size].getStyle()<<" "<<linefit[i]._pointCol[line_size].getsPos()<<" "<<linefit[i]._pointCol[line_size].getzPos()<<endl;
141 leastFit(linefit[i],line_size+1);
142 double chil=linefit[i]._chi;
143 double kl=linefit[i]._k;
144 double bl=linefit[i]._b;
145 if (linefit[i]._pointCol[line_size].getsPos()==-99) {
146 chil = 9999;
147 }
148 if(m_debug>0) cout<<"left "<<line_size<<" chi k b "<<chil<<" "<<kl<<" "<<bl<<endl;
149
150 linefit[i]._pointCol[line_size].setAmb(1);
151 if(m_debug>0) cout<<"Add point right "<<"("<<linefit[i]._pointCol[line_size].getLayerId()<<","<<linefit[i]._pointCol[line_size].getWireId()<<") "<<linefit[i]._pointCol[line_size].getStyle()<<" "<<linefit[i]._pointCol[line_size].getsPos()<<" "<<linefit[i]._pointCol[line_size].getzPos()<<endl;
152 leastFit(linefit[i],line_size+1);
153 double chir=linefit[i]._chi;
154 double kr=linefit[i]._k;
155 double br=linefit[i]._b;
156 if (linefit[i]._pointCol[line_size].getsPos()==-99) {
157 chir = 9999;
158 }
159 if(m_debug>0) cout<<"right "<<line_size<<" chi k b "<<chir<<" "<<kr<<" "<<br<<endl;
160
161 if(chil<chir) {
162 linefit[i]._pointCol[line_size].setAmb(0);
163 linefit[i]._chi=chil;
164 linefit[i]._k=kl;
165 linefit[i]._b=bl;
166 linefit[i]._ambig.push_back(0);
167 //_recStereoHit[j]->setAmb(0);
168 }
169 else linefit[i]._ambig.push_back(1);
170 line_size++;
171
172 //test ambig with truth
173 int same_ambig = 1;
174 for(int ihit=0;ihit<line_size;ihit++){
175 int ambighit = linefit[i]._pointCol[ihit].getAmbig();
176 int ambigTruth = linefit[i]._pointCol[ihit].getLrTruth();
177 int layer =linefit[i]._pointCol[ihit].getLayerId();
178 int wire=linefit[i]._pointCol[ihit].getWireId();
179 //if (ambigTruth == -1) ambigTruth=1;
180 //else if (ambigTruth == 1) ambigTruth=0;
181 //else cout << "wrong in ambig Truth "<<endl;
182
183 //if(ambighit!= ambigTruth ) same_ambig=0;
184 }
185 //if(m_debug >0 ) cout << "same ambig ? "<<same_ambig<<endl;
186
187
188 double dChi = chi_last-linefit[i]._chi;
189 double dChi_n = (linefit[i]._chi)/(line_size)-(chi_last/line_size-1);
190 if(m_debug>0) cout<<"Chi: "<<linefit[i]._chi<<endl;
191 if(m_debug>0) cout<<"dChi: "<<dChi<<endl;
192 if(m_debug>0) cout<<"dChi/n: "<<dChi_n<<endl;
193 if(m_debug>0) cout<<endl;
194 // if( fabs(dChi) > 500.) //FIXME
195 if( fabs(dChi_n) > 25.) //FIXME
196 // linefit[i]._pointCol[j].setAmb(-999);
197 //_recStereoHit[j]->setAmb(-999);
198 //_recStereoHit[j]->setflag(-999);
199 {
200 linefit[i]._pointCol.pop_back();
201 linefit[i]._chi=chi_last;
202 linefit[i]._k=k_last;
203 linefit[i]._b=b_last;
204 linefit[i]._ambig.at(j)=-999;
205 line_size--;
206 }
207 //cout<<linefit[i]._pointCol.size()<<" "<<linefit[i]._ambig.size()<<endl;
208 }
209
210 bool find_best_cluster = false;
211 for(int ilayer=2;ilayer>=0;ilayer--){
212 int ncluster = _cgemCluster[ilayer].size();
213 if(ncluster==0)continue;
214 double chi_last=linefit[i]._chi;
215 double k_last=linefit[i]._k;
216 double b_last=linefit[i]._b;
217 if(m_debug>0) cout<<"last chi k b "<<chi_last<<" "<<k_last<<" "<<b_last<<endl;
218 double chi_temp=9999;
219 bool find_best_cluster = false;
220 HoughRecHit* cluster_temp;
221 for(int icluster=0;icluster<ncluster;icluster++){
222 _cgemCluster[ilayer][icluster]->setflag(-999);
223 Line linefit_temp = linefit[i];
224 linefit_temp._pointCol.push_back(*(_cgemCluster[ilayer][icluster]));
225 leastFit(linefit_temp,line_size+1);
226 if(linefit_temp._chi<chi_temp){
227 chi_temp = linefit_temp._chi;
228 cluster_temp = _cgemCluster[ilayer][icluster];
229 find_best_cluster = true;
230 }
231 }
232 if(!find_best_cluster)continue;
233 linefit[i]._pointCol.push_back(*cluster_temp);
234 linefit[i]._clusterCol.push_back(cluster_temp);
235 leastFit(linefit[i],line_size+1);
236 line_size++;
237
238 double dChi = chi_last-linefit[i]._chi;
239 double dChi_n = (linefit[i]._chi)/(line_size-1)-(chi_last/line_size);
240 if(m_debug>0) cout<<"Chi: "<<linefit[i]._chi<<endl;
241 if(m_debug>0) cout<<"dChi: "<<dChi<<endl;
242 if(m_debug>0) cout<<"dChi/n: "<<dChi_n<<endl;
243 if(m_debug>0) cout<<endl;
244 if( fabs(dChi_n) > 25.) //FIXME
245 {
246 linefit[i]._pointCol.pop_back();
247 linefit[i]._clusterCol.pop_back();
248 linefit[i]._chi=chi_last;
249 linefit[i]._k=k_last;
250 linefit[i]._b=b_last;
251 //linefit[i]._ambig.at(j)=-999;
252 line_size--;
253 }
254 }
255 //reject wrong line
256 if( (linefit[i]._pointCol[0].getsPos()==-99) || (linefit[i]._pointCol[1].getsPos()==-99) || (linefit[i]._pointCol[2].getsPos()==-99) ) {
257 linefit[i]._chi=99999;
258 }
259 }
260 std::sort(linefit,linefit+8,compare_zsfit);
261 for(int i=0;i<8;i++){
262 if(m_debug>0) cout<<"Line :"<<i<<" chis: "<<linefit[i]._chi<<" k,b: "<<linefit[i]._k<<" "<<linefit[i]._b<<endl;
263 int ambig_correct = 0;
264 for(int j=0;j<_hitSize;j++){
265 int ambig = linefit[i]._ambig.at(j);
266 int layer= _recStereoHit.at(j)->getLayerId();
267 int wire= _recStereoHit.at(j)->getWireId();
268 int style= _recStereoHit.at(j)->getStyle();
269 double l=-99;
270 if (ambig!=-999) l= _recStereoHit.at(j)->getsAmb(ambig);
271 double z=-99;
272 if (ambig!=-999) z= _recStereoHit.at(j)->getzAmb(ambig);
273 // cout<<"debug ("<<layer<<" ,"<<wire<<") style "<<style<<" ambig "<<ambig<<" s "<<l<<" z "<<z<<endl;
274 if (l==-99 && z==-99) ambig=-999;
275 if(m_debug>0) cout<<"("<<layer<<" ,"<<wire<<") style "<<style<<" ambig "<<ambig<<" s "<<l<<" z "<<z<<endl;
276 if(i==0) {
277 _recStereoHit[j]->setAmb(ambig); // when line 0 select
278 if( ambig ==-999) _recStereoHit[j]->setflag(-999); // when line 0 select
279 int ambigTruth = _recStereoHit.at(j)->getLrTruth();
280 if (ambigTruth == -1) ambigTruth=1;
281 else if (ambigTruth == 1) ambigTruth=0;
282 if(ambig ==ambigTruth) ambig_correct++;
283 //cout<<"ambig : "<<ambig<<" "<<ambigTruth<<endl;
284 _pro_correct = (double)ambig_correct/(double)_hitSize;
285
286 for(int icluster=0;icluster<linefit[i]._clusterCol.size();icluster++){
287 linefit[i]._clusterCol[icluster]->setflag(0);
288 }
289 }
290 }
291 }
292 _tanl=linefit[0]._k;
293 _z0=linefit[0]._b;
294 if(m_debug>0) cout<<"z0 tanl : "<<_z0<<" "<<_tanl<<endl;
295
296 /*
297 linefit[0]._pointCol.push_back(_vecPoint[_hitSize-1][0]);
298 linefit[0]._pointCol.push_back(_vecPoint[_hitSize-2][0]);
299 linefit[0]._pointCol.push_back(_vecPoint[_hitSize-3][0]);
300
301 linefit[1]._pointCol.push_back(_vecPoint[_hitSize-1][0]);
302 linefit[1]._pointCol.push_back(_vecPoint[_hitSize-2][0]);
303 linefit[1]._pointCol.push_back(_vecPoint[_hitSize-3][1]);
304
305 linefit[2]._pointCol.push_back(_vecPoint[_hitSize-1][0]);
306 linefit[2]._pointCol.push_back(_vecPoint[_hitSize-2][1]);
307 linefit[2]._pointCol.push_back(_vecPoint[_hitSize-3][0]);
308
309 linefit[3]._pointCol.push_back(_vecPoint[_hitSize-1][0]);
310 linefit[3]._pointCol.push_back(_vecPoint[_hitSize-2][1]);
311 linefit[3]._pointCol.push_back(_vecPoint[_hitSize-3][1]);
312
313 linefit[4]._pointCol.push_back(_vecPoint[_hitSize-1][1]);
314 linefit[4]._pointCol.push_back(_vecPoint[_hitSize-2][0]);
315 linefit[4]._pointCol.push_back(_vecPoint[_hitSize-3][0]);
316
317 linefit[5]._pointCol.push_back(_vecPoint[_hitSize-1][1]);
318 linefit[5]._pointCol.push_back(_vecPoint[_hitSize-2][0]);
319 linefit[5]._pointCol.push_back(_vecPoint[_hitSize-3][1]);
320
321 linefit[6]._pointCol.push_back(_vecPoint[_hitSize-1][1]);
322 linefit[6]._pointCol.push_back(_vecPoint[_hitSize-2][1]);
323 linefit[6]._pointCol.push_back(_vecPoint[_hitSize-3][0]);
324
325 linefit[7]._pointCol.push_back(_vecPoint[_hitSize-1][1]);
326 linefit[7]._pointCol.push_back(_vecPoint[_hitSize-2][1]);
327 linefit[7]._pointCol.push_back(_vecPoint[_hitSize-3][1]);
328
329
330 for(int i=0;i<8;i++){
331 linefit[i]._amg=i*pow(10,(_hitSize-3));
332 if(m_debug>0) cout<<"Start line ooooooooooooooooooooooooooooooooooooooooooooooooo"<<i<<endl;
333 for(int j=4;j<_hitSize+1;j++){
334 if(m_debug>0) cout<<"Start add _vecPoint ****************************"<<j<<endl;
335 linefit[i]._pointCol.push_back(_vecPoint[_hitSize-j][0]);
336 leastFit(linefit[i],j);
337 double chil=linefit[i]._chi;
338 double kl=linefit[i]._k;
339 double bl=linefit[i]._b;
340
341 linefit[i]._pointCol[j-1]=_vecPoint[_hitSize-j][1];
342 leastFit(linefit[i],j);
343 double chir=linefit[i]._chi;
344 double kr=linefit[i]._k;
345 double br=linefit[i]._b;
346 if(m_debug>0){
347 cout<<chil<<" "<<kl<<" "<<bl<<endl;
348 cout<<chir<<" "<<kr<<" "<<br<<endl;
349 }
350
351 if(chil<chir) {
352 // cout<<"j: "<<j<<endl;
353 linefit[i]._pointCol[j-1]=_vecPoint[_hitSize-j][0];
354 linefit[i]._chi=chil;
355 linefit[i]._k=kl;
356 linefit[i]._b=bl;
357 linefit[i]._amg+=0*pow(10,(_hitSize-j));
358 }
359 else { linefit[i]._amg+=1*pow(10,(_hitSize-j)); }
360 }
361 if(m_debug>0) cout<<"Start add _vecPoint ****************************"<<"0"<<endl;
362 linefit[i]._pointCol.push_back(point0);
363 leastFit(linefit[i],_hitSize+1);
364 double chi0=linefit[i]._chi;
365 double k0=linefit[i]._k;
366 double b0=linefit[i]._b;
367 if(m_debug>0) cout<<chi0<<" "<<k0<<" "<<b0<<endl;
368}
369std::sort(linefit,linefit+8,compare_zsfit);
370for(int i=0;i<8;i++){
371 if(m_debug>0) cout<<"Line :"<<i<<" chis: "<<linefit[i]._chi<<" k,b: "<<linefit[i]._k<<" "<<linefit[i]._b<<" amg: "<<linefit[i]._amg<<endl;
372}
373_tanl=linefit[0]._k;
374_z0=linefit[0]._b;
375if(m_debug>0) cout<<"z0 tanl : "<<_z0<<" "<<_tanl<<endl;
376*/
377}
compare_zsfit
bool compare_zsfit(const HoughZsFit::Line &a, const HoughZsFit::Line &b)
Definition: HoughZsFit.cxx:416

Referenced by HoughZsFit().

◆ doit() [2/2]

void HoughZsFit::doit ( )

◆ getPro() [1/2]

double HoughZsFit::getPro ( ) const
inline

Definition at line 25 of file InstallArea/include/MdcHoughFinder/MdcHoughFinder/HoughZsFit.h.

25{return _pro_correct;}

Referenced by HoughTrack::fitzs().

◆ getPro() [2/2]

double HoughZsFit::getPro ( ) const
inline

Definition at line 25 of file Reconstruction/MdcHoughFinder/MdcHoughFinder-00-00-12/MdcHoughFinder/HoughZsFit.h.

25{return _pro_correct;}

◆ getTanl() [1/2]

double HoughZsFit::getTanl ( ) const
inline

Definition at line 23 of file InstallArea/include/MdcHoughFinder/MdcHoughFinder/HoughZsFit.h.

23{return _tanl;}

Referenced by HoughTrack::fitzs().

◆ getTanl() [2/2]

double HoughZsFit::getTanl ( ) const
inline

Definition at line 23 of file Reconstruction/MdcHoughFinder/MdcHoughFinder-00-00-12/MdcHoughFinder/HoughZsFit.h.

23{return _tanl;}

◆ getZ0() [1/2]

double HoughZsFit::getZ0 ( ) const
inline

Definition at line 24 of file InstallArea/include/MdcHoughFinder/MdcHoughFinder/HoughZsFit.h.

24{return _z0;}

Referenced by HoughTrack::fitzs().

◆ getZ0() [2/2]

double HoughZsFit::getZ0 ( ) const
inline

Definition at line 24 of file Reconstruction/MdcHoughFinder/MdcHoughFinder-00-00-12/MdcHoughFinder/HoughZsFit.h.

24{return _z0;}

◆ hough() [1/2]

void HoughZsFit::hough ( )

Definition at line 588 of file HoughZsFit.cxx.

588 {
589 int extNbin =10;
590 int binx,biny,binz;
591 double x_max = 0.93/sqrt(1-0.93*0.93);
592 double y_max = 50;
593 int x_bin = 800 ;
594 int y_bin = 800 ;
595 //int m_nPoint3D = 5;
596 double x_step = 2*x_max/x_bin/m_nPoint3D;
597 //cout<<m_nPoint3D<<endl;
598 //double y_step = 2*y_max/y_bin;
599 TH2D* houghSpace = new TH2D("houghSpace","houghSpace",x_bin,-x_max,x_max,y_bin,-y_max,y_max);
600 //TH2D* houghSpace = new TH2D("houghSpace","houghSpace",x_bin,-x_max+x_step/2,x_max+x_step/2,y_bin,-y_max+x_step/2,y_max+y_step/2);
601 for(int i=0;i<_hitSize;i++){
602 //cout<<_recStereoHit[i]->getLayerId()<<endl;
603 double x = -x_max - 0.5*x_step;
604 _recStereoHit[i]->setAmb(0);
605 //cout<<"l: "<<_recStereoHit[i]->getsPos()<<" "<<_recStereoHit[i]->getzPos()<<endl;
606 double yl = -_recStereoHit[i]->getsPos()*x + _recStereoHit[i]->getzPos();
607 double dyl = -_recStereoHit[i]->getsPos()*x_step;
608 _recStereoHit[i]->setAmb(1);
609 //cout<<"r: "<<_recStereoHit[i]->getsPos()<<" "<<_recStereoHit[i]->getzPos()<<endl;
610 double yr = -_recStereoHit[i]->getsPos()*x + _recStereoHit[i]->getzPos();
611 double dyr = -_recStereoHit[i]->getsPos()*x_step;
612 for(int j=0;j<x_bin*m_nPoint3D;j++){
613 x += x_step;
614 yl += dyl;
615 yr += dyr;
616 houghSpace->Fill(x,yl);
617 //if(_recStereoHit[i]->getDetectorType()==MDC)houghSpace->Fill(x,yr);
618 houghSpace->Fill(x,yr);
619 }
620 }
621 houghSpace->GetMaximumBin(binx,biny,binz);
622 _tanl = houghSpace->GetXaxis()->GetBinCenter(binx);
623 _z0 = houghSpace->GetYaxis()->GetBinCenter(biny);
624 //cout<<"map1 "<<_z0<<" "<<_tanl<<" "<<endl;
625/*
626 //double h = houghSpace->GetBinContent(binx,biny);
627 double x_low = houghSpace->GetXaxis()->GetBinLowEdge(binx-extNbin);
628 double x_up = houghSpace->GetXaxis()->GetBinUpEdge(binx +extNbin);
629 double y_low = houghSpace->GetYaxis()->GetBinLowEdge(biny-extNbin);
630 double y_up = houghSpace->GetYaxis()->GetBinUpEdge(biny +extNbin);
631 int x_bin2 = 100 ;
632 int y_bin2 = 100 ;
633 double x_step2 = (x_up-x_low)/x_bin2;
634 TH2D* houghSpace2 = new TH2D("houghSpace2","houghSpace2",x_bin2,x_low,x_up,y_bin2,y_low,y_up);
635 for(int i=0;i<_hitSize;i++){
636 double x = x_low - 0.5*x_step2;
637 _recStereoHit[i]->setAmb(0);
638 double yl = -_recStereoHit[i]->getsPos()*x + _recStereoHit[i]->getzPos();
639 double dyl = -_recStereoHit[i]->getsPos()*x_step2;
640 _recStereoHit[i]->setAmb(1);
641 double yr = -_recStereoHit[i]->getsPos()*x + _recStereoHit[i]->getzPos();
642 double dyr = -_recStereoHit[i]->getsPos()*x_step2;
643 for(int j=0;j<x_bin2;j++){
644 x += x_step2;
645 yl += dyl;
646 yr += dyr;
647 houghSpace2->Fill(x,yl);
648 if(_recStereoHit[i]->getDetectorType()==MDC)houghSpace2->Fill(x,yr);
649 //houghSpace2->Fill(x,yr);
650 }
651 }
652 houghSpace2->GetMaximumBin(binx,biny,binz);
653 _tanl = houghSpace2->GetXaxis()->GetBinCenter(binx);
654 _z0 = houghSpace2->GetYaxis()->GetBinCenter(biny);
655*/
656 int l[3] = {-1,-1,-1};//index of best cluster for each layer
657 double dZ[3] = {999,999,999};
658 int nStereo(0);
659 //cout<<nStereo<<endl;
660 for(int i=0;i<_hitSize;i++){
661 _recStereoHit[i]->setAmb(0);
662 double s0 = _recStereoHit[i]->getsPos();
663 double z0 = _recStereoHit[i]->getzPos();
664 _recStereoHit[i]->setAmb(1);
665 double s1 = _recStereoHit[i]->getsPos();
666 double z1 = _recStereoHit[i]->getzPos();
667 double zl = s0*_tanl + _z0;
668 double zr = s1*_tanl + _z0;
669 double z,s;
670 if(fabs(z0-zl) < fabs(z1-zr)){
671 s = s0;
672 z = z0;
673 _recStereoHit[i]->setAmb(0);
674 }else{
675 s = s1;
676 z = z1;
677 _recStereoHit[i]->setAmb(1);
678 }
679 double deltaZ = (z-(s*_tanl+_z0));
680 //if(fabs(deltaZ)<m_factor*z_cut[_recStereoHit[i]->getLayerId()]){
681 if(fabs(deltaZ)<m_factor*m_cut_deltaZ[_recStereoHit[i]->getLayerId()]){
682 _recStereoHit[i]->setflag(0);
683 nStereo++;
684 }else{
685 _recStereoHit[i]->setflag(-999);
686 }
687 if(_recStereoHit[i]->getDetectorType()==CGEM && _recStereoHit[i]->getflag() ==0){
688 if(fabs(deltaZ) < dZ[_recStereoHit[i]->getLayerId()]){
689 l[_recStereoHit[i]->getLayerId()] = i;
690 dZ[_recStereoHit[i]->getLayerId()] = fabs(deltaZ);
691 }
692 }
693 }
694 for(int i=0;i<_hitSize;i++){
695 if(_recStereoHit[i]->getDetectorType()==CGEM){
696 if(i == l[0] || i == l[1] || i== l[2]){
697 _recStereoHit[i]->setflag(0);
698 nStereo++;
699 }else{
700 _recStereoHit[i]->setflag(-999);
701 }
702 }
703 }
704 //cout<<"map2 "<<_z0<<" "<<_tanl<<" "<<nStereo<<endl;
705 delete houghSpace;
706 //delete houghSpace2;
707}
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
s
XmlRpcServer s
Definition: HelloServer.cpp:11

Referenced by HoughZsFit().

◆ hough() [2/2]

void HoughZsFit::hough ( )

◆ initPoint() [1/2]

void HoughZsFit::initPoint ( )

◆ initPoint() [2/2]

void HoughZsFit::initPoint ( )

◆ leastFit() [1/2]

void HoughZsFit::leastFit ( Line &  linefit,
int  N 
)

Definition at line 395 of file HoughZsFit.cxx.

395 {
396 double *x=new double[N+1];
397 double *y=new double[N+1];
398 for(int i=0;i<N;i++){
399 double a=linefit._pointCol[i].getsPos();
400 double b=linefit._pointCol[i].getzPos();
401 x[i]=a;
402 y[i]=b;
403 //if(m_debug>0) cout<<" x "<<a<<" y "<<b<<endl;
404 }
405 x[N]=0.;
406 y[N]=0.;
407 double k(0),b(0),chi2(0);
408 leastLine(N+1,x,y,k,b,chi2);
409 linefit._k =k;
410 linefit._b =b;
411 linefit._chi=chi2;
412 delete []x;
413 delete []y;
414
415}
HoughZsFit::leastLine
int leastLine(int, double x[], double y[], double &, double &, double &)
Definition: HoughZsFit.cxx:420

Referenced by cgemfit(), and doit().

◆ leastFit() [2/2]

void HoughZsFit::leastFit ( Line &  linefit,
int  N 
)

◆ leastLine() [1/2]

int HoughZsFit::leastLine ( int  nhit,
double  x[],
double  y[],
double &  k,
double &  b,
double &  chi2 
)

Definition at line 420 of file HoughZsFit.cxx.

420 {
421 double x_sum=0;
422 double y_sum=0;
423 double x2_sum=0;
424 double y2_sum=0;
425 double xy_sum=0;
426 for(int i=0;i<nhit;i++){
427 x_sum=x_sum+x[i];
428 y_sum=y_sum+y[i];
429 x2_sum=x2_sum+x[i]*x[i];
430 y2_sum=y2_sum+y[i]*y[i];
431 xy_sum=xy_sum+x[i]*y[i];
432 }
433 b=(x2_sum*y_sum-xy_sum*x_sum)/(nhit*x2_sum-x_sum*x_sum);
434 k=(nhit*xy_sum-x_sum*y_sum)/(nhit*x2_sum-x_sum*x_sum);
435
436 double yE[3];
437 for(int i=0;i<nhit;i++){
438 yE[i]=k*x[i]+b;
439 double chi2_temp;
440 if(yE[i]!=0) chi2_temp=(y[i]-yE[i])*(y[i]-yE[i])/yE[i]*yE[i];
441 else chi2_temp=0;
442 chi2+=chi2_temp;
443 }
444
445 return 1;
446}

Referenced by leastFit().

◆ leastLine() [2/2]

int HoughZsFit::leastLine ( int  ,
double  x[],
double  y[],
double &  ,
double &  ,
double &   
)

◆ print() [1/2]

void HoughZsFit::print ( )

Definition at line 454 of file HoughZsFit.cxx.

454 {
455 for(int i =0;i<_recStereoHit.size();i++){
456 cout<<"("<<_recStereoHit[i]->getLayerId()<<","<<_recStereoHit[i]->getWireId()<<")"<<endl;
457 }
458}

Referenced by HoughZsFit().

◆ print() [2/2]

void HoughZsFit::print ( )

◆ setmethod() [1/2]

static void HoughZsFit::setmethod ( int  m)
inlinestatic

Definition at line 34 of file InstallArea/include/MdcHoughFinder/MdcHoughFinder/HoughZsFit.h.

34{m_recMethod = m;}

Referenced by MdcHoughFinder::initialize().

◆ setmethod() [2/2]

static void HoughZsFit::setmethod ( int  m)
inlinestatic

Definition at line 34 of file Reconstruction/MdcHoughFinder/MdcHoughFinder-00-00-12/MdcHoughFinder/HoughZsFit.h.

34{m_recMethod = m;}

◆ sortHit() [1/2]

void HoughZsFit::sortHit ( )

Definition at line 451 of file HoughZsFit.cxx.

451 {
452 std::sort (_recStereoHit.begin(),_recStereoHit.end(),layer_in_track);
453}
layer_in_track
bool layer_in_track(const HoughRecHit *hita, const HoughRecHit *hitb)
Definition: HoughZsFit.cxx:448

Referenced by HoughZsFit().

◆ sortHit() [2/2]

void HoughZsFit::sortHit ( )

Member Data Documentation

◆ m_cut_deltaZ

vector< double > HoughZsFit::m_cut_deltaZ
static

Definition at line 35 of file InstallArea/include/MdcHoughFinder/MdcHoughFinder/HoughZsFit.h.

Referenced by hough(), and MdcHoughFinder::initialize().

◆ m_debug

int HoughZsFit::m_debug =0
static

Definition at line 28 of file InstallArea/include/MdcHoughFinder/MdcHoughFinder/HoughZsFit.h.

Referenced by cgemfit(), doit(), HoughZsFit(), and MdcHoughFinder::initialize().

◆ m_factor

double HoughZsFit::m_factor =5
static

Definition at line 30 of file InstallArea/include/MdcHoughFinder/MdcHoughFinder/HoughZsFit.h.

Referenced by hough(), and MdcHoughFinder::initialize().

◆ m_nPoint3D

int HoughZsFit::m_nPoint3D =5
static

Definition at line 31 of file InstallArea/include/MdcHoughFinder/MdcHoughFinder/HoughZsFit.h.

Referenced by hough(), and MdcHoughFinder::initialize().

◆ m_recMethod

int HoughZsFit::m_recMethod =1
static

Definition at line 29 of file InstallArea/include/MdcHoughFinder/MdcHoughFinder/HoughZsFit.h.

Referenced by HoughZsFit(), and setmethod().

The documentation for this class was generated from the following files: