BOSS 7.0.4
BESIII Offline Software System
Loading...
Searching...
No Matches
HoughPeak.cxx
Go to the documentation of this file.
1#include "MdcHoughFinder/HoughPeak.h"
2#include<fstream>
3#include<sstream>
4double HoughPeak::m_dSigma_cut;
5double HoughPeak::m_dSigma_cut2;
6string HoughPeak::m_file;
7
8HoughPeak::HoughPeak(){
9}
10HoughPeak::~HoughPeak(){
11// cout<<" HoughPeak destruct" <<endl;
12 for(int i =0;i<_houghPeakHitList.size();i++){
13 // delete _houghPeakHitList[i];
14 _houghPeakHitList[i]=NULL;
15 }
16 _houghPeakHitList.clear();
17}
18HoughPeak& HoughPeak::operator=(const HoughPeak& other)
19{
20 if( this == &other ) return *this;
21 _height=other._height;
22 _theta=other._theta;
23 _rho=other._rho;
24 _thetaBin=other._thetaBin;
25 _rhoBin=other._rhoBin;
26 _isCandiTrack=other._isCandiTrack;
27 _peakNum=other._peakNum;
28 _charge=other._charge;
29 if(_houghPeakHitList.size() != 0 ){
30 for(int i =0;i<_houghPeakHitList.size();i++){
31 _houghPeakHitList[i]=NULL;
32 }
33 _houghPeakHitList.clear();
34 }
35 _houghPeakHitList = other._houghPeakHitList;
36 //for(int i =0;i<other._houghPeakHitList.size();i++){
37 // const HoughHit* p_hit= new HoughHit( *(other._houghPeakHitList[i]));
38 // _houghPeakHitList.push_back(p_hit);
39 //}
40 return *this;
41}
42
43HoughPeak::HoughPeak(const HoughPeak& other)
44{
45 if(_houghPeakHitList.size() != 0 ){
46 for(int i =0;i<_houghPeakHitList.size();i++){
47 _houghPeakHitList[i]=NULL;
48 }
49 _houghPeakHitList.clear();
50 }
51 _height=other._height;
52 _theta=other._theta;
53 _rho=other._rho;
54 _thetaBin=other._thetaBin;
55 _rhoBin=other._rhoBin;
56 _houghPeakHitList = other._houghPeakHitList;
57 _isCandiTrack=other._isCandiTrack;
58 _peakNum=other._peakNum;
59 _charge=other._charge;
60 //for(int i =0;i<other._houghPeakHitList.size();i++){
61 // const HoughHit* p_hit= new HoughHit( *(other._houghPeakHitList[i]));
62 // _houghPeakHitList.push_back(p_hit);
63 //}
64}
65
66HoughPeak::HoughPeak(int itheta ,int irho, double theta,double rho,vector< const HoughHit* >** mapHitList,bool isCandiTrack,int peakNum){
67 _thetaBin=itheta;
68 _rhoBin=irho;
69 _theta=theta;
70 _rho=rho;
71 _houghPeakHitList = mapHitList[itheta][irho];
72 _isCandiTrack =isCandiTrack;
73 _peakNum=peakNum;
74 // int t_size=mapHitList[itheta][irho].size();
75 // for(int i =0;i<t_size;i++){
76 // const HoughHit* p_hit= new HoughHit( *(mapHitList[itheta][irho][i]));
77 // _houghPeakHitList.push_back(p_hit);
78 // }
79}
80
81HoughPeak::HoughPeak(int height,int itheta ,int irho, double theta,double rho,bool isCandiTrack,int peakNum,int charge){
82 _height=height;
83 _thetaBin=itheta;
84 _rhoBin=irho;
85 _theta=theta;
86 _rho=rho;
87// _houghPeakHitList = mapHitList[itheta][irho];
88 _isCandiTrack =isCandiTrack;
89 _peakNum=peakNum;
90 _charge=charge;
91// collectHits(hitlist);
92}
93
94void HoughPeak::printAllHit() const{
95 int size=_houghPeakHitList.size();
96 if( size ==0 ) return ;
97 cout<<"Print hitlist in HoughPeak " <<endl;
98 cout<<"at center ("<<_theta<<" ,"<<_rho<<"): "<<size<<endl;
99 for(int i=0;i<size;i++){
100 int layer= _houghPeakHitList[i]->getLayerId();
101 int wire = _houghPeakHitList[i]->getWireId();
102 std::cout<<"Peak ("<<layer<<","<<wire<<") "<< _houghPeakHitList[i]<<std::endl;
103 }
104}
105
106int HoughPeak::getHitNum(int select ) const{
107 int size=_houghPeakHitList.size();
108 int n0,n1,n2,n3,n4,n5,n6;
109 n0=n1=n2=n3=n4=n5=n6=0;
110 for(int i=0;i<size;i++){
111 int cir= _houghPeakHitList[i]->getCirList();
112 int index = _houghPeakHitList[i]->digi()->getTrackIndex();
113 int style= _houghPeakHitList[i]->getStyle();
114 // cout<<"peak compare ("<<_houghPeakHitList[i]->getLayerId()<<" ," <<_houghPeakHitList[i]->getWireId()<<" )"<<endl;
115 n0++;
116 if( style == 1 ) n1++;
117 if( style == 2 ) n2++;
118 if( index < 0 ) n3++;
119 if( index >=0 && cir<=1 && style==0 ) n4++;
120 if( index <0 || cir>1 ) n5++;
121 if( index >=0 && cir==0 && style==0 ) n6++;
122 }
123 if( select == 0 ) return n0;
124 if( select == 1 ) return n1;
125 if( select == 2 ) return n2;
126 if( select == 3 ) return n3;
127 if( select == 4 ) return n4;
128 if( select == 5 ) return n5;
129 if( select == 6 ) return n6;
130}
131int HoughPeak::getHitNumA(int select ) const{
132 int size=_houghPeakHitList.size();
133 int n0,n1,n2,n3,n4,n5,n6;
134 n0=n1=n2=n3=n4=n5=n6=0;
135 for(int i=0;i<size;i++){
136 int cir= _houghPeakHitList[i]->getCirList();
137 int index = _houghPeakHitList[i]->digi()->getTrackIndex();
138 int type = _houghPeakHitList[i]->getSlayerType();
139 int style= _houghPeakHitList[i]->getStyle();
140 if( type ==0 ) n0++;
141 if( type==0 && style == 1 ) n1++;
142 if( type==0 && style == 2 ) n2++;
143 if( type==0 && index < 0 ) n3++;
144 if( type==0 && index >=0 && cir<=1 && style==0) n4++;
145 if( type==0 && (index <0 || cir>1) ) n5++;
146 if( type==0 && index >=0 && cir==0 && style==0) n6++;
147 }
148 if( select == 0 ) return n0;
149 if( select == 1 ) return n1;
150 if( select == 2 ) return n2;
151 if( select == 3 ) return n3;
152 if( select == 4 ) return n4;
153 if( select == 5 ) return n5;
154 if( select == 6 ) return n6;
155}
156int HoughPeak::getHitNumS(int select ) const{
157 int size=_houghPeakHitList.size();
158 int n0,n1,n2,n3,n4,n5,n6;
159 n0=n1=n2=n3=n4=n5=n6=0;
160 for(int i=0;i<size;i++){
161 int cir= _houghPeakHitList[i]->getCirList();
162 int index = _houghPeakHitList[i]->digi()->getTrackIndex();
163 int type = _houghPeakHitList[i]->getSlayerType();
164 int style= _houghPeakHitList[i]->getStyle();
165 if( type !=0 ) n0++;
166 if( type!=0 && style == 1 ) n1++;
167 if( type!=0 && style == 2 ) n2++;
168 if( type!=0 && index < 0 ) n3++;
169 if( type!=0 && index >=0 && cir<=1 && style==0) n4++;
170 if( type!=0 && (index <0 || cir>1) ) n5++;
171 if( type!=0 && index >=0 && cir==0 && style==0) n6++;
172 }
173 if( select == 0 ) return n0;
174 if( select == 1 ) return n1;
175 if( select == 2 ) return n2;
176 if( select == 3 ) return n3;
177 if( select == 4 ) return n4;
178 if( select == 5 ) return n5;
179 if( select == 6 ) return n6;
180}
181
182int HoughPeak::collectHits(const HoughHitList& hitlist){
183 //read in
184 double par[7][12];
185 ifstream input;
186 input.open(m_file.c_str());
187 string line;
188 int i=0;
189 while ( getline(input,line) ){
190 istringstream stream(line);
191 double sigma;
192 int j=0;
193 while(stream>>sigma){
194 par[i][j]=sigma;
195 j++;
196 }
197 i++;
198 }
199
200 double cirr_track = fabs(1./(_rho));
201 double cirx_track = (1./_rho)*cos(_theta);
202 double ciry_track = (1./_rho)*sin(_theta);
203 std::vector<HoughHit>::const_iterator iter= hitlist.getList().begin();
204 for(; iter!= hitlist.getList().end(); iter++){
205 const HoughHit *hit = &(*iter);
206 int layer = hit->getLayerId();
207 if( hit->getSlayerType()!=0 ) continue;
208 double cirx_hit = hit->getMidX();
209 double ciry_hit = hit->getMidY();
210 double cirr_hit = hit->getDriftDist();
211 if( _charge == 1 && (cirx_track*ciry_hit - ciry_track*cirx_hit>=0) ) continue; // suppose charge -1 FIXME
212 if( _charge ==-1 && (cirx_track*ciry_hit - ciry_track*cirx_hit<=0) ) continue; // suppose charge -1 FIXME
213 double l1l2 = sqrt( (cirx_hit-cirx_track)*(cirx_hit-cirx_track)+(ciry_hit-ciry_track)*(ciry_hit-ciry_track) );
214 double deltaD ;
215 if( l1l2>cirr_track ) deltaD = l1l2-cirr_track-cirr_hit;
216 if( l1l2<=cirr_track ) deltaD = l1l2-cirr_track+cirr_hit;
217
218 //cal flight length
219
220 double theta_temp;
221 double l_temp;
222 if(cirx_track==0||cirx_hit-cirx_track==0){
223 theta_temp=0;
224 }
225 else{
226 theta_temp=M_PI-atan2(ciry_hit-ciry_track,cirx_hit-cirx_track)+atan2(ciry_track,cirx_track);
227 if(theta_temp>2*M_PI){
228 theta_temp=theta_temp-2*M_PI;
229 }
230 if(theta_temp<0){
231 theta_temp=theta_temp+2*M_PI;
232 }
233 }
234 //cout<<" charge "<<_charge <<" "<<theta_temp<<endl;
235 if(_charge==-1) {
236 l_temp=cirr_track*theta_temp;
237 }
238 if(_charge==1) {
239 theta_temp=2*M_PI-theta_temp;
240 l_temp=cirr_track*(theta_temp);
241 }
242 double pt = fabs( (1./_rho)/333.567 );
243 //if( fabs(deltaD)<0.2 && l_temp<=35 ) _houghPeakHitList.push_back(hit); //suppose 60MeV FIXME
244 //if( fabs(deltaD)<0.2 && l_temp>35 && l_temp<=45) _houghPeakHitList.push_back(hit);
245 /*
246 double d_cut1 = m_dcut1;
247 if( pt<0.06 && fabs(deltaD)<d_cut1 && l_temp<=35 ) _houghPeakHitList.push_back(hit);
248 if( 0.06<pt&&pt<0.07 && fabs(deltaD)<d_cut1 && l_temp<=35 ) _houghPeakHitList.push_back(hit);
249 if( 0.07<pt&&pt<0.08 && fabs(deltaD)<d_cut1 && l_temp<=43 ) _houghPeakHitList.push_back(hit);
250 if( 0.08<pt&&pt<0.09 && fabs(deltaD)<d_cut1 && l_temp<=43 ) _houghPeakHitList.push_back(hit);
251 if( 0.09<pt&&pt<0.10 && fabs(deltaD)<d_cut1 && l_temp<=41 ) _houghPeakHitList.push_back(hit);
252 if( 0.10<pt&&pt<0.11 && fabs(deltaD)<d_cut1 && l_temp<=41 ) _houghPeakHitList.push_back(hit);
253 if( 0.11<pt&&pt<0.12 && fabs(deltaD)<d_cut1 && l_temp<=41 ) _houghPeakHitList.push_back(hit);
254 if( pt>0.12 && fabs(deltaD)<d_cut1 ) _houghPeakHitList.push_back(hit);
255
256 double d_cut2 = m_dcut2;
257 if( pt<0.06 && fabs(deltaD)<d_cut2 && l_temp>35&&l_temp<45 ) _houghPeakHitList.push_back(hit);
258 if( 0.06<pt&&pt<0.07 && fabs(deltaD)<d_cut2 && l_temp>35&& l_temp<=45) _houghPeakHitList.push_back(hit);
259 if( 0.07<pt&&pt<0.08 && fabs(deltaD)<d_cut2 && l_temp>43&& l_temp<=50) _houghPeakHitList.push_back(hit);
260 if( 0.10<pt&&pt<0.11 && fabs(deltaD)<d_cut2 && l_temp>80&& l_temp<=105) _houghPeakHitList.push_back(hit);
261 if( 0.11<pt&&pt<0.12 && fabs(deltaD)<d_cut2 && l_temp>80&& l_temp<=105) _houghPeakHitList.push_back(hit);
262 */
263
264
265 if( pt<0.06 && fabs(deltaD)<0.1&& l_temp<=35 ) _houghPeakHitList.push_back(hit);
266 if( 0.06<pt&&pt<0.07 && fabs(deltaD)<0.1 && l_temp<=35 ) _houghPeakHitList.push_back(hit);
267 if( 0.07<pt&&pt<0.08 && fabs(deltaD)<0.1 && l_temp<=43 ) _houghPeakHitList.push_back(hit);
268 if( 0.08<pt&&pt<0.09 && fabs(deltaD)<0.1 && l_temp<=43 ) _houghPeakHitList.push_back(hit);
269 if( 0.09<pt&&pt<0.10 && fabs(deltaD)<0.1 && l_temp<=41 ) _houghPeakHitList.push_back(hit);
270 if( 0.10<pt&&pt<0.11 && fabs(deltaD)<0.1 && l_temp<=41 ) _houghPeakHitList.push_back(hit);
271 if( 0.11<pt&&pt<0.12 && fabs(deltaD)<0.1 && l_temp<=41 ) _houghPeakHitList.push_back(hit);
272
273 if( pt<0.06 && fabs(deltaD)<0.1&& l_temp>35&&l_temp<45 ) _houghPeakHitList.push_back(hit);
274 if( 0.06<pt&&pt<0.07 && fabs(deltaD)<0.1 && l_temp>35&& l_temp<=45) _houghPeakHitList.push_back(hit);
275 if( 0.07<pt&&pt<0.08 && fabs(deltaD)<0.1 && l_temp>43&& l_temp<=50) _houghPeakHitList.push_back(hit);
276 if( 0.08<pt&&pt<0.09 && fabs(deltaD)<0.1 && l_temp>43&& l_temp<=50) _houghPeakHitList.push_back(hit);
277 // if( 0.09<pt&&pt<0.10 && fabs(deltaD)<0. && l_temp>41&& l_temp<=50) _houghPeakHitList.push_back(hit);
278 // if( 0.10<pt&&pt<0.11 && fabs(deltaD)<0. && l_temp>80&& l_temp<=105) _houghPeakHitList.push_back(hit);
279 // if( 0.11<pt&&pt<0.12 && fabs(deltaD)<0. && l_temp>80&& l_temp<=105) _houghPeakHitList.push_back(hit);
280 if( pt>0.12 && fabs(deltaD)<0.1) _houghPeakHitList.push_back(hit);
281
282/*
283 double n = m_dSigma_cut;
284 double n2 = m_dSigma_cut2;
285 if( layer<20){
286 if( 0.05<pt&&pt<0.06 ) {
287 if(layer<16 && fabs(deltaD)<n*par[0][layer-8]) _houghPeakHitList.push_back(hit);
288 if(layer>=16 && fabs(deltaD)<n2*par[0][layer-8]) _houghPeakHitList.push_back(hit);
289 }
290 if( 0.06<pt&&pt<0.07 ) {
291 if(layer<18 && fabs(deltaD)<n*par[1][layer-8]) _houghPeakHitList.push_back(hit);
292 if(layer>=18 && fabs(deltaD)<n2*par[0][layer-8]) _houghPeakHitList.push_back(hit);
293 }
294 if( 0.07<pt&&pt<0.08 ) {
295 if(fabs(deltaD)<n*par[2][layer-8]) _houghPeakHitList.push_back(hit);
296 }
297 if( 0.08<pt&&pt<0.09 ) {
298 if(fabs(deltaD)<n*par[3][layer-8]) _houghPeakHitList.push_back(hit);
299 }
300 if( 0.09<pt&&pt<0.10 ) {
301 if(fabs(deltaD)<n*par[4][layer-8]) _houghPeakHitList.push_back(hit);
302 }
303 if( 0.10<pt&&pt<0.11 ) {
304 if(fabs(deltaD)<n*par[5][layer-8]) _houghPeakHitList.push_back(hit);
305 }
306 if( 0.11<pt&&pt<0.12 ) {
307 if(fabs(deltaD)<n*par[6][layer-8]) _houghPeakHitList.push_back(hit);
308 }
309 if( 0.12<pt || pt<0.05 ) {
310 if(fabs(deltaD)<0.1) _houghPeakHitList.push_back(hit);
311 }
312 }
313 else if ( 0.12<pt || pt<0.05 ) {
314 if(fabs(deltaD)<0.1) _houghPeakHitList.push_back(hit);
315 }
316 */
317
318 }
319 return _houghPeakHitList.size();
320
321}
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:100
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:344
HoughPeak::getHitNumS
int getHitNumS(int) const
Definition: HoughPeak.cxx:156
HoughPeak::getHitNum
int getHitNum(int) const
Definition: HoughPeak.cxx:106
HoughPeak::operator=
HoughPeak & operator=(const HoughPeak &other)
Definition: HoughPeak.cxx:18
HoughPeak::printAllHit
void printAllHit() const
Definition: HoughPeak.cxx:94
HoughPeak::collectHits
int collectHits(const HoughHitList &)
Definition: HoughPeak.cxx:182
HoughPeak::getHitNumA
int getHitNumA(int) const
Definition: HoughPeak.cxx:131