BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
HoughHitList.cxx
Go to the documentation of this file.
1#include "HoughHitList.h"
2#include "GaudiKernel/Bootstrap.h"
3#include "GaudiKernel/IDataProviderSvc.h"
4#include "GaudiKernel/IInterface.h"
5#include "GaudiKernel/ISvcLocator.h"
6#include "GaudiKernel/Kernel.h"
7#include "GaudiKernel/Service.h"
8#include "GaudiKernel/SmartDataPtr.h"
9#include "Identifier/Identifier.h"
10#include "Identifier/MdcID.h"
11#include "McTruth/MdcMcHit.h"
12
14
16
18
20 std::vector<HoughHit> emptyHitList;
21 emptyHitList.swap( _houghHitList );
22}
23
24// int HoughHitList::nMdcHit() const{
25// int nMdcHit = 0;
26// std::vector<HoughHit>::const_iterator iter = _houghHitList.begin();
27// for (;iter != _houghHitList.end(); iter++) {
28// if(*iter->detectorType()==HoughHit::detectorType::MDC) nMdcHit++;
29// }
30// return nMdcHit;
31// }
32//
33// int HoughHitList::nCgemHit() const{
34// int nCgemHit = 0;
35// std::vector<HoughHit>::const_iterator iter = _houghHitList.begin();
36// for (;iter != _houghHitList.end(); iter++) {
37// if(*iter->detectorType()==HoughHit::detectorType::CGEM) nCgemHit++;
38// }
39// return nCgemHit;
40// }
41
42// FIXME
43void HoughHitList::addHit( HoughHit& a ) { _houghHitList.push_back( a ); }
44
46 int nHitAdd = 0;
47 MdcDigiVec::iterator iter = mdcDigiVec.begin();
48 for ( ; iter != mdcDigiVec.end(); iter++ )
49 {
50 const MdcDigi* const digi = ( *iter );
51 HoughHit hit( digi ); // yzhang memory leakage?
52 if ( hit.driftTime() > 1000 ) continue; // to big hit ,maybe error hit
53 if ( hit.driftTime() <= 0 ) continue; // to big hit ,maybe error hit
54 // int stype= hit.getSlayerType();
55 _houghHitList.push_back( hit );
56 nHitAdd++;
57 }
58 // continousHit();
59 return nHitAdd;
60}
61
62bool small_layer( const HoughHit& a, const HoughHit& b ) {
63 return a.getWireId() < b.getWireId();
64}
65
66// judge continous -- remove noise before
68 const int maxCell[43] = { 40, 44, 48, 56, 64, 72, 80, 80, 76, 76, 88,
69 88, 100, 100, 112, 112, 128, 128, 140, 140, 160, 160,
70 160, 160, 176, 176, 176, 176, 208, 208, 208, 208, 240,
71 240, 240, 240, 256, 256, 256, 256, 288, 288, 288 };
72 vector<HoughHit> vec_HoughHit[43];
73 vector<HoughHit>::const_iterator iter = _houghHitList.begin();
74 for ( int order = 0; iter != _houghHitList.end(); iter++, order++ )
75 {
76 int layer = ( *iter ).getLayerId();
77 // int wire= (*iter).getWireId();
78 vec_HoughHit[layer].push_back( *iter );
79 }
80 vector<HoughHit> vec_HoughHit_del;
81 // find continuous hits 5 hits or more
82 for ( int i = 0; i < 43; i++ )
83 {
84 // cout<<" layer "<<i<<" collect "<<vec_HoughHit[i].size()<<" hits "<<endl;
85 if ( i < 8 ) continue; // don't do for inner chamber
86 vector<int> vec_seeds;
87 std::sort( vec_HoughHit[i].begin(), vec_HoughHit[i].end(), small_layer );
88 if ( vec_HoughHit[i].size() <= 4 ) continue;
89 // for(unsigned int j=0;j<vec_HoughHit[i].size();j++){
90 // cout<<"("<<i<<","<<vec_HoughHit[i][j].getWireId()<<")";
91 // }
92 for ( unsigned int j = 0; j < vec_HoughHit[i].size() - 4; j++ )
93 {
94 vector<int>::iterator iter_hit =
95 find( vec_seeds.begin(), vec_seeds.end(), vec_HoughHit[i][j].getWireId() );
96 if ( iter_hit != vec_seeds.end() ) continue;
97 int wire_last = vec_HoughHit[i][j].getWireId();
98 int wire = -999;
99 // bool conti=true;
100 int seeds = 1;
101 if ( wire_last == 0 )
102 {
103 // cout<<"maybe end conect to begin "<<endl;
104 for ( unsigned int k = vec_HoughHit[i].size() - j - 1; k > 0; k-- )
105 {
106 wire = vec_HoughHit[i][j + k].getWireId();
107 int charge = vec_HoughHit[i][j + k].getCharge();
108 int driftTime = vec_HoughHit[i][j + k].driftTime();
109 // cout<<"wirelast ("<<i<<","<<wire_last<<")"<<endl;
110 // cout<<"wire ("<<i<<","<<wire<<")"<<endl;
111 // if( (wire-maxCell[i]+1)!=wire_last || ( driftTime>0 && driftTime<800) ) break;
112 if ( ( wire - maxCell[i] + 1 ) != wire_last ||
113 ( charge > 0 && driftTime > 0 && driftTime < 800 ) )
114 break;
115 // cout<<"conti hit "<<"("<<i<<","<<wire<<")"<<endl;
116 wire_last = wire - maxCell[i];
117 seeds++;
118 // cout<<"seeds "<<seeds<<endl;
119 if ( seeds == 5 )
120 {
121 vec_seeds.push_back( wire + 4 - maxCell[i] );
122 vec_seeds.push_back( wire + 3 );
123 vec_seeds.push_back( wire + 2 );
124 vec_seeds.push_back( wire + 1 );
125 vec_seeds.push_back( wire );
126 }
127 if ( seeds > 5 ) vec_seeds.push_back( wire );
128 }
129 wire_last = 0;
130
131 for ( unsigned int k = 1; k < vec_HoughHit[i].size() - j; k++ )
132 {
133 vector<int>::iterator iter_hit0 =
134 find( vec_seeds.begin(), vec_seeds.end(), vec_HoughHit[i][j].getWireId() );
135 if ( iter_hit0 != vec_seeds.end() ) continue;
136 wire = vec_HoughHit[i][j + k].getWireId();
137 int charge = vec_HoughHit[i][j + k].getCharge();
138 int driftTime = vec_HoughHit[i][j + k].driftTime();
139 // cout<<"wirelast ("<<i<<","<<wire_last<<")"<<endl;
140 // cout<<"wire ("<<i<<","<<wire<<")"<<endl;
141 if ( wire <= maxCell[i] )
142 {
143 // if( (wire-1)!=wire_last || ( driftTime>0 && driftTime<800)) break;
144 if ( ( wire - 1 ) != wire_last ||
145 ( charge > 0 && driftTime > 0 && driftTime < 800 ) )
146 break;
147 // cout<<"conti hit "<<"("<<i<<","<<wire<<")"<<endl;
148 wire_last = wire;
149 seeds++;
150 // cout<<"seeds "<<seeds<<endl;
151 if ( seeds == 5 )
152 {
153 vec_seeds.push_back( wire - 4 );
154 vec_seeds.push_back( wire - 3 );
155 vec_seeds.push_back( wire - 2 );
156 vec_seeds.push_back( wire - 1 );
157 vec_seeds.push_back( wire );
158 }
159 if ( seeds > 5 ) vec_seeds.push_back( wire );
160 }
161 }
162 }
163
164 else
165 {
166 for ( unsigned int k = 1; k < vec_HoughHit[i].size() - j; k++ )
167 {
168 wire = vec_HoughHit[i][j + k].getWireId();
169 int charge = vec_HoughHit[i][j + k].getCharge();
170 int driftTime = vec_HoughHit[i][j + k].driftTime();
171 // cout<<"wirelast ("<<i<<","<<wire_last<<")"<<endl;
172 // cout<<"wire ("<<i<<","<<wire<<") tdc
173 //"<<vec_HoughHit[i][j+k].driftTime()<<" adc
174 //"<<vec_HoughHit[i][j+k].getCharge()<<endl;
175 if ( wire <= maxCell[i] )
176 {
177 // if( (wire-1)!=wire_last || ( driftTime>0 && driftTime<800)) break;
178 if ( ( wire - 1 ) != wire_last ||
179 ( charge > 0 && driftTime > 0 && driftTime < 800 ) )
180 break;
181 // cout<<"conti hit "<<"("<<i<<","<<wire<<")"<<endl;
182 wire_last = wire;
183 seeds++;
184 // cout<<"seeds "<<seeds<<endl;
185 if ( seeds == 5 )
186 {
187 vec_seeds.push_back( wire - 4 );
188 vec_seeds.push_back( wire - 3 );
189 vec_seeds.push_back( wire - 2 );
190 vec_seeds.push_back( wire - 1 );
191 vec_seeds.push_back( wire );
192 }
193 if ( seeds > 5 ) vec_seeds.push_back( wire );
194 }
195 }
196 }
197 }
198 for ( unsigned int ihit = 0; ihit < vec_seeds.size(); ihit++ )
199 {
200 // cout<<"vec_seeds "<<"("<<i<<","<<vec_seeds[ihit]<<")"<<endl;
201 for ( unsigned int jhit = 0; jhit < vec_HoughHit[i].size(); jhit++ )
202 {
203 if ( vec_HoughHit[i][jhit].getWireId() == vec_seeds[ihit] )
204 vec_HoughHit_del.push_back( vec_HoughHit[i][jhit] );
205 }
206 }
207 }
208 for ( unsigned int ihit = 0; ihit < vec_HoughHit_del.size(); ihit++ )
209 {
210 // cout<<"("<<vec_HoughHit_del[ihit].getLayerId()<<","<<vec_HoughHit_del[ihit].getWireId()<<")"<<"
211 // tdc adc
212 //"<<vec_HoughHit_del[ihit].driftTime()<<","<<vec_HoughHit_del[ihit].getCharge()<<endl;
213 remove( vec_HoughHit_del[ihit] );
214 }
215}
216
217void HoughHitList::remove( const HoughHit& hit ) {
218 vector<HoughHit>::iterator iter = _houghHitList.begin();
219 for ( ; iter != _houghHitList.end(); iter++ )
220 {
221 if ( hit.digi() == ( *iter ).digi() )
222 {
223 _houghHitList.erase( iter );
224 iter--;
225 // cout<<"remove ("<<hit.getLayerId()<<","<<hit.getWireId()<<")"<<endl;
226 }
227 }
228}
229
230int HoughHitList::addTruthInfo( std::map<int, const HepVector> mcTkPars ) {
231 // initialize digi info.
232 const MdcMcHit* truthHitPtr[43][288];
233 for ( int i = 0; i < 43; i++ )
234 {
235 for ( int j = 0; j < 288; j++ ) { truthHitPtr[i][j] = NULL; }
236 }
237 IDataProviderSvc* eventSvc = NULL;
238 Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
239 SmartDataPtr<MdcMcHitCol> mdcMcHitCol( eventSvc, "/Event/MC/MdcMcHitCol" );
240 if ( mdcMcHitCol )
241 {
242 MdcMcHitCol::iterator iter_mchit = mdcMcHitCol->begin();
243 for ( ; iter_mchit != mdcMcHitCol->end(); iter_mchit++ )
244 {
245 const Identifier id = ( *iter_mchit )->identify();
246 int l = MdcID::layer( id );
247 int w = MdcID::wire( id );
248 const MdcMcHit* mcHit = ( *iter_mchit );
249 truthHitPtr[l][w] = mcHit;
250 }
251 }
252 else
253 {
254 std::cout << __FILE__ << " Could not get MdcMcHitCol " << std::endl;
255 return 0;
256 }
257
258 // add truth info. to HoughHit in list
259 int nHitTruthAdd = 0;
260 for ( std::vector<HoughHit>::iterator iter = _houghHitList.begin();
261 iter != _houghHitList.end(); iter++ )
262 {
263 int l = ( *iter ).getLayerId();
264 int w = ( *iter ).getWireId();
265 int mcTkId = ( *iter ).digi()->getTrackIndex();
266 if ( mcTkId >= 1000 ) mcTkId -= 1000;
267 if ( mcTkId >= 0 && NULL != truthHitPtr[l][w] )
268 {
269 ( *iter ).setTruthInfo( truthHitPtr[l][w] );
270 // get truth track param to calc. flight length
271 if ( mcTkPars.find( mcTkId ) != mcTkPars.end() )
272 {
273 ( *iter ).setStyle( 0 ); // truth hit no decay no deltae
274 // HepVector helix = mcTkPars.find(mcTkId)->second;
275 // fltLen = (z_pos - dz)/tanl
276 // double rtemp = 333.567 * (1./fabs(helix[2]));
277 // double flt = (truthHitPtr[l][w]->getPositionZ()/10.-helix[3])/helix[4];
278 // (*iter).setFltLenTruth( flt );
279 /*
280 double cR = (-333.567)/helix[2];
281 double cX = (cR+helix[0])*cos(helix[1]);
282 double cY = (cR+helix[0])*sin(helix[1]);
283 double x= truthHitPtr[l][w]->getPositionX()/10;
284 double y= truthHitPtr[l][w]->getPositionY()/10;
285 //cout<<"circle ("<<cX<<","<<cY<<" )"<<cR<<endl;
286 //cout<<"hit ("<<x<<","<<y<<" )"<<endl;
287 double theta_temp;
288 if( cX==0 || x-cX==0 ){
289 theta_temp=0;
290 }
291 else{
292 theta_temp=M_PI-atan2(y-cY,x-cX)+atan2(cY,cX);
293 if(theta_temp>2*M_PI){
294 theta_temp=theta_temp-2*M_PI;
295 }
296 if(theta_temp<0){
297 theta_temp=theta_temp+2*M_PI;
298 }
299 if( fabs(helix[2])/helix[2] == 1)
300 theta_temp = 2*M_PI-theta_temp;
301 }
302 //cout<<theta_temp<<" cir "<<(*iter).getCirList()<<endl;
303 if( theta_temp > M_PI && (*iter).getCirList() ==0 ) {
304 (*iter).setCirList(1);
305 //cout<<"("<<l<<" ,"<<w<<" ) "<<" theta "<<theta_temp<<(*iter).getCirList()<<endl;
306 }
307 */
308 }
309 else
310 {
311 ( *iter ).setStyle( 1 ); // decay
312 }
313 nHitTruthAdd++;
314 }
315 else
316 {
317 ( *iter ).setStyle( 2 ); // delta
318 }
319 if ( mcTkId < 0 ) ( *iter ).setStyle( -999 ); // noise
320 // if( mcTkId<0 ) (*iter).setOuter(0); // outer not influence noise
321 // if (m_debug>0) cout<<"style "<<(*iter).getStyle()<<endl;
322 // cout<<" hit ("<<l<<" ,"<<w<<" )"<<(*iter).getCirList()<<" style "<<(*iter).getStyle()<<"
323 // tq
324 // "<<(*iter).digi()->getTimeChannel()<<" "<<(*iter).digi()->getChargeChannel()<<endl;
325 // cout<<" hit ("<<l<<"
326 // ,"<<w<<" )"<<(*iter).getStyle()<<" trackindex "<<(*iter).digi()->getTrackIndex()<<endl;
327 // cout<<" hit in list
328 // "<<(*iter).getStyle()<<endl;
329 }
330
331 // set cir
332 int cirlist = 0;
333 if ( _houghHitList.size() != 0 ) _houghHitList[0].setCirList( 0 );
334 for ( unsigned int id = 1; id < _houghHitList.size(); id++ )
335 {
336 // if( ((_houghHitList[id].getLayerId() < _houghHitList[id-1].getLayerId()) ||
337 //(_houghHitList[id].getLayerId() == _houghHitList[id-1].getLayerId()&&
338 //_houghHitList[id-1].getLayerId() ==_houghHitList[id+1].getLayerId()) )&&
339 //(cirlist==0||cirlist==2) ) cirlist++;
340 if ( ( ( _houghHitList[id].getLayerId() < _houghHitList[id - 1].getLayerId() ) ) &&
341 ( cirlist == 0 || cirlist == 2 ) )
342 cirlist++;
343 if ( _houghHitList[id].getLayerId() > _houghHitList[id - 1].getLayerId() &&
344 ( cirlist == 1 || cirlist == 3 ) )
345 cirlist++;
346 if ( cirlist > 3 ) cirlist = 999;
347 _houghHitList[id].setCirList( cirlist );
348 }
349 int label_out_circle = 999;
350 int label_num = 0;
351 int label_fist = 0;
352 int max_layer = 0;
353 for ( unsigned int id = 0; id < _houghHitList.size(); id++ )
354 {
355 if ( _houghHitList[id].getStyle() == -999 ) continue;
356 if ( _houghHitList[id].getLayerId() > max_layer )
357 max_layer = _houghHitList[id].getLayerId();
358 }
359 for ( unsigned int id = 1; id < _houghHitList.size(); id++ )
360 {
361 if ( ( _houghHitList[id - 1].getLayerId() == _houghHitList[id].getLayerId() ) &&
362 ( _houghHitList[id].getLayerId() == _houghHitList[id + 1].getLayerId() ) &&
363 ( _houghHitList[id + 1].getLayerId() == _houghHitList[id + 2].getLayerId() ) &&
364 ( _houghHitList[id + 2].getLayerId() == _houghHitList[id + 3].getLayerId() ) &&
365 _houghHitList[id].getLayerId() == max_layer )
366 label_out_circle = _houghHitList[id - 1].getLayerId();
367 }
368 for ( int unsigned id = 0; id < _houghHitList.size(); id++ )
369 {
370 if ( _houghHitList[id].getLayerId() == label_out_circle ) label_num++;
371 if ( label_num == 1 ) label_fist = id;
372 }
373 for ( int id = 0; id < label_num / 2; id++ )
374 {
375 if ( label_num % 2 == 0 ) _houghHitList[label_fist + label_num / 2 + id].setCirList( 1 );
376 if ( label_num % 2 == 1 )
377 _houghHitList[label_fist + label_num / 2 + id + 1].setCirList( 1 );
378 }
379
380 return nHitTruthAdd;
381}
382
384 print();
385 for ( std::vector<HoughHit>::iterator iter = _houghHitList.begin();
386 iter != _houghHitList.end(); iter++ )
387 { ( *iter ).printAll(); }
388 std::cout << std::endl;
389}
391 print();
392 for ( std::vector<HoughHit>::iterator iter = _houghHitList.begin();
393 iter != _houghHitList.end(); iter++ )
394 { ( *iter ).printTruth(); }
395 std::cout << std::endl;
396}
397
399 std::cout << "MdcHoughFinder hit list: nHit=" << _houghHitList.size()
400 << std::endl; //<<" type "<<type()<<std::endl;
401}
402
403// get hit list
404const std::vector<HoughHit>& HoughHitList::getList() const { return _houghHitList; }
405std::vector<HoughHit>& HoughHitList::getList() { return _houghHitList; }
406int HoughHitList::getHitNum( int select ) const {
407 int size = _houghHitList.size();
408 int n0, n1, n2, n3, n4, n5, n6;
409 n0 = n1 = n2 = n3 = n4 = n5 = n6 = 0;
410 for ( int i = 0; i < size; i++ )
411 {
412 int cir = _houghHitList[i].getCirList();
413 int index = _houghHitList[i].digi()->getTrackIndex();
414 int style = _houghHitList[i].getStyle();
415 n0++;
416 if ( style == 1 ) n1++;
417 if ( style == 2 ) n2++;
418 if ( index < 0 ) n3++;
419 if ( index >= 0 && cir <= 1 && style == 0 ) n4++;
420 if ( index < 0 || cir > 1 ) n5++;
421 if ( index >= 0 && cir == 0 && style == 0 ) n6++;
422 }
423 if ( select == 0 ) return n0;
424 else if ( select == 1 ) return n1;
425 else if ( select == 2 ) return n2;
426 else if ( select == 3 ) return n3;
427 else if ( select == 4 ) return n4;
428 else if ( select == 5 ) return n5;
429 else if ( select == 6 ) return n6;
430 else return -999;
431}
432int HoughHitList::getHitNumA( int select ) const {
433 int size = _houghHitList.size();
434 int n0, n1, n2, n3, n4, n5, n6;
435 n0 = n1 = n2 = n3 = n4 = n5 = n6 = 0;
436 for ( int i = 0; i < size; i++ )
437 {
438 int cir = _houghHitList[i].getCirList();
439 int index = _houghHitList[i].digi()->getTrackIndex();
440 int type = _houghHitList[i].getSlayerType();
441 int style = _houghHitList[i].getStyle();
442 if ( type == 0 ) n0++;
443 if ( type == 0 && style == 1 ) n1++;
444 if ( type == 0 && style == 2 ) n2++;
445 if ( type == 0 && index < 0 ) n3++;
446 if ( type == 0 && index >= 0 && cir <= 1 && style == 0 ) n4++;
447 if ( type == 0 && ( index < 0 || cir > 1 ) ) n5++;
448 if ( type == 0 && index >= 0 && cir == 0 && style == 0 ) n6++;
449 }
450 if ( select == 0 ) return n0;
451 else if ( select == 1 ) return n1;
452 else if ( select == 2 ) return n2;
453 else if ( select == 3 ) return n3;
454 else if ( select == 4 ) return n4;
455 else if ( select == 5 ) return n5;
456 else if ( select == 6 ) return n6;
457 else return -999;
458}
459int HoughHitList::getHitNumS( int select ) const {
460 int size = _houghHitList.size();
461 int n0, n1, n2, n3, n4, n5, n6;
462 n0 = n1 = n2 = n3 = n4 = n5 = n6 = 0;
463 for ( int i = 0; i < size; i++ )
464 {
465 int cir = _houghHitList[i].getCirList();
466 int index = _houghHitList[i].digi()->getTrackIndex();
467 int type = _houghHitList[i].getSlayerType();
468 int style = _houghHitList[i].getStyle();
469 if ( type != 0 ) n0++;
470 if ( type != 0 && style == 1 ) n1++;
471 if ( type != 0 && style == 2 ) n2++;
472 if ( type != 0 && index < 0 ) n3++;
473 if ( type != 0 && index >= 0 && cir <= 1 && style == 0 ) n4++;
474 if ( type != 0 && ( index < 0 || cir > 1 ) ) n5++;
475 if ( type != 0 && index >= 0 && cir == 0 && style == 0 ) n6++;
476 }
477 if ( select == 0 ) return n0;
478 else if ( select == 1 ) return n1;
479 else if ( select == 2 ) return n2;
480 else if ( select == 3 ) return n3;
481 else if ( select == 4 ) return n4;
482 else if ( select == 5 ) return n5;
483 else if ( select == 6 ) return n6;
484 else return -999;
485}
double w
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
bool small_layer(const HoughHit &a, const HoughHit &b)
std::vector< MdcDigi * > MdcDigiVec
int n2
Definition SD0Tag.cxx:59
int n1
Definition SD0Tag.cxx:58
int addMdcDigiList(MdcDigiVec mdcDigiVec)
void remove(const HoughHit &hit)
int getHitNumA(int) const
void addHit(HoughHit &a)
HoughHitType type() const
void continousHit()
int addTruthInfo(std::map< int, const HepVector > mcTkPars)
const std::vector< HoughHit > & getList() const
void clearHoughHitList()
int getHitNum(int) const
int getHitNumS(int) const
const MdcDigi * digi() const
Definition HoughHit.h:53
double driftTime() const
Definition HoughHit.cxx:170
int getWireId() const
Definition HoughHit.h:61
static int layer(const Identifier &id)
Values of different levels (failure returns 0).
Definition MdcID.cxx:47
static int wire(const Identifier &id)
Definition MdcID.cxx:52