BOSS 7.0.5
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;
10
11HoughZsFit::HoughZsFit(vector<HoughRecHit>* recHitCol)
12{
13 vector<HoughRecHit>::iterator iter=recHitCol->begin();
14 for(;iter!=recHitCol->end();iter++){
15 //if( fabs((*iter).getzAmb(0))>10 && fabs((*iter).getzAmb(1))>10 ) continue;
16 if( (*iter).getSlayerType() !=0 && (*iter).getflag()==0 && (*iter).getLayerId()<8 ) _recStereoHit.push_back(&(*iter));
17 }
18 _hitSize=_recStereoHit.size();
19 _pro_correct =0;
20 sortHit();
21 if(m_debug>0) print();
22 if(m_debug>0) cout<<"size of rechit in zs fit "<<_hitSize<<endl;
23
24// _vecPoint = new Point*[_hitSize];
25// for(int i=0;i<_hitSize;i++){
26// _vecPoint[i] = new Point[2];
27// }
28// initPoint();
29 if(_hitSize>=3) doit();
30 else {_tanl=0;_z0=0;}
31}
32
33HoughZsFit::~HoughZsFit(){
34}
35
36void HoughZsFit::doit(){
37
38// Point point0;
39// point0.first=0.;
40// point0.second=0.;
41
42 Line linefit[8];
43
44 for(int iline=0;iline<8;iline++){
45 linefit[iline]._pointCol.push_back( *(_recStereoHit[0]) );
46 linefit[iline]._pointCol.push_back( *(_recStereoHit[1]) );
47 linefit[iline]._pointCol.push_back( *(_recStereoHit[2]) );
48 }
49 if(m_debug>0) {
50 for(int i =0;i<3;i++){
51 cout<<" the first 3 hits ("<<i<<" "<< _recStereoHit[i]->getLayerId()<<","<<_recStereoHit[i]->getWireId()<<")"<<endl;
52 cout<<" left "<<(_recStereoHit[i])->getsAmb(0)<<" "<<(_recStereoHit[i])->getzAmb(0)<<endl;
53 cout<<" right "<<(_recStereoHit[i])->getsAmb(1)<<" "<<(_recStereoHit[i])->getzAmb(1)<<endl;
54 }
55 }
56
57 linefit[0]._pointCol[0].setAmb(0);
58 linefit[0]._pointCol[1].setAmb(0);
59 linefit[0]._pointCol[2].setAmb(0);
60
61 linefit[1]._pointCol[0].setAmb(0);
62 linefit[1]._pointCol[1].setAmb(0);
63 linefit[1]._pointCol[2].setAmb(1);
64
65 linefit[2]._pointCol[0].setAmb(0);
66 linefit[2]._pointCol[1].setAmb(1);
67 linefit[2]._pointCol[2].setAmb(0);
68
69 linefit[3]._pointCol[0].setAmb(0);
70 linefit[3]._pointCol[1].setAmb(1);
71 linefit[3]._pointCol[2].setAmb(1);
72
73
74 linefit[4]._pointCol[0].setAmb(1);
75 linefit[4]._pointCol[1].setAmb(0);
76 linefit[4]._pointCol[2].setAmb(0);
77
78 linefit[5]._pointCol[0].setAmb(1);
79 linefit[5]._pointCol[1].setAmb(0);
80 linefit[5]._pointCol[2].setAmb(1);
81
82 linefit[6]._pointCol[0].setAmb(1);
83 linefit[6]._pointCol[1].setAmb(1);
84 linefit[6]._pointCol[2].setAmb(0);
85
86 linefit[7]._pointCol[0].setAmb(1);
87 linefit[7]._pointCol[1].setAmb(1);
88 linefit[7]._pointCol[2].setAmb(1);
89
90 for(int i=0;i<8;i++){
91 linefit[i]._ambig.push_back( linefit[i]._pointCol[0].getAmbig() );
92 linefit[i]._ambig.push_back( linefit[i]._pointCol[1].getAmbig() );
93 linefit[i]._ambig.push_back( linefit[i]._pointCol[2].getAmbig() );
94 if(m_debug>0) cout<<" ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^begin line "<<i<<endl;
95 leastFit(linefit[i],3);
96 //judst if there is noise in the 3 outer
97 int line_size=3;
98 for(int j=3;j<_hitSize;j++){
99 double chi_last=linefit[i]._chi;
100 double k_last=linefit[i]._k;
101 double b_last=linefit[i]._b;
102 if(m_debug>0) cout<<"last "<<j<<" chi k b "<<chi_last<<" "<<k_last<<" "<<b_last<<endl;
103 linefit[i]._pointCol.push_back( *(_recStereoHit[j]) );
104 //int layer =linefit[i]._pointCol.back().getLayerId();
105 //int wire=linefit[i]._pointCol.back().getWireId();
106 //cout<<"???ambig ("<<layer<<","<<wire<<") "<<linefit[i]._pointCol.back().getLrTruth()<<endl;
107 linefit[i]._pointCol[line_size].setAmb(0);
108 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;
109 leastFit(linefit[i],line_size+1);
110 double chil=linefit[i]._chi;
111 double kl=linefit[i]._k;
112 double bl=linefit[i]._b;
113 if (linefit[i]._pointCol[line_size].getsPos()==-99) {
114 chil = 9999;
115 }
116 if(m_debug>0) cout<<"left "<<line_size<<" chi k b "<<chil<<" "<<kl<<" "<<bl<<endl;
117
118 linefit[i]._pointCol[line_size].setAmb(1);
119 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;
120 leastFit(linefit[i],line_size+1);
121 double chir=linefit[i]._chi;
122 double kr=linefit[i]._k;
123 double br=linefit[i]._b;
124 if (linefit[i]._pointCol[line_size].getsPos()==-99) {
125 chir = 9999;
126 }
127 if(m_debug>0) cout<<"right "<<line_size<<" chi k b "<<chir<<" "<<kr<<" "<<br<<endl;
128
129 if(chil<chir) {
130 linefit[i]._pointCol[line_size].setAmb(0);
131 linefit[i]._chi=chil;
132 linefit[i]._k=kl;
133 linefit[i]._b=bl;
134 linefit[i]._ambig.push_back(0);
135 //_recStereoHit[j]->setAmb(0);
136 }
137 else linefit[i]._ambig.push_back(1);
138 line_size++;
139
140 //test ambig with truth
141 int same_ambig = 1;
142 for(int ihit=0;ihit<line_size;ihit++){
143 int ambighit = linefit[i]._pointCol[ihit].getAmbig();
144 int ambigTruth = linefit[i]._pointCol[ihit].getLrTruth();
145 int layer =linefit[i]._pointCol[ihit].getLayerId();
146 int wire=linefit[i]._pointCol[ihit].getWireId();
147 //if (ambigTruth == -1) ambigTruth=1;
148 //else if (ambigTruth == 1) ambigTruth=0;
149 //else cout << "wrong in ambig Truth "<<endl;
150
151 //if(ambighit!= ambigTruth ) same_ambig=0;
152 }
153 //if(m_debug >0 ) cout << "same ambig ? "<<same_ambig<<endl;
154
155
156 double dChi = chi_last-linefit[i]._chi;
157 double dChi_n = (linefit[i]._chi)/(line_size-1)-(chi_last/line_size);
158 if(m_debug>0) cout<<"dChi: "<<dChi<<endl;
159 if(m_debug>0) cout<<"dChi/n: "<<dChi_n<<endl;
160 if(m_debug>0) cout<<endl;
161 // if( fabs(dChi) > 500.) //FIXME
162 if( fabs(dChi_n) > 25.) //FIXME
163 // linefit[i]._pointCol[j].setAmb(-999);
164 //_recStereoHit[j]->setAmb(-999);
165 //_recStereoHit[j]->setflag(-999);
166 {
167 linefit[i]._pointCol.pop_back();
168 linefit[i]._chi=chi_last;
169 linefit[i]._k=k_last;
170 linefit[i]._b=b_last;
171 linefit[i]._ambig.at(j)=-999;
172 line_size--;
173 }
174 }
175 //reject wrong line
176 if( (linefit[i]._pointCol[0].getsPos()==-99) || (linefit[i]._pointCol[1].getsPos()==-99) || (linefit[i]._pointCol[2].getsPos()==-99) ) {
177 linefit[i]._chi=99999;
178 }
179 }
180 std::sort(linefit,linefit+8,compare_zsfit);
181 for(int i=0;i<8;i++){
182 if(m_debug>0) cout<<"Line :"<<i<<" chis: "<<linefit[i]._chi<<" k,b: "<<linefit[i]._k<<" "<<linefit[i]._b<<endl;
183 int ambig_correct = 0;
184 for(int j=0;j<_hitSize;j++){
185 int ambig = linefit[i]._ambig.at(j);
186 int layer= _recStereoHit.at(j)->getLayerId();
187 int wire= _recStereoHit.at(j)->getWireId();
188 int style= _recStereoHit.at(j)->getStyle();
189 double l=-99;
190 if (ambig!=-999) l= _recStereoHit.at(j)->getsAmb(ambig);
191 double z=-99;
192 if (ambig!=-999) z= _recStereoHit.at(j)->getzAmb(ambig);
193 // cout<<"debug ("<<layer<<" ,"<<wire<<") style "<<style<<" ambig "<<ambig<<" s "<<l<<" z "<<z<<endl;
194 if (l==-99 && z==-99) ambig=-999;
195 if(m_debug>0) cout<<"("<<layer<<" ,"<<wire<<") style "<<style<<" ambig "<<ambig<<" s "<<l<<" z "<<z<<endl;
196 if(i==0) {
197 _recStereoHit[j]->setAmb(ambig); // when line 0 select
198 if( ambig ==-999) _recStereoHit[j]->setflag(-999); // when line 0 select
199 int ambigTruth = _recStereoHit.at(j)->getLrTruth();
200 if (ambigTruth == -1) ambigTruth=1;
201 else if (ambigTruth == 1) ambigTruth=0;
202 if(ambig ==ambigTruth) ambig_correct++;
203 //cout<<"ambig : "<<ambig<<" "<<ambigTruth<<endl;
204 _pro_correct = (double)ambig_correct/(double)_hitSize;
205 }
206 }
207 }
208 _tanl=linefit[0]._k;
209 _z0=linefit[0]._b;
210 if(m_debug>0) cout<<"z0 tanl : "<<_z0<<" "<<_tanl<<endl;
211
212 /*
213 linefit[0]._pointCol.push_back(_vecPoint[_hitSize-1][0]);
214 linefit[0]._pointCol.push_back(_vecPoint[_hitSize-2][0]);
215 linefit[0]._pointCol.push_back(_vecPoint[_hitSize-3][0]);
216
217 linefit[1]._pointCol.push_back(_vecPoint[_hitSize-1][0]);
218 linefit[1]._pointCol.push_back(_vecPoint[_hitSize-2][0]);
219 linefit[1]._pointCol.push_back(_vecPoint[_hitSize-3][1]);
220
221 linefit[2]._pointCol.push_back(_vecPoint[_hitSize-1][0]);
222 linefit[2]._pointCol.push_back(_vecPoint[_hitSize-2][1]);
223 linefit[2]._pointCol.push_back(_vecPoint[_hitSize-3][0]);
224
225 linefit[3]._pointCol.push_back(_vecPoint[_hitSize-1][0]);
226 linefit[3]._pointCol.push_back(_vecPoint[_hitSize-2][1]);
227 linefit[3]._pointCol.push_back(_vecPoint[_hitSize-3][1]);
228
229 linefit[4]._pointCol.push_back(_vecPoint[_hitSize-1][1]);
230 linefit[4]._pointCol.push_back(_vecPoint[_hitSize-2][0]);
231 linefit[4]._pointCol.push_back(_vecPoint[_hitSize-3][0]);
232
233 linefit[5]._pointCol.push_back(_vecPoint[_hitSize-1][1]);
234 linefit[5]._pointCol.push_back(_vecPoint[_hitSize-2][0]);
235 linefit[5]._pointCol.push_back(_vecPoint[_hitSize-3][1]);
236
237 linefit[6]._pointCol.push_back(_vecPoint[_hitSize-1][1]);
238 linefit[6]._pointCol.push_back(_vecPoint[_hitSize-2][1]);
239 linefit[6]._pointCol.push_back(_vecPoint[_hitSize-3][0]);
240
241 linefit[7]._pointCol.push_back(_vecPoint[_hitSize-1][1]);
242 linefit[7]._pointCol.push_back(_vecPoint[_hitSize-2][1]);
243 linefit[7]._pointCol.push_back(_vecPoint[_hitSize-3][1]);
244
245
246 for(int i=0;i<8;i++){
247 linefit[i]._amg=i*pow(10,(_hitSize-3));
248 if(m_debug>0) cout<<"Start line ooooooooooooooooooooooooooooooooooooooooooooooooo"<<i<<endl;
249 for(int j=4;j<_hitSize+1;j++){
250 if(m_debug>0) cout<<"Start add _vecPoint ****************************"<<j<<endl;
251 linefit[i]._pointCol.push_back(_vecPoint[_hitSize-j][0]);
252 leastFit(linefit[i],j);
253 double chil=linefit[i]._chi;
254 double kl=linefit[i]._k;
255 double bl=linefit[i]._b;
256
257 linefit[i]._pointCol[j-1]=_vecPoint[_hitSize-j][1];
258 leastFit(linefit[i],j);
259 double chir=linefit[i]._chi;
260 double kr=linefit[i]._k;
261 double br=linefit[i]._b;
262 if(m_debug>0){
263 cout<<chil<<" "<<kl<<" "<<bl<<endl;
264 cout<<chir<<" "<<kr<<" "<<br<<endl;
265 }
266
267 if(chil<chir) {
268 // cout<<"j: "<<j<<endl;
269 linefit[i]._pointCol[j-1]=_vecPoint[_hitSize-j][0];
270 linefit[i]._chi=chil;
271 linefit[i]._k=kl;
272 linefit[i]._b=bl;
273 linefit[i]._amg+=0*pow(10,(_hitSize-j));
274 }
275 else { linefit[i]._amg+=1*pow(10,(_hitSize-j)); }
276 }
277 if(m_debug>0) cout<<"Start add _vecPoint ****************************"<<"0"<<endl;
278 linefit[i]._pointCol.push_back(point0);
279 leastFit(linefit[i],_hitSize+1);
280 double chi0=linefit[i]._chi;
281 double k0=linefit[i]._k;
282 double b0=linefit[i]._b;
283 if(m_debug>0) cout<<chi0<<" "<<k0<<" "<<b0<<endl;
284}
285std::sort(linefit,linefit+8,compare_zsfit);
286for(int i=0;i<8;i++){
287 if(m_debug>0) cout<<"Line :"<<i<<" chis: "<<linefit[i]._chi<<" k,b: "<<linefit[i]._k<<" "<<linefit[i]._b<<" amg: "<<linefit[i]._amg<<endl;
288}
289_tanl=linefit[0]._k;
290_z0=linefit[0]._b;
291if(m_debug>0) cout<<"z0 tanl : "<<_z0<<" "<<_tanl<<endl;
292*/
293}
294
295
296//void HoughZsFit::initPoint(){
297// for(int i=0;i<_hitSize;i++){
298// if(m_debug>0){
299// // cout<<"left s: "<<_recStereoHit[i].getsAmb(0)<<endl;
300// // cout<<"left z: "<<_recStereoHit[i].getzAmb(0)<<endl;
301// // cout<<"right s: "<<_recStereoHit[i].getsAmb(1)<<endl;
302// // cout<<"right z: "<<_recStereoHit[i].getzAmb(1)<<endl;
303// }
304// _vecPoint[i][0].first =_recStereoHit[i]->getsAmb(0);
305// _vecPoint[i][0].second=_recStereoHit[i]->getzAmb(0);
306// _vecPoint[i][1].first =_recStereoHit[i]->getsAmb(1);
307// _vecPoint[i][1].second=_recStereoHit[i]->getzAmb(1);
308// }
309//
310//}
311void HoughZsFit::leastFit(Line &linefit ,int N){
312 double *x=new double[N+1];
313 double *y=new double[N+1];
314 for(int i=0;i<N;i++){
315 double a=linefit._pointCol[i].getsPos();
316 double b=linefit._pointCol[i].getzPos();
317 x[i]=a;
318 y[i]=b;
319 if(m_debug>0) cout<<" x "<<a<<" y "<<b<<endl;
320 }
321 x[N]=0.;
322 y[N]=0.;
323 double k(0),b(0),chi2(0);
324 leastLine(N+1,x,y,k,b,chi2);
325 linefit._k =k;
326 linefit._b =b;
327 linefit._chi=chi2;
328 delete []x;
329 delete []y;
330
331}
332bool compare_zsfit(const HoughZsFit::Line& a,const HoughZsFit::Line& b){
333 return abs(a._chi)<abs(b._chi);
334}
335
336int HoughZsFit::leastLine(int nhit,double x[],double y[],double &k,double &b,double& chi2){
337 double x_sum=0;
338 double y_sum=0;
339 double x2_sum=0;
340 double y2_sum=0;
341 double xy_sum=0;
342 for(int i=0;i<nhit;i++){
343 x_sum=x_sum+x[i];
344 y_sum=y_sum+y[i];
345 x2_sum=x2_sum+x[i]*x[i];
346 y2_sum=y2_sum+y[i]*y[i];
347 xy_sum=xy_sum+x[i]*y[i];
348 }
349 b=(x2_sum*y_sum-xy_sum*x_sum)/(nhit*x2_sum-x_sum*x_sum);
350 k=(nhit*xy_sum-x_sum*y_sum)/(nhit*x2_sum-x_sum*x_sum);
351
352 double yE[3];
353 for(int i=0;i<nhit;i++){
354 yE[i]=k*x[i]+b;
355 double chi2_temp;
356 if(yE[i]!=0) chi2_temp=(y[i]-yE[i])*(y[i]-yE[i])/yE[i]*yE[i];
357 else chi2_temp=0;
358 chi2+=chi2_temp;
359 }
360
361 return 1;
362}
363
364bool layer_in_track(const HoughRecHit* hita,const HoughRecHit* hitb){
365 return (*hita).getLayerId() > (*hitb).getLayerId();
366}
367void HoughZsFit::sortHit(){
368 std::sort (_recStereoHit.begin(),_recStereoHit.end(),layer_in_track);
369}
370void HoughZsFit::print(){
371 for(int i =0;i<_recStereoHit.size();i++){
372 cout<<"("<<_recStereoHit[i]->getLayerId()<<","<<_recStereoHit[i]->getWireId()<<")"<<endl;
373 }
374}
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:100
layer_in_track
bool layer_in_track(const HoughRecHit *hita, const HoughRecHit *hitb)
Definition: HoughZsFit.cxx:364
compare_zsfit
bool compare_zsfit(const HoughZsFit::Line &a, const HoughZsFit::Line &b)
Definition: HoughZsFit.cxx:332
recHitCol
std::vector< HoughRecHit > recHitCol
Definition: InstallArea/include/MdcHoughFinder/MdcHoughFinder/Hough2D.h:22
layer_in_track
bool layer_in_track(const HoughRecHit *hita, const HoughRecHit *hitb)
Definition: HoughZsFit.cxx:364
compare_zsfit
bool compare_zsfit(const HoughZsFit::Line &a, const HoughZsFit::Line &b)
Definition: HoughZsFit.cxx:332
HoughZsFit::HoughZsFit
HoughZsFit(vector< HoughRecHit > *recHitCol)
Definition: HoughZsFit.cxx:11
HoughZsFit::leastLine
int leastLine(int, double x[], double y[], double &, double &, double &)
Definition: HoughZsFit.cxx:336
HoughZsFit::doit
void doit()
Definition: HoughZsFit.cxx:36
HoughZsFit::print
void print()
Definition: HoughZsFit.cxx:370
HoughZsFit::leastFit
void leastFit(Line &linefit, int N)
Definition: HoughZsFit.cxx:311
HoughZsFit::sortHit
void sortHit()
Definition: HoughZsFit.cxx:367
y
double y[1000]
Definition: draw_charge_space_coarse_II.cxx:8
b
const double b
Definition: slope.cxx:9