CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
HoughZsFit.cxx
Go to the documentation of this file.
1#include "MdcHoughFinder/HoughZsFit.h"
2#include "TTree.h"
3#include "TH2D.h"
4#include "TF1.h"
5#include "TGraph.h"
6#include <algorithm>
7#include <assert.h>
8
9int HoughZsFit::m_debug=0;
10int HoughZsFit::m_recMethod =1;
11double HoughZsFit::m_factor=5;
12int HoughZsFit::m_nPoint3D=5;
13//double z_cut[43] = {//200MeV
14 ////0.255, 0.405, 0.455,//0.995
15 ////0.355, 0.485, 0.495,//0.999
16 //0.0354506, 0.0344371, 0.047115,
17 //0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
18 ////9.66, 10.65, 11.04, 12.41, 12.42, 14.25, 14.6, 15.75, 16.5, 17.25, 18.96, 20.25, 21.58, 22.41, 23.8, 24.65//0.995
19 ////12.74, 14.25, 13.92, 15.13, 15.66, 17.29, 18.2, 19.53, 20.9, 21.39, 22.8, 23.25, 25.22, 26.19, 27.16, 28.13//0.999
20 //0.930282, 1.07567, 1.1647, 1.23004, 1.45074, 1.42505, 1.46573, 1.46672, 1.66033, 1.74953, 1.65215, 1.74229, 1.58215, 1.69198, 1.73003, 1.78065,
21 //0,0,0,0,0,0,0
22 //};
23vector<double> HoughZsFit::m_cut_deltaZ(43,0);
24
25//double z_cut[36] = {300MeV
26 //0.145, 0.255, 0.345,//0.990
27 //0.185, 0.305, 0.405,//0.995
28 //0.285, 0.435, 0.485,//0.999
29 //0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
30 //7.42, 7.95, 7.84, 8.67, 9.54, 10.07, 11, 11.97, 12.54, 13.11, 13.68, 15.25, 16.38, 17.55, 18.2, 18.85//0.990
31 //9.1, 9.45, 10.08, 10.37, 10.98, 11.59, 13, 13.65, 14.74, 15.41, 15.12, 17.25, 18.46, 19.71, 21, 20.01//0.995
32 //12.74, 13.35, 14.88, 15.47, 14.94, 17.29, 18.6, 19.11, 20.9, 21.39, 21.36, 22.75, 23.66, 24.57, 26.04, 26.39//0.999
33 //};
34
35HoughZsFit::HoughZsFit(vector<HoughRecHit>* recHitCol)
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}
64
65HoughZsFit::~HoughZsFit(){
66}
67
68void HoughZsFit::doit(){
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}
378
379
380//void HoughZsFit::initPoint(){
381// for(int i=0;i<_hitSize;i++){
382// if(m_debug>0){
383// // cout<<"left s: "<<_recStereoHit[i].getsAmb(0)<<endl;
384// // cout<<"left z: "<<_recStereoHit[i].getzAmb(0)<<endl;
385// // cout<<"right s: "<<_recStereoHit[i].getsAmb(1)<<endl;
386// // cout<<"right z: "<<_recStereoHit[i].getzAmb(1)<<endl;
387// }
388// _vecPoint[i][0].first =_recStereoHit[i]->getsAmb(0);
389// _vecPoint[i][0].second=_recStereoHit[i]->getzAmb(0);
390// _vecPoint[i][1].first =_recStereoHit[i]->getsAmb(1);
391// _vecPoint[i][1].second=_recStereoHit[i]->getzAmb(1);
392// }
393//
394//}
395void HoughZsFit::leastFit(Line &linefit ,int N){
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}
416bool compare_zsfit(const HoughZsFit::Line& a,const HoughZsFit::Line& b){
417 return abs(a._chi)<abs(b._chi);
418}
419
420int HoughZsFit::leastLine(int nhit,double x[],double y[],double &k,double &b,double& chi2){
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}
447
448bool layer_in_track(const HoughRecHit* hita,const HoughRecHit* hitb){
449 return (*hita).getLayerId() > (*hitb).getLayerId();
450}
451void HoughZsFit::sortHit(){
452 std::sort (_recStereoHit.begin(),_recStereoHit.end(),layer_in_track);
453}
454void HoughZsFit::print(){
455 for(int i =0;i<_recStereoHit.size();i++){
456 cout<<"("<<_recStereoHit[i]->getLayerId()<<","<<_recStereoHit[i]->getWireId()<<")"<<endl;
457 }
458}
459
460int HoughZsFit::cgemfit(){
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}
587
588void HoughZsFit::hough(){
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
abs
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:101
s
XmlRpcServer s
Definition: HelloServer.cpp:11
layer_in_track
bool layer_in_track(const HoughRecHit *hita, const HoughRecHit *hitb)
Definition: HoughZsFit.cxx:448
compare_zsfit
bool compare_zsfit(const HoughZsFit::Line &a, const HoughZsFit::Line &b)
Definition: HoughZsFit.cxx:416
recHitCol
std::vector< HoughRecHit > recHitCol
Definition: InstallArea/include/MdcHoughFinder/MdcHoughFinder/Hough2D.h:24
layer_in_track
bool layer_in_track(const HoughRecHit *hita, const HoughRecHit *hitb)
Definition: HoughZsFit.cxx:448
compare_zsfit
bool compare_zsfit(const HoughZsFit::Line &a, const HoughZsFit::Line &b)
Definition: HoughZsFit.cxx:416
HoughZsFit::HoughZsFit
HoughZsFit(vector< HoughRecHit > *recHitCol)
Definition: HoughZsFit.cxx:35
HoughZsFit::leastLine
int leastLine(int, double x[], double y[], double &, double &, double &)
Definition: HoughZsFit.cxx:420
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::leastFit
void leastFit(Line &linefit, int N)
Definition: HoughZsFit.cxx:395
HoughZsFit::sortHit
void sortHit()
Definition: HoughZsFit.cxx:451