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