CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
HoughPeak.cxx
Go to the documentation of this file.
1#include "MdcHoughFinder/HoughPeak.h"
2double HoughPeak::m_dcut1;
3double HoughPeak::m_dcut2;
4double HoughPeak::m_factor=5.0;
5bool HoughPeak::m_debug=0;
6//double cut_deltaD[43] = {
7////0.33, 0.47, 0.45, //0.999
8////0.27, 0.39, 0.37, //0.995
9////0.23, 0.33, 0.33, //0.990
10//0.144294, 0.125804, 0.104357,
11//0,0,0,0,0,
12////0.297, 0.273, 0.273, 0.273, 0.324, 0.324, 0.34, 0.356, 0.435, 0.435, 0.475, 0.475, //0.999
13////0.261, 0.237, 0.225, 0.225, 0.236, 0.236, 0.244, 0.268, 0.315, 0.315, 0.355, 0.405, //0.995
14////0.237, 0.213, 0.195, 0.201, 0.204, 0.204, 0.204, 0.236, 0.245, 0.255, 0.295, 0.355, //0.990
15//0.0846565, 0.0700097, 0.056961, 0.0474704, 0.0394827, 0.0363385, 0.0343435, 0.0348542, 0.0385759, 0.0419172, 0.0473385, 0.0565087,
16//0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
17////4.05, 4.05, 4.15, 4.45, 4.35, 4.55, 4.55 //0.999
18////2.95, 2.95, 3.25, 3.45, 3.45, 3.75, 3.85 //0.995
19////2.35, 2.55, 2.75, 2.85, 3.05, 3.25, 3.35 /0.990
20//0.0998953, 0.0988695, 0.105915, 0.109863, 0.112257, 0.128139, 0.148522
21//};
22vector<double> HoughPeak::m_cut_deltaD(43,0);
23
24HoughPeak::HoughPeak(){
25}
26HoughPeak::~HoughPeak(){
27 // cout<<" HoughPeak destruct" <<endl;
28 for(int i =0;i<_houghPeakHitList.size();i++){
29 // delete _houghPeakHitList[i];
30 _houghPeakHitList[i]=NULL;
31 }
32 _houghPeakHitList.clear();
33}
34HoughPeak& HoughPeak::operator=(const HoughPeak& other)
35{
36 if( this == &other ) return *this;
37 _height=other._height;
38 _theta=other._theta;
39 _rho=other._rho;
40 _thetaBin=other._thetaBin;
41 _rhoBin=other._rhoBin;
42 _isCandiTrack=other._isCandiTrack;
43 _peakNum=other._peakNum;
44 _charge=other._charge;
45 if(_houghPeakHitList.size() != 0 ){
46 for(int i =0;i<_houghPeakHitList.size();i++){
47 _houghPeakHitList[i]=NULL;
48 }
49 _houghPeakHitList.clear();
50 }
51 _houghPeakHitList = other._houghPeakHitList;
52 //for(int i =0;i<other._houghPeakHitList.size();i++){
53 // const HoughHit* p_hit= new HoughHit( *(other._houghPeakHitList[i]));
54 // _houghPeakHitList.push_back(p_hit);
55 //}
56 return *this;
57}
58
59HoughPeak::HoughPeak(const HoughPeak& other)
60{
61 if(_houghPeakHitList.size() != 0 ){
62 for(int i =0;i<_houghPeakHitList.size();i++){
63 _houghPeakHitList[i]=NULL;
64 }
65 _houghPeakHitList.clear();
66 }
67 _height=other._height;
68 _theta=other._theta;
69 _rho=other._rho;
70 _thetaBin=other._thetaBin;
71 _rhoBin=other._rhoBin;
72 _houghPeakHitList = other._houghPeakHitList;
73 _isCandiTrack=other._isCandiTrack;
74 _peakNum=other._peakNum;
75 _charge=other._charge;
76 //for(int i =0;i<other._houghPeakHitList.size();i++){
77 // const HoughHit* p_hit= new HoughHit( *(other._houghPeakHitList[i]));
78 // _houghPeakHitList.push_back(p_hit);
79 //}
80}
81
82HoughPeak::HoughPeak(int itheta ,int irho, double theta,double rho,vector< const HoughHit* >** mapHitList,bool isCandiTrack,int peakNum){
83 _thetaBin=itheta;
84 _rhoBin=irho;
85 _theta=theta;
86 _rho=rho;
87 _houghPeakHitList = mapHitList[itheta][irho];
88 _isCandiTrack =isCandiTrack;
89 _peakNum=peakNum;
90 // int t_size=mapHitList[itheta][irho].size();
91 // for(int i =0;i<t_size;i++){
92 // const HoughHit* p_hit= new HoughHit( *(mapHitList[itheta][irho][i]));
93 // _houghPeakHitList.push_back(p_hit);
94 // }
95}
96
97HoughPeak::HoughPeak(int height,int itheta ,int irho, double theta,double rho,bool isCandiTrack,int peakNum,int charge){
98 _height=height;
99 _thetaBin=itheta;
100 _rhoBin=irho;
101 _theta=theta;
102 _rho=rho;
103 // _houghPeakHitList = mapHitList[itheta][irho];
104 _isCandiTrack =isCandiTrack;
105 _peakNum=peakNum;
106 _charge=charge;
107 // collectHits(hitlist);
108}
109
110int HoughPeak::getHitNum(int select ) const{
111 int size=_houghPeakHitList.size();
112 int n0,n1,n2,n3,n4,n5,n6;
113 n0=n1=n2=n3=n4=n5=n6=0;
114 for(int i=0;i<size;i++){
115 int cir= _houghPeakHitList[i]->getCirList();
116 int style= _houghPeakHitList[i]->getStyle();
117 // cout<<"peak compare ("<<_houghPeakHitList[i]->getLayerId()<<" ," <<_houghPeakHitList[i]->getWireId()<<" )"<<endl;
118 n0++;
119 if( style == 1 ) n1++;
120 if( style == 2 ) n2++;
121 if(_houghPeakHitList[i]->getDetectorType()==MDC){
122 int index = _houghPeakHitList[i]->digi()->getTrackIndex();
123 if( index < 0 ) n3++;
124 if( index >=0 && cir<=1 && style==0 ) n4++;
125 if( index <0 || cir>1 ) n5++;
126 if( index >=0 && cir==0 && style==0 ) n6++;
127 }
128 }
129 if( select == 0 ) return n0;
130 if( select == 1 ) return n1;
131 if( select == 2 ) return n2;
132 if( select == 3 ) return n3;
133 if( select == 4 ) return n4;
134 if( select == 5 ) return n5;
135 if( select == 6 ) return n6;
136}
137int HoughPeak::getHitNumA(int select ) const{
138 int size=_houghPeakHitList.size();
139 int n0,n1,n2,n3,n4,n5,n6;
140 n0=n1=n2=n3=n4=n5=n6=0;
141 for(int i=0;i<size;i++){
142 int cir= _houghPeakHitList[i]->getCirList();
143 int type = _houghPeakHitList[i]->getSlayerType();
144 int style= _houghPeakHitList[i]->getStyle();
145 if( type ==0 ) n0++;
146 if( type==0 && style == 1 ) n1++;
147 if( type==0 && style == 2 ) n2++;
148 if(_houghPeakHitList[i]->getDetectorType()==MDC){
149 int index = _houghPeakHitList[i]->digi()->getTrackIndex();
150 if( type==0 && index < 0 ) n3++;
151 if( type==0 && index >=0 && cir<=1 && style==0) n4++;
152 if( type==0 && (index <0 || cir>1) ) n5++;
153 if( type==0 && index >=0 && cir==0 && style==0) n6++;
154 }
155 }
156 if( select == 0 ) return n0;
157 if( select == 1 ) return n1;
158 if( select == 2 ) return n2;
159 if( select == 3 ) return n3;
160 if( select == 4 ) return n4;
161 if( select == 5 ) return n5;
162 if( select == 6 ) return n6;
163}
164int HoughPeak::getHitNumS(int select ) const{
165 int size=_houghPeakHitList.size();
166 int n0,n1,n2,n3,n4,n5,n6;
167 n0=n1=n2=n3=n4=n5=n6=0;
168 for(int i=0;i<size;i++){
169 int cir= _houghPeakHitList[i]->getCirList();
170 int type = _houghPeakHitList[i]->getSlayerType();
171 int style= _houghPeakHitList[i]->getStyle();
172 if( type !=0 ) n0++;
173 if( type!=0 && style == 1 ) n1++;
174 if( type!=0 && style == 2 ) n2++;
175 if(_houghPeakHitList[i]->getDetectorType()==MDC){
176 int index = _houghPeakHitList[i]->digi()->getTrackIndex();
177 if( type!=0 && index < 0 ) n3++;
178 if( type!=0 && index >=0 && cir<=1 && style==0) n4++;
179 if( type!=0 && (index <0 || cir>1) ) n5++;
180 if( type!=0 && index >=0 && cir==0 && style==0) n6++;
181 }
182 }
183 if( select == 0 ) return n0;
184 if( select == 1 ) return n1;
185 if( select == 2 ) return n2;
186 if( select == 3 ) return n3;
187 if( select == 4 ) return n4;
188 if( select == 5 ) return n5;
189 if( select == 6 ) return n6;
190}
191
192//int HoughPeak::collectHits(const HoughHitList& hitlist){
193int HoughPeak::collectHits( HoughHitList& hitlist){
194 if(m_debug>0)cout<<endl<<"in "<<__FILE__<<" collectHits:"<<endl;
195 //mean
196 //double p0[3] = {0.198127,0.124575,0.147637};
197 //double p1[3] = {-16.3613,-8.02201,-11.4986};
198 //double p2[3] = {0.235927,0.255211,0.269964};
199 //x3y4
200 //double p0[3] = {0.124424,0.0803608,0.514702};
201 //double p1[3] = {-2.41956,0.050727,-42.9705};
202 //double p2[3] = {0.35053,0.297292,0.355571};
203 double cirr_track = fabs(1./(_rho));
204 double cirx_track = (1./_rho)*cos(_theta);
205 double ciry_track = (1./_rho)*sin(_theta);
206 double pt = fabs( (1./_rho)/333.567 );
207 if(m_debug>0)cout<<"peakNum:"<<_peakNum<<" peak height:"<<_height<<" rho,theta:"<<_rho<<","<<_theta<<" circle center(x,y):("<<cirx_track<<","<<ciry_track<<")"<<endl;
208 //std::vector<HoughHit>::const_iterator iter= hitlist.getList().begin();
209 std::vector<HoughHit>::iterator iter= hitlist.getList().begin();
210 //int XVcluster(0);
211 for(; iter!= hitlist.getList().end(); iter++){
212 //const HoughHit *hit = &(*iter);
213 HoughHit *hit = &(*iter);
214 int layer = hit->getLayerId();
215 int wire = hit->getWireId();
216 if(hit->getDetectorType() == MDC && hit->getSlayerType()!=0 ) continue;
217 if(hit->getDetectorType() == CGEM && (hit->getCluster())->getflag() ==1) continue;
218 double cirx_hit = hit->getMidX();
219 double ciry_hit = hit->getMidY();
220 double cirr_hit = hit->getDriftDist();
221 if( _charge == 1 && (cirx_track*ciry_hit - ciry_track*cirx_hit>=0) ) continue; // suppose charge -1 FIXME
222 if( _charge ==-1 && (cirx_track*ciry_hit - ciry_track*cirx_hit<=0) ) continue; // suppose charge -1 FIXME
223 double deltaD ;
224 if(hit->getDetectorType() == MDC){
225 double l1l2 = sqrt( (cirx_hit-cirx_track)*(cirx_hit-cirx_track)+(ciry_hit-ciry_track)*(ciry_hit-ciry_track) );
226 if( l1l2>cirr_track ) deltaD = l1l2-cirr_track-cirr_hit;
227 if( l1l2<=cirr_track ) deltaD = l1l2-cirr_track+cirr_hit;
228 }else if(hit->getDetectorType() == CGEM){
229 double r_cgem = hit->getCgemGeomSvc()->getCgemLayer(layer)->getMiddleROfGapD()/10.;
230 double phi_cross = intersect_cylinder(_charge,cirx_track,ciry_track,r_cgem);
231 double phih_cluster = hit->getCluster()->getrecphi();
232 double dphi = phih_cluster - phi_cross;
233 if(dphi>M_PI)dphi-=2*M_PI;
234 if(dphi<-1*M_PI)dphi+=2*M_PI;
235 //double cut_factor = (1/(p0[layer]*pt*1000+p1[layer])+p2[layer])/cut_deltaD[layer];
236 //cout<<layer<<" "<<pt<<" "<<cut_factor<<endl;
237 //deltaD = r_cgem*dphi/cut_factor;
238 deltaD = r_cgem*dphi;
239 //hit->setDeltaD(deltaD);
240 }
241
242 //cal flight length
243
244 double theta_temp;
245 double l_temp;
246 //if(cirx_track==0||cirx_hit-cirx_track==0){
247 //theta_temp=0;
248 //}
249 //else{
250 theta_temp=M_PI-atan2(ciry_hit-ciry_track,cirx_hit-cirx_track)+atan2(ciry_track,cirx_track);
251 if(theta_temp>2*M_PI){
252 theta_temp=theta_temp-2*M_PI;
253 }
254 if(theta_temp<0){
255 theta_temp=theta_temp+2*M_PI;
256 }
257 //}
258 //cout<<" charge "<<_charge <<" "<<theta_temp<<endl;
259 if(_charge==-1) {
260 l_temp=cirr_track*theta_temp;
261 }
262 if(_charge==1) {
263 theta_temp=2*M_PI-theta_temp;
264 l_temp=cirr_track*(theta_temp);
265 }
266 if(m_debug>0)cout<<"hit: ("<<layer<<","<<wire<<") "<<" pt: "<<pt<<" deltaD: "<<deltaD<<" ,l: "<<l_temp<<endl;
267 //if( fabs(deltaD)<0.2 && l_temp<=35 ) _houghPeakHitList.push_back(hit); //suppose 60MeV FIXME
268 //if( fabs(deltaD)<0.2 && l_temp>35 && l_temp<=45) _houghPeakHitList.push_back(hit);
269 /*
270 double d_cut1 = m_dcut1;
271 if( pt<0.06 && fabs(deltaD)<d_cut1 && l_temp<=35 ) _houghPeakHitList.push_back(hit);
272 if( 0.06<pt&&pt<0.07 && fabs(deltaD)<d_cut1 && l_temp<=35 ) _houghPeakHitList.push_back(hit);
273 if( 0.07<pt&&pt<0.08 && fabs(deltaD)<d_cut1 && l_temp<=43 ) _houghPeakHitList.push_back(hit);
274 if( 0.08<pt&&pt<0.09 && fabs(deltaD)<d_cut1 && l_temp<=43 ) _houghPeakHitList.push_back(hit);
275 if( 0.09<pt&&pt<0.10 && fabs(deltaD)<d_cut1 && l_temp<=41 ) _houghPeakHitList.push_back(hit);
276 if( 0.10<pt&&pt<0.11 && fabs(deltaD)<d_cut1 && l_temp<=41 ) _houghPeakHitList.push_back(hit);
277 if( 0.11<pt&&pt<0.12 && fabs(deltaD)<d_cut1 && l_temp<=41 ) _houghPeakHitList.push_back(hit);
278 if( pt>0.12 && fabs(deltaD)<d_cut1 ) _houghPeakHitList.push_back(hit);
279
280 double d_cut2 = m_dcut2;
281 if( pt<0.06 && fabs(deltaD)<d_cut2 && l_temp>35&&l_temp<45 ) _houghPeakHitList.push_back(hit);
282 if( 0.06<pt&&pt<0.07 && fabs(deltaD)<d_cut2 && l_temp>35&& l_temp<=45) _houghPeakHitList.push_back(hit);
283 if( 0.07<pt&&pt<0.08 && fabs(deltaD)<d_cut2 && l_temp>43&& l_temp<=50) _houghPeakHitList.push_back(hit);
284 if( 0.10<pt&&pt<0.11 && fabs(deltaD)<d_cut2 && l_temp>80&& l_temp<=105) _houghPeakHitList.push_back(hit);
285 if( 0.11<pt&&pt<0.12 && fabs(deltaD)<d_cut2 && l_temp>80&& l_temp<=105) _houghPeakHitList.push_back(hit);
286 */
287 if( pt<0.06 && fabs(deltaD)<0.1&& l_temp<=35 ) _houghPeakHitList.push_back(hit);
288 if( 0.06<pt&&pt<0.07 && fabs(deltaD)<0.1 && l_temp<=35 ) _houghPeakHitList.push_back(hit);
289 if( 0.07<pt&&pt<0.08 && fabs(deltaD)<0.1 && l_temp<=43 ) _houghPeakHitList.push_back(hit);
290 if( 0.08<pt&&pt<0.09 && fabs(deltaD)<0.1 && l_temp<=43 ) _houghPeakHitList.push_back(hit);
291 if( 0.09<pt&&pt<0.10 && fabs(deltaD)<0.1 && l_temp<=41 ) _houghPeakHitList.push_back(hit);
292 if( 0.10<pt&&pt<0.11 && fabs(deltaD)<0.1 && l_temp<=41 ) _houghPeakHitList.push_back(hit);
293 if( 0.11<pt&&pt<0.12 && fabs(deltaD)<0.1 && l_temp<=41 ) _houghPeakHitList.push_back(hit);
294
295 if( pt<0.06 && fabs(deltaD)<0.1&& l_temp>35&&l_temp<45 ) _houghPeakHitList.push_back(hit);
296 if( 0.06<pt&&pt<0.07 && fabs(deltaD)<0.1 && l_temp>35&& l_temp<=45) _houghPeakHitList.push_back(hit);
297 if( 0.07<pt&&pt<0.08 && fabs(deltaD)<0.1 && l_temp>43&& l_temp<=50) _houghPeakHitList.push_back(hit);
298 if( 0.08<pt&&pt<0.09 && fabs(deltaD)<0.1 && l_temp>43&& l_temp<=50) _houghPeakHitList.push_back(hit);
299 // if( 0.09<pt&&pt<0.10 && fabs(deltaD)<0. && l_temp>41&& l_temp<=50) _houghPeakHitList.push_back(hit);
300 // if( 0.10<pt&&pt<0.11 && fabs(deltaD)<0. && l_temp>80&& l_temp<=105) _houghPeakHitList.push_back(hit);
301 // if( 0.11<pt&&pt<0.12 && fabs(deltaD)<0. && l_temp>80&& l_temp<=105) _houghPeakHitList.push_back(hit);
302 //if( pt>0.12 && fabs(deltaD)<0.1) _houghPeakHitList.push_back(hit);
303 //if( pt>0.12 && fabs(deltaD)<m_factor*cut_deltaD[layer]) _houghPeakHitList.push_back(hit);
304 if( pt>0.12 && fabs(deltaD)<m_factor*m_cut_deltaD[layer]) _houghPeakHitList.push_back(hit);
305 //_houghPeakHitList.push_back(hit);
306 }
307 return _houghPeakHitList.size();
308
309}
310
311//calculate phi value, only valid for the first half circle
312double HoughPeak::intersect_cylinder(double r_cylinder, double r_center, double phi_center, int charge)
313{
314 double phi;
315 if(charge==0)return 9999;
316 while(phi_center<0) phi_center += 2*M_PI;
317 while(phi_center>2*M_PI) phi_center -= 2*M_PI;
318 if(r_center<r_cylinder/2) return -9999;
319 double dphi = acos( r_cylinder/(2*r_center) );
320 if(charge<0) phi = phi_center + dphi;
321 else phi = phi_center - dphi;
322 while(phi<0) phi += 2*M_PI;
323 while(phi>2*M_PI) phi -= 2*M_PI;
324 return phi;
325}
326double HoughPeak::intersect_cylinder(int charge, double x_center, double y_center, double r_cylinder)
327{
328 double phi;
329 double phi_center = atan2(y_center,x_center);
330 double r_center = sqrt(x_center*x_center+y_center*y_center);
331 if(charge==0)return 9999;
332 while(phi_center<0) phi_center += 2*M_PI;
333 while(phi_center>2*M_PI) phi_center -= 2*M_PI;
334 if(r_center<r_cylinder/2) return -9999;
335 double dphi = acos( r_cylinder/(2*r_center) );
336 if(charge<0) phi = phi_center + dphi;
337 else phi = phi_center - dphi;
338 while(phi<0) phi += 2*M_PI;
339 while(phi>2*M_PI) phi -= 2*M_PI;
340 return phi;
341}
342
343void HoughPeak::print()
344{
345 //cout<<"........................................print peak "<<_peakNum<<"........................................"<<endl;
346 //cout<<"charge:"<<_charge<<endl;
347 //cout<<"isCandiTrack:"<<_isCandiTrack<<endl;
348 //cout<<"height:"<<_height<<endl;
349 //cout<<"thetaBin:"<<_thetaBin<<" theta:"<<_theta<<endl;
350 //cout<<"rhoBin:"<<_rhoBin<<" rho:"<<_rho<<endl;
351 cout<<"peak "<<_peakNum
352 <<" (theta,rho): "<<"("<<_theta<<", "<<_rho<<"),"
353 <<" height = "<<_height
354 <<", nHit = "<<getHitNumA(0)
355 <<", isCandiTrack:"<<_isCandiTrack
356 <<endl;
357 //printAllHit();
358}
359
360void HoughPeak::printAllHit() const{
361 int size=_houghPeakHitList.size();
362 if( size ==0 ) return ;
363 cout<<"Print hitlist in HoughPeak " <<endl;
364 cout<<"at center ("<<_theta<<" ,"<<_rho<<"): "<<size<<endl;
365 for(int i=0;i<size;i++){
366 int layer= _houghPeakHitList[i]->getLayerId();
367 int wire = _houghPeakHitList[i]->getWireId();
368 std::cout<<_houghPeakHitList[i]->getHitId()<<" ("<<layer<<","<<wire<<"), type "<<_houghPeakHitList[i]->getSlayerType()<<std::endl;
369 }
370}
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:101
sin
double sin(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
n2
int n2
Definition: SD0Tag.cxx:55
n1
int n1
Definition: SD0Tag.cxx:54
M_PI
#define M_PI
Definition: TConstant.h:4
HoughHitList::getList
const std::vector< HoughHit > & getList() const
Definition: HoughHitList.cxx:469
HoughPeak::getHitNumS
int getHitNumS(int) const
Definition: HoughPeak.cxx:164
HoughPeak::getHitNum
int getHitNum(int) const
Definition: HoughPeak.cxx:110
HoughPeak::intersect_cylinder
double intersect_cylinder(double r_cylinder, double r_center, double phi_center, int charge)
Definition: HoughPeak.cxx:312
HoughPeak::operator=
HoughPeak & operator=(const HoughPeak &other)
Definition: HoughPeak.cxx:34
HoughPeak::printAllHit
void printAllHit() const
Definition: HoughPeak.cxx:360
HoughPeak::getHitNumA
int getHitNumA(int) const
Definition: HoughPeak.cxx:137
HoughPeak::collectHits
int collectHits(HoughHitList &)
Definition: HoughPeak.cxx:193
HoughPeak::print
void print()
Definition: HoughPeak.cxx:343