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