CGEM BOSS 6.6.5.f
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#include "McTruth/CgemMcHit.h"
13#include "TArrayI.h"
14#include "Identifier/CgemID.h"
15
16HoughHitList::HoughHitList(){
17}
18
19HoughHitList::HoughHitList(MdcDigiVec mdcDigiVec){
20 addMdcDigiList(mdcDigiVec);
21}
22
23HoughHitList::HoughHitList(RecCgemClusterCol* recCgemClusterCol){
24 addCgemClusterList(recCgemClusterCol);
25}
26
27HoughHitList::~HoughHitList(){
28 clearHoughHitList();
29}
30
31void HoughHitList::clearHoughHitList(){
32 std::vector<HoughHit> emptyHitList;
33 emptyHitList.swap(_houghHitList);
34}
35
36
37//int HoughHitList::nMdcHit() const{
38// int nMdcHit = 0;
39// std::vector<HoughHit>::const_iterator iter = _houghHitList.begin();
40// for (;iter != _houghHitList.end(); iter++) {
41// if(*iter->detectorType()==HoughHit::detectorType::MDC) nMdcHit++;
42// }
43// return nMdcHit;
44//}
45//
46//int HoughHitList::nCgemHit() const{
47// int nCgemHit = 0;
48// std::vector<HoughHit>::const_iterator iter = _houghHitList.begin();
49// for (;iter != _houghHitList.end(); iter++) {
50// if(*iter->detectorType()==HoughHit::detectorType::CGEM) nCgemHit++;
51// }
52// return nCgemHit;
53//}
54
55//FIXME
56void HoughHitList::addHit(HoughHit& a){
57 _houghHitList.push_back(a);
58}
59
60int HoughHitList::addMdcDigiList(MdcDigiVec mdcDigiVec){
61 int nHitAdd = 0;
62 int nHit = _houghHitList.size();
63 MdcDigiVec::iterator iter = mdcDigiVec.begin();
64 for (;iter != mdcDigiVec.end(); iter++) {
65 const MdcDigi* const digi = (*iter);
66 HoughHit hit(digi);//yzhang memory leakage?
67 if( hit.driftTime()>1000 ) continue; //to big hit ,maybe error hit
68 if( hit.driftTime()<=0) continue; //to big hit ,maybe error hit
69 // int stype= hit.getSlayerType();
70 hit.setHitId(nHit+nHitAdd);
71 _houghHitList.push_back(hit);
72 nHitAdd++;
73 }
74 //continousHit();
75 return nHitAdd;
76}
77
78int HoughHitList::addCgemClusterList(RecCgemClusterCol* cgemClusterVec){
79 int nHitAdd = 0;
80 int nHit = _houghHitList.size();
81 RecCgemClusterCol::iterator iter = cgemClusterVec->begin();
82 for (;iter != cgemClusterVec->end(); iter++) {
83 RecCgemCluster* cluster = (*iter);
84 HoughHit hit(cluster);
85 hit.setHitId(nHit+nHitAdd);
86 _houghHitList.push_back(hit);
87 nHitAdd++;
88 //std::cout<<__FILE__<<" "<<__LINE__<<" "<<cluster<<" id "<<cluster->getclusterid()<<" flag "<<cluster->getflag()<<" layer "<<(cluster)->getlayerid()<<std::endl;
89 }
90 //printAll();
91 //continousHit();
92 return nHitAdd;
93}
94
95int HoughHitList::addCgemClusterList(){
96 int nHitAdd = 0;
97 int nHit = _houghHitList.size();
98 //-- Get CGEM cluster
99 IDataProviderSvc* eventSvc = NULL;
100 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
101 SmartDataPtr<RecCgemClusterCol> cgemClusterVec(eventSvc, "/Event/Recon/RecCgemClusterCol");
102 RecCgemClusterCol::iterator iter = cgemClusterVec->begin();
103 for (;iter != cgemClusterVec->end(); iter++) {
104 RecCgemCluster* cluster = (*iter);
105 HoughHit hit(cluster);
106 hit.setHitId(nHit+nHitAdd);
107 _houghHitList.push_back(hit);
108 nHitAdd++;
109 }
110 //printAll();
111 //continousHit();
112 return nHitAdd;
113}
114
115bool small_layer(const HoughHit& a,const HoughHit& b){
116 return a.getWireId() < b.getWireId();
117}
118
119//judge continous -- remove noise before
120int HoughHitList::continousHit(){
121 const int maxCell[43] = {
122 40,44,48,56, 64,72,80,80,
123 76,76,88,88, 100,100,112,112, 128,128,140,140,
124 160,160,160,160, 176,176,176,176, 208,208,208,208, 240,240,240,240,
125 256,256,256,256, 288,288,288 };
126 vector<HoughHit> vec_HoughHit[43];
127 vector<HoughHit>::const_iterator iter=_houghHitList.begin();
128 for(int order=0; iter!=_houghHitList.end();iter++,order++){
129 int layer = (*iter).getLayerId();
130 //int wire= (*iter).getWireId();
131 vec_HoughHit[layer].push_back(*iter);
132 }
133 vector<HoughHit> vec_HoughHit_del;
134 //find continuous hits 5 hits or more
135 for(int i=0;i<43;i++){
136 //cout<<" layer "<<i<<" collect "<<vec_HoughHit[i].size()<<" hits "<<endl;
137 if(i<8) continue; //don't do for inner chamber
138 vector<int> vec_seeds;
139 std::sort(vec_HoughHit[i].begin(),vec_HoughHit[i].end(),small_layer);
140 if(vec_HoughHit[i].size()<=4) continue;
141 // for(unsigned int j=0;j<vec_HoughHit[i].size();j++){
142 // cout<<"("<<i<<","<<vec_HoughHit[i][j].getWireId()<<")";
143 // }
144 for(unsigned int j=0;j<vec_HoughHit[i].size()-4;j++){
145 vector<int>::iterator iter_hit = find(vec_seeds.begin(),vec_seeds.end(),vec_HoughHit[i][j].getWireId());
146 if( iter_hit!=vec_seeds.end() ) continue;
147 int wire_last=vec_HoughHit[i][j].getWireId();
148 int wire=-999;
149 // bool conti=true;
150 int seeds=1;
151 if(wire_last==0) {
152 // cout<<"maybe end conect to begin "<<endl;
153 for(unsigned int k=vec_HoughHit[i].size()-j-1;k>0;k--){
154 wire=vec_HoughHit[i][j+k].getWireId();
155 int charge = vec_HoughHit[i][j+k].getCharge();
156 int driftTime = vec_HoughHit[i][j+k].driftTime();
157 // cout<<"wirelast ("<<i<<","<<wire_last<<")"<<endl;
158 // cout<<"wire ("<<i<<","<<wire<<")"<<endl;
159 //if( (wire-maxCell[i]+1)!=wire_last || ( driftTime>0 && driftTime<800) ) break;
160 if( (wire-maxCell[i]+1)!=wire_last || (charge>0 && driftTime>0 && driftTime<800) ) break;
161 //cout<<"conti hit "<<"("<<i<<","<<wire<<")"<<endl;
162 wire_last= wire-maxCell[i];
163 seeds++;
164 // cout<<"seeds "<<seeds<<endl;
165 if(seeds==5) {
166 vec_seeds.push_back(wire+4-maxCell[i]);
167 vec_seeds.push_back(wire+3);
168 vec_seeds.push_back(wire+2);
169 vec_seeds.push_back(wire+1);
170 vec_seeds.push_back(wire);
171 }
172 if(seeds>5) vec_seeds.push_back(wire);
173 }
174 wire_last=0;
175
176 for(unsigned int k=1;k<vec_HoughHit[i].size()-j;k++){
177 vector<int>::iterator iter_hit0 = find(vec_seeds.begin(),vec_seeds.end(),vec_HoughHit[i][j].getWireId());
178 if( iter_hit0!=vec_seeds.end() ) continue;
179 wire=vec_HoughHit[i][j+k].getWireId();
180 int charge = vec_HoughHit[i][j+k].getCharge();
181 int driftTime = vec_HoughHit[i][j+k].driftTime();
182 // cout<<"wirelast ("<<i<<","<<wire_last<<")"<<endl;
183 // cout<<"wire ("<<i<<","<<wire<<")"<<endl;
184 if( wire<=maxCell[i] ) {
185 //if( (wire-1)!=wire_last || ( driftTime>0 && driftTime<800)) break;
186 if( (wire-1)!=wire_last || (charge>0 && driftTime>0 && driftTime<800)) break;
187 // cout<<"conti hit "<<"("<<i<<","<<wire<<")"<<endl;
188 wire_last= wire;
189 seeds++;
190 // cout<<"seeds "<<seeds<<endl;
191 if(seeds==5) {
192 vec_seeds.push_back(wire-4);
193 vec_seeds.push_back(wire-3);
194 vec_seeds.push_back(wire-2);
195 vec_seeds.push_back(wire-1);
196 vec_seeds.push_back(wire);
197 }
198 if(seeds>5) vec_seeds.push_back(wire);
199 }
200 }
201 }
202
203 else {
204 for(unsigned int k=1;k<vec_HoughHit[i].size()-j;k++){
205 wire=vec_HoughHit[i][j+k].getWireId();
206 int charge = vec_HoughHit[i][j+k].getCharge();
207 int driftTime = vec_HoughHit[i][j+k].driftTime();
208 // cout<<"wirelast ("<<i<<","<<wire_last<<")"<<endl;
209 // cout<<"wire ("<<i<<","<<wire<<") tdc "<<vec_HoughHit[i][j+k].driftTime()<<" adc "<<vec_HoughHit[i][j+k].getCharge()<<endl;
210 if( wire<=maxCell[i] ) {
211 //if( (wire-1)!=wire_last || ( driftTime>0 && driftTime<800)) break;
212 if( (wire-1)!=wire_last || (charge>0 && driftTime>0 && driftTime<800)) break;
213 // cout<<"conti hit "<<"("<<i<<","<<wire<<")"<<endl;
214 wire_last= wire;
215 seeds++;
216 // cout<<"seeds "<<seeds<<endl;
217 if(seeds==5) {
218 vec_seeds.push_back(wire-4);
219 vec_seeds.push_back(wire-3);
220 vec_seeds.push_back(wire-2);
221 vec_seeds.push_back(wire-1);
222 vec_seeds.push_back(wire);
223 }
224 if(seeds>5) vec_seeds.push_back(wire);
225 }
226 }
227 }
228
229 }
230 for(unsigned int ihit=0;ihit<vec_seeds.size();ihit++){
231 // cout<<"vec_seeds "<<"("<<i<<","<<vec_seeds[ihit]<<")"<<endl;
232 for(unsigned int jhit=0;jhit<vec_HoughHit[i].size();jhit++){
233 if(vec_HoughHit[i][jhit].getWireId()==vec_seeds[ihit]) vec_HoughHit_del.push_back(vec_HoughHit[i][jhit]);
234 }
235 }
236 }
237 for(unsigned int ihit=0;ihit<vec_HoughHit_del.size();ihit++){
238 // cout<<"("<<vec_HoughHit_del[ihit].getLayerId()<<","<<vec_HoughHit_del[ihit].getWireId()<<")"<<" tdc adc "<<vec_HoughHit_del[ihit].driftTime()<<","<<vec_HoughHit_del[ihit].getCharge()<<endl;
239 remove(vec_HoughHit_del[ihit]);
240 }
241}
242
243void HoughHitList::remove(const HoughHit& hit){
244 vector<HoughHit>::iterator iter = _houghHitList.begin();
245 for(;iter!=_houghHitList.end();iter++){
246 //if((hit.getDetectorType()==MDC && hit.digi()==(*iter).digi()) || (hit.getDetectorType()==CGEM && hit.getCluster()==(*iter).getCluster())) {
247 if(hit.getHitId()==(*iter).getHitId()){
248 _houghHitList.erase(iter);iter--;
249 // cout<<"remove ("<<hit.getLayerId()<<","<<hit.getWireId()<<")"<<endl;
250 }
251 }
252}
253
254int HoughHitList::addTruthInfo(std::map<int, const HepVector > mcTkPars){
255 //initialize digi info.
256 const MdcMcHit* truthHitPtr[43][288];
257 for(int i=0;i<43;i++){ for(int j=0;j<288;j++){ truthHitPtr[i][j]=NULL; } }
258 IDataProviderSvc* eventSvc = NULL;
259 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
260 SmartDataPtr<MdcMcHitCol> mdcMcHitCol(eventSvc,"/Event/MC/MdcMcHitCol");
261 if(mdcMcHitCol){
262 MdcMcHitCol::iterator iter_mchit = mdcMcHitCol->begin();
263 for (;iter_mchit != mdcMcHitCol->end(); iter_mchit++) {
264 const Identifier id= (*iter_mchit)->identify();
265 int l = MdcID::layer(id);
266 int w = MdcID::wire(id);
267 const MdcMcHit* mcHit = (*iter_mchit);
268 truthHitPtr[l][w] = mcHit;
269 //cout<<"MDC trackID: "<<(*iter_mchit)->getTrackIndex()<<" layer: "<<l<<" wire: "<<w<<endl;
270 }
271 }else{
272 std::cout<<__FILE__<<" Could not get MdcMcHitCol "<<std::endl;
273 return 0;
274 }
275 //bool cgem_true = true;
276 //if(cgem_true)
277 //{
278 eventSvc = NULL;
279 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
280 SmartDataPtr<Event::CgemMcHitCol> cgemMcHitCol(eventSvc, "/Event/MC/CgemMcHitCol");
281 if (!cgemMcHitCol){
282 std::cout<<__FILE__<<" Could not get CgemMcHitCol "<<std::endl;
283 return 0;
284 }
285 //}
286 //add truth info. to HoughHit in list
287 int nHitTruthAdd = 0;
288 for(std::vector<HoughHit>::iterator iter = _houghHitList.begin();iter!= _houghHitList.end(); iter++){
289 if((*iter).getDetectorType()==MDC){
290 int l = (*iter).getLayerId();
291 int w = (*iter).getWireId();
292 int mcTkId = (*iter).digi()->getTrackIndex();
293 if( mcTkId>=1000 ) mcTkId-=1000;
294 if(mcTkId>=0 && NULL!=truthHitPtr[l][w]) {
295 (*iter).setTruthInfo(truthHitPtr[l][w]);
296 //get truth track param to calc. flight length
297 if(mcTkPars.find(mcTkId)!= mcTkPars.end()){
298 (*iter).setStyle(0); // truth hit no decay no deltae
299 //HepVector helix = mcTkPars.find(mcTkId)->second;
300 //fltLen = (z_pos - dz)/tanl
301 //double rtemp = 333.567 * (1./fabs(helix[2]));
302 //double flt = (truthHitPtr[l][w]->getPositionZ()/10.-helix[3])/helix[4];
303 // (*iter).setFltLenTruth( flt );
304 /*
305 double cR = (-333.567)/helix[2];
306 double cX = (cR+helix[0])*cos(helix[1]);
307 double cY = (cR+helix[0])*sin(helix[1]);
308 double x= truthHitPtr[l][w]->getPositionX()/10;
309 double y= truthHitPtr[l][w]->getPositionY()/10;
310 //cout<<"circle ("<<cX<<","<<cY<<" )"<<cR<<endl;
311 //cout<<"hit ("<<x<<","<<y<<" )"<<endl;
312 double theta_temp;
313 if( cX==0 || x-cX==0 ){
314 theta_temp=0;
315 }
316 else{
317 theta_temp=M_PI-atan2(y-cY,x-cX)+atan2(cY,cX);
318 if(theta_temp>2*M_PI){
319 theta_temp=theta_temp-2*M_PI;
320 }
321 if(theta_temp<0){
322 theta_temp=theta_temp+2*M_PI;
323 }
324 if( fabs(helix[2])/helix[2] == 1)
325 theta_temp = 2*M_PI-theta_temp;
326 }
327 //cout<<theta_temp<<" cir "<<(*iter).getCirList()<<endl;
328 if( theta_temp > M_PI && (*iter).getCirList() ==0 ) {
329 (*iter).setCirList(1);
330 //cout<<"("<<l<<" ,"<<w<<" ) "<<" theta "<<theta_temp<<(*iter).getCirList()<<endl;
331 }
332 */
333 }else{
334 (*iter).setStyle(1); //decay
335 }
336 nHitTruthAdd++;
337 }
338 else {
339 (*iter).setStyle(2); //delta
340 }
341 if( mcTkId<0 ) (*iter).setStyle(-999); //noise
342 // if( mcTkId<0 ) (*iter).setOuter(0); // outer not influence noise
343 // if (m_debug>0) cout<<"style "<<(*iter).getStyle()<<endl;
344 //cout<<" hit ("<<l<<" ,"<<w<<" )"<<(*iter).getCirList()<<" style "<<(*iter).getStyle()<<" tq "<<(*iter).digi()->getTimeChannel()<<" "<<(*iter).digi()->getChargeChannel()<<endl;
345 //cout<<" hit ("<<l<<" ,"<<w<<" )"<<(*iter).getStyle()<<" trackindex "<<(*iter).digi()->getTrackIndex()<<endl;
346 //cout<<" hit in list "<<(*iter).getStyle()<<endl;
347 }else if((*iter).getDetectorType()==CGEM){
348 double dphi_min(M_PI);
349 CgemMcHit* cgemMcHit=NULL;
350 Event::CgemMcHitCol::iterator iter_truth = cgemMcHitCol->begin();
351 //cout<<cgemMcHitCol->size()<<endl;
352 bool findbegin(false),findend(false);
353 for(; iter_truth != cgemMcHitCol->end(); ++iter_truth){
354 double midx=((*iter_truth)->GetPositionXOfPrePoint()+(*iter_truth)->GetPositionXOfPostPoint())/2;
355 double midy=((*iter_truth)->GetPositionYOfPrePoint()+(*iter_truth)->GetPositionYOfPostPoint())/2;
356 double midz=((*iter_truth)->GetPositionZOfPrePoint()+(*iter_truth)->GetPositionZOfPostPoint())/2;
357 //cout<<"CGEM trackID: "<<(*iter_truth)->GetTrackID()
358 //<<" pid: "<<(*iter_truth)->GetPDGCode()
359 //<<" parentID: "<<(*iter_truth)->GetParentID()
360 //<<" layer: "<<(*iter_truth)->GetLayerID()
361 //<<" (x,y,z): ("<<midx<<", "<<midy<<", "<<midz<<")"<<endl;
362 if((*iter_truth)->GetLayerID()!=(*iter).getLayerId())continue;
363 double phi_true=atan2(midy,midx);
364 while(phi_true>2*M_PI)phi_true-=2*M_PI;
365 while(phi_true<0)phi_true+=2*M_PI;
366 double dphi=(*iter).getCluster()->getrecphi()-phi_true;
367 if(fabs(dphi)<dphi_min){
368 dphi_min=fabs(dphi);
369 cgemMcHit=(*iter_truth);
370 }
371 //cout<<"layer:"<<layer<<" phi_true:"<<phi_true/M_PI*180<<endl;
372 /*
373 TArrayI ID_Identifier = (*iter_truth)->GetIdentifier();
374 //cout<<(*iter_truth)->GetLayerID()<<" "<<ID_Identifier.GetSize()<<endl;;
375 for(int i=0;i<ID_Identifier.GetSize();i++){
376 int stripid = CgemID::strip(Identifier(ID_Identifier.At(i)));
377 int XV = CgemID::is_xstrip(Identifier(ID_Identifier.At(i)));
378 //cout<<stripid<<" XV:"<<XV<<endl;
379 if((*iter).getCluster()->getflag() == XV && (*iter).getCluster()->getclusterflagb() == stripid) findbegin = true;
380 if((*iter).getCluster()->getflag() == XV && (*iter).getCluster()->getclusterflage() == stripid) findend = true;
381 }
382 if(findbegin && findend){
383 //cgemMcHit=(*iter_truth);
384 (*iter).setTruthInfo(*iter_truth);
385 if(mcTkPars.find(cgemMcHit->GetTrackID())!= mcTkPars.end()){
386 (*iter).setStyle(0); // truth hit no decay no deltae FIXME
387 }else{
388 (*iter).setStyle(1); //decay FIXME
389 }
390 break;
391 }
392 //else{
393 //(*iter).setStyle(2); //delta
394 //}
395 //cout<<(*iter).getCluster()->getclusterflagb()<<"--"<<(*iter).getCluster()->getclusterflage()<<endl;
396 */
397 }
398 //if(!((findbegin && findend))(*iter).setStyle(2);
399 if(NULL!=cgemMcHit){
400 (*iter).setTruthInfo(cgemMcHit);
401 if(mcTkPars.find(cgemMcHit->GetTrackID())!= mcTkPars.end()){// FIXME should use TrackID of cluster(CgemDigi) instead of CgemMcHit
402 (*iter).setStyle(0); // truth hit no decay no deltae FIXME
403 }else{
404 (*iter).setStyle(1); //decay FIXME
405 }
406 }else{
407 (*iter).setStyle(2); //delta
408 }
409 }
410 }
411 //set cir
412 int cirlist=0;
413 if( _houghHitList.size()!=0) _houghHitList[0].setCirList(0);
414 for(unsigned int id = 1; id <_houghHitList.size(); id++){
415 // 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++;
416 if( ((_houghHitList[id].getLayerId() < _houghHitList[id-1].getLayerId()) )&& (cirlist==0||cirlist==2) ) cirlist++;
417 if(_houghHitList[id].getLayerId() > _houghHitList[id-1].getLayerId() && (cirlist==1||cirlist==3) ) cirlist++;
418 if(cirlist>3) cirlist=999;
419 _houghHitList[id].setCirList(cirlist);
420 }
421 int label_out_circle=999;//outermost layer
422 int label_num=0;//hit number of outermost layer
423 int label_fist=0;//first hit of outermost layer
424 int max_layer=0;
425 for(unsigned int id = 0; id <_houghHitList.size(); id++){
426 if(_houghHitList[id].getStyle()==-999) continue;
427 if( _houghHitList[id].getLayerId()>max_layer) max_layer=_houghHitList[id].getLayerId();
428 }
429 for(unsigned int id = 1; id <_houghHitList.size(); id++){
430 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();
431 }
432 for(int unsigned id = 0; id <_houghHitList.size(); id++){
433 if ( _houghHitList[id].getLayerId() == label_out_circle ) label_num++;
434 if( label_num==1 ) label_fist = id;
435 }
436 for(int id = 0; id <label_num/2; id++){
437 if( label_num %2 == 0 ) _houghHitList[label_fist+label_num/2+id].setCirList(1);
438 if( label_num %2 == 1 ) _houghHitList[label_fist+label_num/2+id+1].setCirList(1);
439 }
440
441 return nHitTruthAdd;
442}
443
444void HoughHitList::printAll(){
445 //print();
446 std::cout<<"MdcHoughFinder hit list: nHit="<<_houghHitList.size()<<std::endl;//<<" type "<<type()<<std::endl;
447 for(std::vector<HoughHit>::iterator iter = _houghHitList.begin();iter!= _houghHitList.end(); iter++){
448 (*iter).printAll();
449 }
450 std::cout<<std::endl;
451}
452void HoughHitList::printAll(int){
453 print();
454 for(std::vector<HoughHit>::iterator iter = _houghHitList.begin();iter!= _houghHitList.end(); iter++){
455 (*iter).printTruth();
456 }
457 std::cout<<std::endl;
458}
459
460void HoughHitList::print(){
461 std::cout<<"MdcHoughFinder hit list: nHit="<<_houghHitList.size()<<std::endl;//<<" type "<<type()<<std::endl;
462 for(std::vector<HoughHit>::iterator iter = _houghHitList.begin();iter!= _houghHitList.end(); iter++){
463 (*iter).print();
464 }
465 std::cout<<std::endl;
466}
467
468// get hit list
469const std::vector<HoughHit>& HoughHitList::getList() const{
470 return _houghHitList;
471}
472std::vector<HoughHit>& HoughHitList::getList() {
473 return _houghHitList;
474}
475int HoughHitList::getHitNum(int select ) const{
476 int size=_houghHitList.size();
477 int n0,n1,n2,n3,n4,n5,n6;
478 n0=n1=n2=n3=n4=n5=n6=0;
479 for(int i=0;i<size;i++){
480 int cir= _houghHitList[i].getCirList();
481 int style= _houghHitList[i].getStyle();
482 n0++;
483 if( style == 1 ) n1++;//decay
484 if( style == 2 ) n2++;//delta
485 if(_houghHitList[i].getDetectorType()==MDC){
486 int index = _houghHitList[i].digi()->getTrackIndex();
487 if( index < 0 ) n3++;//noise
488 if( index >=0 && cir<=1 && style==0) n4++;//HOT with truth
489 if( index <0 || cir>1 ) n5++;//noise or multi-circle
490 if( index >=0 && cir==0 && style==0) n6++;//HOT of first half circle with truth
491 }
492 }
493 if( select == 0 ) return n0;
494 else if( select == 1 ) return n1;
495 else if( select == 2 ) return n2;
496 else if( select == 3 ) return n3;
497 else if( select == 4 ) return n4;
498 else if( select == 5 ) return n5;
499 else if( select == 6 ) return n6;
500 else return -999;
501}
502int HoughHitList::getHitNumA(int select ) const{
503 int size=_houghHitList.size();
504 int n0,n1,n2,n3,n4,n5,n6;
505 n0=n1=n2=n3=n4=n5=n6=0;
506 for(int i=0;i<size;i++){
507 int cir= _houghHitList[i].getCirList();
508 int type = _houghHitList[i].getSlayerType();
509 int style= _houghHitList[i].getStyle();
510 if( type ==0 ) n0++;
511 if( type==0 && style == 1 ) n1++;
512 if( type==0 && style == 2 ) n2++;
513 if(_houghHitList[i].getDetectorType()==MDC){
514 int index = _houghHitList[i].digi()->getTrackIndex();
515 if( type==0 && index < 0 ) n3++;
516 if( type==0 && index >=0 && cir<=1 && style ==0) n4++;
517 if( type==0 && (index <0 || cir>1) ) n5++;
518 if( type==0 && index >=0 && cir==0 && style ==0) n6++;
519 }
520 }
521 if ( select == 0 ) return n0;
522 else if ( select == 1 ) return n1;
523 else if ( select == 2 ) return n2;
524 else if ( select == 3 ) return n3;
525 else if ( select == 4 ) return n4;
526 else if ( select == 5 ) return n5;
527 else if ( select == 6 ) return n6;
528 else return -999;
529}
530int HoughHitList::getHitNumS(int select ) const{
531 int size=_houghHitList.size();
532 int n0,n1,n2,n3,n4,n5,n6;
533 n0=n1=n2=n3=n4=n5=n6=0;
534 for(int i=0;i<size;i++){
535 int cir= _houghHitList[i].getCirList();
536 int type = _houghHitList[i].getSlayerType();
537 int style= _houghHitList[i].getStyle();
538 if( type !=0 ) n0++;
539 if( type!=0 && style == 1 ) n1++;
540 if( type!=0 && style == 2 ) n2++;
541 if(_houghHitList[i].getDetectorType()==MDC){
542 int index = _houghHitList[i].digi()->getTrackIndex();
543 if( type!=0 && index < 0 ) n3++;
544 if( type!=0 && index >=0 && cir<=1 && style==0) n4++;
545 if( type!=0 && (index <0 || cir>1) ) n5++;
546 if( type!=0 && index >=0 && cir==0 && style==0) n6++;
547 }
548 }
549 if( select == 0 ) return n0;
550 else if( select == 1 ) return n1;
551 else if( select == 2 ) return n2;
552 else if( select == 3 ) return n3;
553 else if( select == 4 ) return n4;
554 else if( select == 5 ) return n5;
555 else if( select == 6 ) return n6;
556 else return -999;
557
558}
RecCgemClusterCol
ObjectVector< RecCgemCluster > RecCgemClusterCol
Definition: Cgem/CgemRecEvent/CgemRecEvent-00-00-16/CgemRecEvent/RecCgemCluster.h:99
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:101
small_layer
bool small_layer(const HoughHit &a, const HoughHit &b)
Definition: HoughHitList.cxx:115
n2
int n2
Definition: SD0Tag.cxx:55
n1
int n1
Definition: SD0Tag.cxx:54
M_PI
#define M_PI
Definition: TConstant.h:4
HoughHitList::addMdcDigiList
int addMdcDigiList(MdcDigiVec mdcDigiVec)
Definition: HoughHitList.cxx:60
HoughHitList::remove
void remove(const HoughHit &hit)
Definition: HoughHitList.cxx:243
HoughHitList::getHitNumA
int getHitNumA(int) const
Definition: HoughHitList.cxx:502
HoughHitList::addHit
void addHit(HoughHit &a)
Definition: HoughHitList.cxx:56
HoughHitList::addTruthInfo
int addTruthInfo(std::map< int, const HepVector > mcTkPars)
Definition: HoughHitList.cxx:254
HoughHitList::getList
const std::vector< HoughHit > & getList() const
Definition: HoughHitList.cxx:469
HoughHitList::clearHoughHitList
void clearHoughHitList()
Definition: HoughHitList.cxx:31
HoughHitList::addCgemClusterList
int addCgemClusterList()
Definition: HoughHitList.cxx:95
HoughHitList::getHitNum
int getHitNum(int) const
Definition: HoughHitList.cxx:475
HoughHitList::getHitNumS
int getHitNumS(int) const
Definition: HoughHitList.cxx:530
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition: bak_Identifier-01-02-14/src/MdcID.cxx:49
MdcID::wire
static int wire(const Identifier &id)
Definition: bak_Identifier-01-02-14/src/MdcID.cxx:54