BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
Hough2D.cxx
Go to the documentation of this file.
1#include "Hough2D.h"
2
3int Hough2D::m_debug = 0;
4TrkContextEv* Hough2D::_context = NULL;
5const IMdcCalibFunSvc* Hough2D::_mdcCalibFunSvc = NULL;
6int Hough2D::m_qualityFactor = 99999;
7double Hough2D::m_dropTrkDrCut = 1;
8// double Hough2D::m_dropTrkDrCut=999; //no limit of dr
9double Hough2D::m_dropTrkDzCut = 10;
10double Hough2D::m_dropTrkNhitCut = 5;
11double Hough2D::m_dropTrkChi2Cut = 99999.;
12double Hough2D::m_dropTrkChi2NdfCut = 99999.;
13
14Hough2D::Hough2D() {
15 //_m_gm = MdcDetector::instance(0);
16 //_m_gm = MdcDetector::instance(0);
17}
18
19// Hough2D::~Hough2D(){
20// for(int i =0;i<_MdcHitCol.size();i++){
21// cout<<" destruct " <<i<<endl;
22// delete _MdcHitCol[i];
23// }
24// }
25
26Hough2D::Hough2D( const Hough2D& other ) {
27 _circleR = other._circleR;
28 _circleX = other._circleX;
29 _circleY = other._circleY;
30 _d0 = other._d0;
31 _phi0 = other._phi0;
32 _omega = other._omega;
33 _charge = other._charge;
34 _pt = other._pt;
35 _nfit = other._nfit;
36 _chi2_aver2D = other._chi2_aver2D;
37 _bunchT0 = other._bunchT0;
38 //_m_gm = other._m_gm;
39 // _recHitVec=other._recHitVec;
40 _recHitVec = other._recHitVec;
41
42 // for(int i =0;i<other._MdcHitCol.size();i++){
43 // MdcHit* p_hit= new MdcHit( *(other._MdcHitCol[i]));
44 // _MdcHitCol.push_back(p_hit);
45 // }
46}
47
48Hough2D::Hough2D( recHitCol hitCol, double bunchtime )
49 : _recHitVec( hitCol ), _bunchT0( bunchtime ) {
50 //_m_gm = MdcDetector::instance(0);
51}
52
53// void Hough2D::print(){
54// cout<<"Print Hough2D "<<endl;
55// cout<<"d0: "<<_d0<<endl;
56// cout<<"phi0: "<<_phi0<<endl;
57// cout<<"omega: "<<_omega<<endl;
58// cout<<"charge: "<<_charge<<endl;
59//
60// int size=_recHitVec.size();
61// cout<<"Hit Candi size: "<<size<<endl;
62// for(int i=0;i<size;i++){
63// int layer= _recHitVec[i].getLayerId();
64// int wire = _recHitVec[i].getWireId();
65// std::cout<<"("<<layer<<","<<wire<<") "<<std::endl;
66// }
67// }
68
69int Hough2D::fitLeast() {
70 double x_sum = 0;
71 double y_sum = 0;
72 double x2_sum = 0;
73 double y2_sum = 0;
74 double x3_sum = 0;
75 double y3_sum = 0;
76 double x2y_sum = 0;
77 double xy2_sum = 0;
78 double xy_sum = 0;
79 double a1 = 0;
80 double a2 = 0;
81 double b1 = 0;
82 double b2 = 0;
83 double c1 = 0;
84 double c2 = 0;
85 // double x_aver,y_aver,r2;
86 int hitCandiSize = _recHitVec.size();
87 if ( hitCandiSize >= 3 )
88 {
89 for ( int i = 0; i < hitCandiSize; i++ )
90 {
91 // if( _recHitVec[i].getflag()!=0 ) continue;
92 // cout<<"size of hit in least2d "<<hitCandiSize<<endl;
93 // if(_recHitVec[i]->getSlayerType() !=0 ) continue;
94 x_sum = x_sum + ( _recHitVec[i].getMidX() );
95 y_sum = y_sum + ( _recHitVec[i].getMidY() );
96 x2_sum = x2_sum + ( _recHitVec[i].getMidX() ) * ( _recHitVec[i].getMidX() );
97 y2_sum = y2_sum + ( _recHitVec[i].getMidY() ) * ( _recHitVec[i].getMidY() );
98 x3_sum = x3_sum + ( _recHitVec[i].getMidX() ) * ( _recHitVec[i].getMidX() ) *
99 ( _recHitVec[i].getMidX() );
100 y3_sum = y3_sum + ( _recHitVec[i].getMidY() ) * ( _recHitVec[i].getMidY() ) *
101 ( _recHitVec[i].getMidY() );
102 x2y_sum = x2y_sum + ( _recHitVec[i].getMidX() ) * ( _recHitVec[i].getMidX() ) *
103 ( _recHitVec[i].getMidY() );
104 xy2_sum = xy2_sum + ( _recHitVec[i].getMidX() ) * ( _recHitVec[i].getMidY() ) *
105 ( _recHitVec[i].getMidY() );
106 xy_sum = xy_sum + ( _recHitVec[i].getMidX() ) * ( _recHitVec[i].getMidY() );
107 }
108 // hitCandiSize+=1000; //0 constrain
109 a1 = x2_sum - x_sum * x_sum / hitCandiSize;
110 a2 = xy_sum - x_sum * y_sum / hitCandiSize;
111 b1 = a2;
112 b2 = y2_sum - y_sum * y_sum / hitCandiSize;
113 c1 = ( x3_sum + xy2_sum - x_sum * ( x2_sum + y2_sum ) / hitCandiSize ) / 2.;
114 c2 = ( y3_sum + x2y_sum - y_sum * ( x2_sum + y2_sum ) / hitCandiSize ) / 2.;
115 // cout<<"a1 a2 b1 b2 c1 c2 "<<a1<<" "<<a2<<" "<<b1<<" "<<b2<<" "<<c1<<" "<<c2<<endl;
116
117 // x_aver=x_sum/hitCandiSize;
118 // y_aver=y_sum/hitCandiSize;
119
120 _circleX = ( b1 * c2 - b2 * c1 ) / ( b1 * b1 - a1 * b2 );
121 _circleY = ( b1 * c1 - a1 * c2 ) / ( b1 * b1 - a1 * b2 );
122 _circleR = sqrt( ( x2_sum + y2_sum - 2 * _circleX * x_sum - 2 * _circleY * y_sum ) /
123 hitCandiSize +
124 _circleX * _circleX + _circleY * _circleY );
125 // double r_temp=sqrt(r2);
126 // cout<<"circle: ("<<_circleX<<","<<_circleY<<","<<_circleR<<endl;
127
128 _d0 = sqrt( _circleX * _circleX + _circleY * _circleY ) - _circleR;
129 _phi0 = atan2( _circleY, _circleX ) + M_PI / 2.;
130 _omega = 1 / _circleR; // according with Q
131 double pt_temp = _circleR / 333.567;
132 if ( m_debug > 0 ) cout << " pt: " << pt_temp << endl;
133 if ( m_debug > 0 ) cout << "par: (" << _d0 << "," << _phi0 << "," << _omega << ")" << endl;
134 _pt = pt_temp;
135
136 _circleX_least = ( b1 * c2 - b2 * c1 ) / ( b1 * b1 - a1 * b2 );
137 _circleY_least = ( b1 * c1 - a1 * c2 ) / ( b1 * b1 - a1 * b2 );
138 _circleR_least = sqrt( ( x2_sum + y2_sum - 2 * _circleX * x_sum - 2 * _circleY * y_sum ) /
139 hitCandiSize +
140 _circleX * _circleX + _circleY * _circleY );
141 _pt_least = pt_temp;
142 return 1;
143 }
144 else { return 0; }
145}
146
147int Hough2D::fit() {
148 // replace in houghTrack
149 // if(_charge==-1){
150 // _d0=_d0;
151 // _omega=-1.*fabs(_omega);
152 // }
153 // if(_charge==1){
154 // _d0=-_d0;
155 // _omega=1.*fabs(_omega);
156 // _phi0=_phi0-M_PI;
157 // }
158 TrkExchangePar tt( _d0, _phi0, _omega, 0, 0 );
159 if ( m_debug > 0 )
160 cout << "q d0 phi0 omega: " << _charge << " " << _d0 << " " << _phi0 << " " << _omega
161 << endl;
162 TrkCircleMaker circlefactory;
163 float chisum = 0.;
164 newTrack = circlefactory.makeTrack( tt, chisum, *_context, _bunchT0 );
165 bool permitFlips = true;
166 bool lPickHits = true;
167 circlefactory.setFlipAndDrop( *newTrack, permitFlips, lPickHits );
168 // int nDigi = digiToHots(newTrack);
169 TrkHitList* trkHitList = newTrack->hits();
170 int digiId = 0;
171 vector<MdcHit*> t_hitCol;
172 std::vector<HoughRecHit>::iterator iter = _recHitVec.begin();
173 for ( ; iter != _recHitVec.end(); iter++, digiId++ )
174 {
175 if ( ( *iter ).getflag() != 0 ) continue;
176 // if( (*iter).getOuter()==1) continue;
177 // const MdcDetector* mgm = Global::m_gm;
178 const MdcDigi* aDigi = ( *iter ).digi();
179 // MdcHit *hit = new MdcHit(aDigi, _m_gm);
180 // MdcHit *hit = new MdcHit(aDigi, mgm);
181 MdcHit* hit = new MdcHit( aDigi, Global::m_gm );
182 // MdcHit mdcHit(aDigi,_m_gm);
183 // MdcHit *hit = &mdcHit;
184 t_hitCol.push_back( hit );
185 Identifier id = aDigi->identify();
186 int layer = MdcID::layer( id );
187 int wire = MdcID::wire( id );
188 // if(layer>=12) continue;
189 // if ((layer>=0&&layer<=7)||(layer>=20&&layer<=35)) {continue;}
190 hit->setCalibSvc( _mdcCalibFunSvc );
191 hit->setCountPropTime( true );
192 // new fltLen and ambig
193 int newAmbig = 0;
194 // calc. position of this point
195 MdcRecoHitOnTrack temp( *hit, newAmbig, _bunchT0 * 1.e9 ); // m_bunchT0 nano second here
196 MdcHitOnTrack* newhot = &temp;
197 newhot->setFltLen( ( *iter ).getRet().second ); // caculate by mc
198 double distoTrack = ( *iter ).getDisToTrack();
199
200 // double ddCut=1.0;
201 double ddCut = 1.0;
202 int use_in2d = 1;
203 if ( hit->driftTime( _bunchT0, 0 ) > 1000 ) use_in2d = 0; // >1000or800?
204 if ( hit->driftDist( _bunchT0, 0 ) > ddCut ) use_in2d = 0;
205 // if(m_debug>0) cout<<"flt : "<<(*iter).getRet().second <<endl;
206 if ( m_debug > 0 )
207 std::cout << " (" << layer << "," << wire << ",rt " << hit->rawTime() << ",dt "
208 << hit->driftTime( _bunchT0, 0 ) << ",dd " << hit->driftDist( _bunchT0, 0 )
209 << ") T0 " << _mdcCalibFunSvc->getT0( layer, wire ) << " timewalk "
210 << _mdcCalibFunSvc->getTimeWalk( layer, aDigi->getChargeChannel() )
211 << "dist: " << hit->driftDist( _bunchT0, 0 ) << " disToTrk " << distoTrack
212 << " use?: " << use_in2d << endl;
213 if ( use_in2d == 0 ) continue;
214 trkHitList->appendHot( newhot );
215 }
216 // newTrack->printAll(std::cout);
217
218 TrkHitList* qhits = newTrack->hits();
219 TrkErrCode err = qhits->fit();
220 const TrkFit* newFit = newTrack->fitResult();
221 int nActiveHit = newTrack->hots()->nActive();
222 int fit_stat = false;
223 double chi2 = -999.;
224 if ( !newFit || ( newFit->nActive() < 3 ) || err.failure() != 0 )
225 {
226 if ( m_debug > 0 )
227 {
228 cout << " global 2d fit failed ";
229 if ( newFit ) cout << " nAct " << newFit->nActive();
230 cout << "ERR1 failure =" << err.failure() << endl;
231 fit_stat = 0;
232 }
233 }
234 else
235 {
236 TrkExchangePar par = newFit->helix( 0. );
237 if ( abs( 1 / ( par.omega() ) ) > 0.03 )
238 {
239 fit_stat = 1;
240 chi2 = newFit->chisq();
241 if ( m_debug > 0 ) cout << "chi2 " << chi2 << endl;
242 }
243 else
244 {
245 fit_stat = 0;
246 chi2 = -999;
247 }
248
249 bool badQuality = false;
250 if ( fabs( chi2 ) > m_qualityFactor * m_dropTrkChi2Cut )
251 {
252 if ( m_debug > 0 )
253 {
254 std::cout << __FILE__ << " " << __LINE__ << " drop track by chi2 " << chi2 << " > "
255 << m_qualityFactor << " * " << m_dropTrkChi2Cut << std::endl;
256 }
257 badQuality = 1;
258 }
259 if ( fabs( par.d0() ) > m_qualityFactor * m_dropTrkDrCut )
260 {
261 if ( m_debug > 0 )
262 {
263 std::cout << __FILE__ << " " << __LINE__ << " drop track by d0 " << par.d0() << " > "
264 << m_qualityFactor << " * " << m_dropTrkDrCut << std::endl;
265 }
266 badQuality = 1;
267 }
268 if ( ( fabs( chi2 ) / qhits->nHit() ) > m_qualityFactor * m_dropTrkChi2NdfCut )
269 {
270 if ( m_debug > 0 )
271 {
272 std::cout << __FILE__ << " " << __LINE__ << " drop track by chi2/ndf"
273 << ( fabs( chi2 ) / qhits->nHit() ) << " > " << m_qualityFactor << " * "
274 << m_dropTrkChi2NdfCut << std::endl;
275 }
276 badQuality = 1;
277 }
278 if ( nActiveHit <= m_dropTrkNhitCut )
279 {
280 if ( m_debug > 0 )
281 {
282 std::cout << __FILE__ << " " << __LINE__ << " drop track by nhit" << nActiveHit
283 << " <= " << m_dropTrkNhitCut << std::endl;
284 }
285 badQuality = 1;
286 }
287 if ( badQuality ) fit_stat = 0;
288 }
289 if ( fit_stat == 1 )
290 {
291 if ( m_debug > 0 )
292 {
293 cout << " global 2d fit success" << endl;
294 cout << __FILE__ << __LINE__ << "AFTER fit 2d " << endl;
295 newTrack->print( std::cout );
296 }
297 TrkExchangePar par = newFit->helix( 0. );
298 double d0 = par.d0();
299 double phi0 = par.phi0();
300 double omega = par.omega();
301 double r_temp = -1. / par.omega();
302 double x_cirtemp = sin( par.phi0() ) * ( par.d0() + 1 / par.omega() ) *
303 -1.; // z axis verse,x_babar = -x_belle
304 double y_cirtemp = -1. * cos( par.phi0() ) * ( par.d0() + 1 / par.omega() ) * -1; //???
305 if ( m_debug > 0 )
306 {
307 cout << " circle after globle 2d: "
308 << "(" << x_cirtemp << "," << y_cirtemp << "," << r_temp << ")" << endl;
309 cout << "pt_rec: " << 1. / omega / 333.567 << endl;
310 }
311 _pt = 1. / omega / 333.567;
312 _circleX = x_cirtemp;
313 _circleY = y_cirtemp;
314 _circleR = _charge / omega;
315 _d0 = par.d0();
316 _phi0 = par.phi0();
317 _omega = par.omega();
318
319 int nfit2d = 0;
320 if ( m_debug > 1 ) cout << " hot list:" << endl;
321 TrkHotList::hot_iterator hotIter = qhits->hotList().begin();
322 int lay = ( (MdcHit*)( hotIter->hit() ) )->layernumber();
323 int wir = ( (MdcHit*)( hotIter->hit() ) )->wirenumber();
324 int hittemp = lay * 1000 + wir;
325 while ( hotIter != qhits->hotList().end() )
326 {
327 if ( m_debug > 1 )
328 {
329 cout << "(" << ( (MdcHit*)( hotIter->hit() ) )->layernumber() << ","
330 << ( (MdcHit*)( hotIter->hit() ) )->wirenumber() << ":" << hotIter->isActive()
331 << ") ";
332 }
333 if ( hotIter->isActive() == 1 ) nfit2d++;
334 hotIter++;
335 }
336 _nfit = nfit2d;
337 }
338 _chi2_aver2D = chi2 / _nfit;
339
340 for ( int i = 0; i < t_hitCol.size(); i++ ) { delete t_hitCol[i]; }
341 delete newTrack;
342
343 if ( m_debug > 0 ) cout << " in 2D fit number of Active hits : " << _nfit << endl;
344 return fit_stat;
345}
346
347/*
348int Hough2D::digiToHots(TrkRecoTrk* newTrack){
349 TrkHitList* trkHitList = newTrack->hits();
350 int digiId=0;
351 std::vector<HoughRecHit>::iterator iter = _recHitVec.begin();
352 for (;iter != _recHitVec.end(); iter++, digiId++) {
353 if( (*iter).getflag()!=0) continue;
354 //if( (*iter).getOuter()==1) continue;
355 const MdcDigi* aDigi = (*iter).digi();
356 //cout<<"Digi in hot "<<aDigi<<endl;
357 MdcHit *hit = new MdcHit(aDigi, _m_gm);
358 _MdcHitCol.push_back(hit);
359 Identifier id = aDigi->identify();
360 int layer = MdcID::layer(id);
361 int wire = MdcID::wire(id);
362 //if ((layer>=0&&layer<=7)||(layer>=20&&layer<=35)) {continue;}
363 hit->setCalibSvc(_mdcCalibFunSvc);
364 hit->setCountPropTime(true);
365 // new fltLen and ambig
366 int newAmbig = 0;
367 // calc. position of this point
368 MdcRecoHitOnTrack temp(*hit, newAmbig, _bunchT0*1.e9);//m_bunchT0 nano second here
369 MdcHitOnTrack* newhot = &temp;
370 newhot->setFltLen( (*iter).getRet().second); //caculate by mc
371
372 double ddCut;
373 if(layer<8) ddCut=1.0;
374 else ddCut=1.0;
375 int use_in2d=1;
376 if(hit->driftTime(_bunchT0,0)>800) use_in2d=0;;
377 if(hit->driftDist(_bunchT0,0)>ddCut) use_in2d=0;
378 //if(m_debug>0) cout<<"flt : "<<(*iter).getRet().second <<endl;
379 if(m_debug>0) std::cout<<" ("<<layer<<","<<wire<<",rt "<<hit->rawTime()<<",dt
380"<<hit->driftTime(_bunchT0,0)<<") T0 " << _mdcCalibFunSvc->getT0(layer,wire) <<" timewalk "<<
381_mdcCalibFunSvc->getTimeWalk(layer, aDigi->getChargeChannel())<< "dist:
382"<<hit->driftDist(_bunchT0,0)<<" use?: "<<use_in2d<<endl; if(use_in2d==0) continue;
383 trkHitList->appendHot(newhot);
384 }
385 return trkHitList->nHit();
386}
387*/
388void Hough2D::print() {
389 cout << " print hough2d " << endl;
390 cout << "par: " << _d0 << " " << _phi0 << " " << _omega << endl;
391}
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
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
TrkCircleMaker
TrkSimpleMaker< TrkCircleRep > TrkCircleMaker
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkFitter/TrkCircleMaker.h:14
M_PI
#define M_PI
Definition TConstant.h:4
Hough2D::m_dropTrkChi2NdfCut
static double m_dropTrkChi2NdfCut
Definition Hough2D.h:78
Hough2D::m_dropTrkDzCut
static double m_dropTrkDzCut
Definition Hough2D.h:76
Hough2D::_context
static TrkContextEv * _context
Definition Hough2D.h:69
Hough2D::fitLeast
int fitLeast()
Definition Hough2D.cxx:69
Hough2D::Hough2D
Hough2D()
Definition Hough2D.cxx:14
Hough2D::m_qualityFactor
static int m_qualityFactor
Definition Hough2D.h:73
Hough2D::m_debug
static int m_debug
Definition Hough2D.h:68
Hough2D::m_dropTrkChi2Cut
static double m_dropTrkChi2Cut
Definition Hough2D.h:74
Hough2D::print
void print()
Definition Hough2D.cxx:388
Hough2D::fit
virtual int fit()
Definition Hough2D.cxx:147
Hough2D::_mdcCalibFunSvc
static const IMdcCalibFunSvc * _mdcCalibFunSvc
Definition Hough2D.h:70
Hough2D::m_dropTrkNhitCut
static double m_dropTrkNhitCut
Definition Hough2D.h:77
Hough2D::m_dropTrkDrCut
static double m_dropTrkDrCut
Definition Hough2D.h:75
MdcHit::driftTime
double driftTime(double tof, double z) const
Definition MdcHit.cxx:142
MdcHit::driftDist
double driftDist(double, int, double, double, double) const
Definition MdcHit.cxx:157
MdcHit::setCalibSvc
void setCalibSvc(const IMdcCalibFunSvc *calibSvc)
Definition MdcHit.cxx:136
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0).
Definition MdcID.cxx:47
MdcID::wire
static int wire(const Identifier &id)
Definition MdcID.cxx:52
RawData::identify
virtual Identifier identify() const
Definition RawData.cxx:15
RawData::getChargeChannel
unsigned int getChargeChannel() const
Definition RawData.cxx:35
TrkAbsFit::chisq
virtual double chisq() const =0
TrkFit::nActive
virtual int nActive() const =0
TrkFit::helix
virtual TrkExchangePar helix(double fltL) const =0
TrkHitList::fit
TrkErrCode fit()
Definition TrkHitList.cxx:43
TrkHitList::appendHot
TrkHitOnTrk * appendHot(const TrkHitOnTrk *theHot)
Definition TrkHitList.cxx:73
TrkHitList::hotList
const TrkHotList & hotList() const
Definition TrkHitList.cxx:37
TrkHotList::hot_iterator
TrkHitOnTrkIter< TrkHotList::const_iterator_traits > hot_iterator
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkBase/TrkHotList.h:43
TrkSimpleMaker::setFlipAndDrop
bool setFlipAndDrop(TrkRecoTrk &, bool allowFlips, bool allowDrops) const
TrkSimpleMaker::makeTrack
TrkRecoTrk * makeTrack(const TrkExchangePar &helix, const double chi2, const TrkContext &, double trackT0) const
Global::m_gm
const MdcDetector * m_gm
Definition HoughGlobal.cxx:3