CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
McTruth.cxx
Go to the documentation of this file.
1#include "PrintMcInfo/PrintMcInfo.h"
2using namespace std;
3
4void PrintMcInfo::mkmap()
5{
6 string str,str1,str2,str3,str4;
7 int id;
8
9 char* pdt_path=getenv("BESEVTGENROOT");
10 if(0 == pdt_path)
11 std::cerr<<"ERROR : Can't get env $BESEVTGENROOT"<<std::endl;
12 std::string pdtname=std::string(pdt_path) + "/share/pdt.table";
13 ifstream fin(pdtname.c_str());
14
15 while(getline(fin,str))
16 {
17 fin>>str;
18 if(str == "add")
19 {
20 fin>>str1>>str2>>str3>>str4;
21 id = atoi(str4.c_str());
22 map_pid.insert(map<int, string>::value_type(id,str3));
23 }
24 }
25}
26
27void PrintMcInfo::printTitle(ofstream& os,int m_OutputLevel)
28{
29 if(m_OutputLevel==1)
30 {
31 os.setf(ios::left);
32 os<<" "<<setw(27)<<"cosTeta phi"<<setw(50)<<"P4 : momentum Energy"
33 <<setw(10)<<" "<<setw(40)<<"initialPosition"
34 <<setw(15)<<"finalPosition"<<endl;
35 os.setf(ios::left);
36 for (int m = 0; m < 29; m++) {os<< " ";}
37 os<<"Px, Py, Pz , E (Gev)";
38 for (int m = 0; m < 35; m++) {os<< " ";}
39 os<<"x(cm) y(cm) z(cm)";
40 for (int m = 0; m < 25; m++) {os<< " ";}
41 os <<"x(cm) y(cm) z(cm)"<<endl;
42 }
43 /* if(m_OutputLevel==2)
44 {
45 os<<'\t'<<"MdcHit"<<'\t'<<'\t'<<"TofHit"<<'\t'<<'\t'<<'\t'<<"EmcHit"<<'\t'<<'\t'<<'\t'<<"MucHit"
46 <<endl;
47 os<<'\t'<<"(layer wire)"<<" "<<"(b/ec layer phi_module)"<<" "<<"(b/ec theta phi)"<<" "<<"(b/ec segment layer strip)"<<endl;
48 }
49 */
50}
51
52
53
54
55void PrintMcInfo::printTree(ofstream& os,Event::McParticle* pmcp,int m_OutputLevel,int tab=0)
56{
57 bool decay = (pmcp->daughterList()).size()==0?false:true;
58 if(decay)
59 {
60 for (int m = 0; m < tab; m++) {os<< '\t';}
61
62 m_pid = (*pmcp).particleProperty();
63 m_trkIndex = (*pmcp).trackIndex();
64 map<int,string>::iterator iterPar;
65 iterPar = map_pid.find(m_pid);
66 if (iterPar!=map_pid.end())
67 os<<iterPar->second<<"["<<m_trkIndex<<"]"<<" "<<"->";
68
69 //loop: to get the daughters of the decayed particle
70 SmartRefVector<Event::McParticle>::const_iterator begin = (pmcp->daughterList()).begin();
71 SmartRefVector<Event::McParticle>::const_iterator end = (pmcp->daughterList()).end();
72 SmartRefVector<Event::McParticle>::const_iterator it;
73 for(it = begin;it != end;it++)
74 {
75 m_trkIndex = (*it)->trackIndex();
76 map<int,string>::iterator iter;
77 iter = map_pid.find((*it)->particleProperty());
78 if(iter!=map_pid.end())
79 {
80 daughters = iter->second;
81 m_trkIndex = (*it)->trackIndex();
82 os.setf(ios::left);
83 os<<daughters<<"["<<m_trkIndex<<"]"<<" ";
84 daughters.erase();
85 }
86 else
87 cout<<"can not find the daughter particle in pdt_table"<<endl;
88 }
89 os<<endl;
90
91 for(it = begin;it != end;it++)
92 {
93 SmartRef<Event::McParticle> pdMc = (*it);
94 printTree(os,pdMc,m_OutputLevel,tab+1);
95 }
96 }
97 if((!decay)&&(pmcp->primaryParticle()))
98 {
99 m_pid = (*pmcp).particleProperty();
100 map<int,string>::iterator iterPar;
101 iterPar = map_pid.find(m_pid);
102 if (iterPar!=map_pid.end())
103 os<<endl
104 <<"********************************************************"<<endl
105 <<" in this tree, there's only one primary particle : "<<iterPar->second<<endl
106 <<"********************************************************"<<endl;
107 }
108}
109
110
111
112
113void PrintMcInfo::printPartInf(ofstream& os,Event::McParticle* pMc,int m_OutputLevel,int tab=0)
114{
115 bool decay = (pMc->daughterList()).size()==0?false:true;
116
117 SmartDataPtr<Event::MdcMcHitCol> mdcMcHitCol(eventSvc(), "/Event/MC/MdcMcHitCol");
118 SmartDataPtr<Event::TofMcHitCol> tofMcHitCol(eventSvc(), "/Event/MC/TofMcHitCol");
119 SmartDataPtr<Event::EmcMcHitCol> emcMcHitCol(eventSvc(), "/Event/MC/EmcMcHitCol");
120 SmartDataPtr<Event::MucMcHitCol> mucMcHitCol(eventSvc(), "/Event/MC/MucMcHitCol");
121
122 if((!decay)&&(pMc->primaryParticle()))
123 {
124 os.setf(ios::left);
125 map<int,string>::iterator iter_map;
126 iter_map = map_pid.find(pMc->particleProperty());
127 if(iter_map!=map_pid.end())
128 {
129 if(m_OutputLevel==1)
130 {
131 string name = iter_map->second;
132 os<<"["<<pMc->trackIndex()<<"]"<<name<<endl;
133 HepLorentzVector p4 = (pMc)->initialFourMomentum();
134 os<<" "<<"["<<setw(12)<<p4.vect().cosTheta()<<setw(12)<<p4.vect().phi()<<"] ";
135 os<<"["<<setw(12)<<p4.x()<<setw(12)<<p4.y()
136 <<setw(12)<<p4.z()<<" "<<setw(10)<<p4.t()<<"] ";
137
138 HepLorentzVector ipst = pMc->initialPosition();
139 os<<"["<<setw(12)<<ipst.x();
140 os<<setw(12)<<ipst.y();
141 os<<setw(10)<<ipst.z()<<"] ";
142
143 HepLorentzVector fpst = (pMc)->finalPosition();
144 os<<"["<<setw(12)<<fpst.x();
145 os<<setw(12)<<fpst.y();
146 os<<setw(10)<<fpst.z()<<"] ";
147 os<<endl;
148 }
149 if(m_OutputLevel==2)
150 {
151 int trkIdx = pMc->trackIndex();
152 PrintMcInfo::printHit(os,mdcMcHitCol,tofMcHitCol,emcMcHitCol,mucMcHitCol,trkIdx);
153 }
154
155 }
156 }
157 if(decay)
158 {
159 SmartRefVector<Event::McParticle>::const_iterator begin = (pMc->daughterList()).begin();
160 SmartRefVector<Event::McParticle>::const_iterator end = (pMc->daughterList()).end();
161 SmartRefVector<Event::McParticle>::const_iterator it;
162 for(it = begin;it != end;it++)
163 {
164 m_trkIndex = (*it)->trackIndex();
165 map<int,string>::iterator iter;
166 iter = map_pid.find((*it)->particleProperty());
167 if(iter!=map_pid.end())
168 {
169 daughters = iter->second;
170 if(m_OutputLevel ==1)
171 {
172 os.setf(ios::left);
173 os<<"["<<m_trkIndex<<"]"<<setw(20)<<daughters<<endl;
174
175 HepLorentzVector p4 = (*it)->initialFourMomentum();
176 os<<" "<<"["<<setw(12)<<p4.vect().cosTheta()<<setw(12)<<p4.vect().phi()<<"] ";
177 os<<"["<<setw(12)<<p4.x()<<setw(12)<<p4.y()
178 <<setw(12)<<p4.z()<<" "<<setw(10)<<p4.t()<<"] ";
179
180 //os<<"["<<setw(4)<<(*it)->vertexIndex0();
181 HepLorentzVector ipst = (*it)->initialPosition();
182 os<<"["<<setw(12)<<ipst.x();
183 os<<setw(12)<<ipst.y();
184 os<<setw(10)<<ipst.z()<<"] ";
185 //os<<setw(8)<<ipst.t()<<"] ";
186
187 //os<<"["<<setw(4)<<(*it)->vertexIndex1();
188 HepLorentzVector fpst = (*it)->finalPosition();
189 os<<"["<<setw(12)<<fpst.x();
190 os<<setw(12)<<fpst.y();
191 os<<setw(10)<<fpst.z()<<"] ";
192 //os<<setw(8)<<fpst.t()<<"]";
193 }
194 if(m_OutputLevel==2)
195 {
196 if(((*mdcMcHitCol).size()==0&&(*tofMcHitCol).size()==0&&(*emcMcHitCol).size()==0&&(*mucMcHitCol).size()==0)&&tab==0)
197 os<<" --------------------------------------------------------"<<endl
198 <<" mdcMcHitCol,tofMcHitCol,emcMcHitCol,mucMcHitCol all empty"<<endl
199 <<" --------------------------------------------------------"<<endl;
200 if(!((*mdcMcHitCol).size()==0&&(*tofMcHitCol).size()==0&&(*emcMcHitCol).size()==0&&(*mucMcHitCol).size()==0))
201
202 {
203 os.setf(ios::left);
204 os<<"["<<m_trkIndex<<"]"<<setw(12)<<daughters<<endl;
205 }
206
207 int trkIdx = (*it)->trackIndex();
208 PrintMcInfo::printHit(os,mdcMcHitCol,tofMcHitCol,emcMcHitCol,mucMcHitCol,trkIdx);
209 }
210 os<<endl;
211 }
212 else cout<<"can not find the daughter particle in pdt_table"<<endl;
213 }
214 daughters.erase();
215
216 for(it = begin;it != end;it++)
217 {
218 SmartRef<Event::McParticle> pdMc = (*it);
219 printPartInf(os,pdMc,m_OutputLevel,tab+1);
220 }
221 }
222}
223
224
225
226void PrintMcInfo::printHit(ofstream& os,Event::MdcMcHitCol& mdcCol,Event::TofMcHitCol& tofCol,Event::EmcMcHitCol& emcCol,Event::MucMcHitCol& mucCol,int& trk_Idx)
227{
228 if(!(mdcCol.size()==0&&tofCol.size()==0&&emcCol.size()==0&&mucCol.size()==0))
229 {
230 os<< '\t' <<setw(66)<<"MdcHit"
231 <<setw(25)<<"TofHit"
232 <<setw(17)<<"EmcHit"
233 <<"MucHit"<<endl;
234 os<< '\t' <<"(layer, wire, position(x,y,z)(mm) ,drift_d(mm), DeDx(MeV/m))"<<" "<<"(B/EC layer phi_module)"<<" "<<"(B/EC theta phi)"<<" "<<"(B/EC segment layer strip)"<<endl;
235 }
236 Event::MdcMcHitCol::const_iterator it_mdc = mdcCol.begin();
237 Event::TofMcHitCol::const_iterator it_tof = tofCol.begin();
238 Event::EmcMcHitCol::const_iterator it_emc = emcCol.begin();
239 Event::MucMcHitCol::const_iterator it_muc = mucCol.begin();
240 vector<Event::MdcMcHitCol::const_iterator> vmdc;
241 vector<Event::TofMcHitCol::const_iterator> vtof;
242 vector<Event::EmcMcHitCol::const_iterator> vemc;
243 vector<Event::MucMcHitCol::const_iterator> vmuc;
244 if(mdcCol.size() != 0)
245 {
246 for(;it_mdc != mdcCol.end();it_mdc++)
247 {
248 int trkIndexmdc = (*it_mdc)->getTrackIndex();
249 if(trkIndexmdc == trk_Idx)
250 {
251 vmdc.push_back(it_mdc);
252 }
253 }
254 }
255 if(tofCol.size() != 0)
256 {
257 for(;it_tof != tofCol.end();it_tof++)
258 {
259 int trkIndextof = (*it_tof)->getTrackIndex();
260 if(trkIndextof==trk_Idx) vtof.push_back(it_tof);
261 }
262 }
263
264 if(emcCol.size() != 0)
265 {
266 for(;it_emc != emcCol.end();it_emc++)
267 {
268 int trkIndexemc = (*it_emc)->getTrackIndex();
269 if(trkIndexemc==trk_Idx) vemc.push_back(it_emc);
270 }
271 }
272
273 if(mucCol.size() != 0)
274 {
275 for(;it_muc != mucCol.end();it_muc++)
276 {
277 int trkIndexmuc = (*it_muc)->getTrackIndex();
278 if(trkIndexmuc==trk_Idx) vmuc.push_back(it_muc);
279 }
280 }
281
282
283 for(int i=0;;i++)
284 {
285 bool bmdc=(i>vmdc.size())?true:false;
286 bool btof=(i>vtof.size())?true:false;
287 bool bemc=(i>vemc.size())?true:false;
288 bool bmuc=(i>vmuc.size())?true:false;
289 if((i>=vmdc.size())&&(i>=vtof.size())&&(i>=vemc.size())&&(i>=vmuc.size())) break;
290 os<< '\t' ;
291 if(vmdc.size()>0)
292 {
293 if(i<vmdc.size())
294 {
295 const Identifier ident_mdc = (*vmdc[i])->identify();
296 os<<setw(5)<<MdcID::layer(ident_mdc);
297 os<<setw(5)<<MdcID::wire(ident_mdc);
298 os<<"[ ";
299 os<<setw(10)<<(*vmdc[i])->getPositionX();
300 os<<setw(10)<<(*vmdc[i])->getPositionY();
301 os<<setw(10)<<(*vmdc[i])->getPositionZ();
302 os<<"] ";
303 os<<setw(10)<<(*vmdc[i])->getDriftDistance();
304 os<<setw(10)<<(*vmdc[i])->getDepositEnergy();
305 }else{for (int m = 0; m < 64; m++) os<<" ";}
306 }else{for (int m = 0; m < 64; m++) os<<" ";}
307
308
309 if(vtof.size()>0)
310 {
311 if(i<vtof.size())
312 {
313 const Identifier ident_tof =(*vtof[i])->identify();
314 if(TofID::barrel_ec(ident_tof)==1){os<<"B "<<" ";}
315 if(TofID::barrel_ec(ident_tof)==0){os<<"E_E"<<" ";}
316 if(TofID::barrel_ec(ident_tof)==2){os<<"W_E"<<" ";}
317 os<<setw(3)<<TofID::layer(ident_tof);
318 os<<setw(17)<<TofID::phi_module(ident_tof);
319 }else{for (int m = 0; m < 25; m++) os<<" ";}
320 }else{for (int m = 0; m < 25; m++) os<<" ";}
321
322
323 if(vemc.size()>0)
324 {
325 if(i<vemc.size())
326 {
327 const Identifier ident_emc = (*vemc[i])->identify();
328 if(EmcID::barrel_ec(ident_emc)==1){os<<"B "<<" ";}
329 if(EmcID::barrel_ec(ident_emc)==0){os<<"E_E"<<" ";}
330 if(EmcID::barrel_ec(ident_emc)==2){os<<"W_E"<<" ";}
331 os<<setw(5)<<EmcID::theta_module(ident_emc);
332 os<<setw(9)<<EmcID::phi_module(ident_emc);
333 }else{for (int m = 0; m < 19; m++) os<<" ";}
334 }else{for (int m = 0; m < 19; m++) os<<" ";}
335
336 if(vmuc.size()>0)
337 {
338 if(i<vmuc.size())
339 {
340 const Identifier ident_muc = (*vmuc[i])->identify();
341 if(MucID::barrel_ec(ident_muc)==1){os<<"B "<<" ";}
342 if(MucID::barrel_ec(ident_muc)==0){os<<"E_E"<<" ";}
343 if(MucID::barrel_ec(ident_muc)==2){os<<"W_E"<<" ";}
344 os<<setw(3)<<MucID::segment(ident_muc);
345 os<<setw(3)<<MucID::layer(ident_muc);
346 os<<setw(4)<<MucID::strip(ident_muc);
347 os<<endl;
348 }else{ os<<endl;}
349 }else{os<<endl;}
350
351
352 //i++;
353 }
354 vmdc.clear();
355 vtof.clear();
356 vemc.clear();
357 vmuc.clear();
358
359}
360
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:101
EmcID::barrel_ec
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
Definition EmcID.cxx:38
EmcID::theta_module
static unsigned int theta_module(const Identifier &id)
Definition EmcID.cxx:43
EmcID::phi_module
static unsigned int phi_module(const Identifier &id)
Definition EmcID.cxx:48
Event::McParticle::primaryParticle
bool primaryParticle() const
Retrieve whether this is a primary particle.
Definition McParticle.cxx:13
Event::McParticle::initialPosition
const HepLorentzVector & initialPosition() const
Retrieve pointer to the start, end vertex positions.
Definition McParticle.cxx:95
Event::McParticle::trackIndex
int trackIndex() const
Definition McParticle.h:131
Event::McParticle::daughterList
const SmartRefVector< McParticle > & daughterList() const
access the process name
Definition McParticle.h:164
Event::McParticle::particleProperty
StdHepId particleProperty() const
Retrieve particle property.
Definition McParticle.cxx:7
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
MucID::barrel_ec
static int barrel_ec(const Identifier &id)
Values of different levels.
Definition MucID.cxx:41
MucID::layer
static int layer(const Identifier &id)
Definition MucID.cxx:61
MucID::segment
static int segment(const Identifier &id)
Definition MucID.cxx:51
MucID::strip
static int strip(const Identifier &id)
Definition MucID.cxx:76
PrintMcInfo::printPartInf
void printPartInf(ofstream &, Event::McParticle *, int, int)
Definition McTruth.cxx:113
PrintMcInfo::printTree
void printTree(ofstream &, Event::McParticle *, int, int)
Definition McTruth.cxx:55
PrintMcInfo::printHit
void printHit(ofstream &, Event::MdcMcHitCol &, Event::TofMcHitCol &, Event::EmcMcHitCol &, Event::MucMcHitCol &, int &)
Definition McTruth.cxx:226
PrintMcInfo::printTitle
void printTitle(ofstream &os, int)
Definition McTruth.cxx:27
PrintMcInfo::mkmap
void mkmap()
Definition McTruth.cxx:4
TofID::phi_module
static int phi_module(const Identifier &id)
Definition TofID.cxx:117
TofID::barrel_ec
static int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
Definition TofID.cxx:95
TofID::layer
static int layer(const Identifier &id)
Definition TofID.cxx:109
Event::MdcMcHitCol
ObjectVector< MdcMcHit > MdcMcHitCol
Definition MdcMcHit.h:151
Event::MucMcHitCol
ObjectVector< MucMcHit > MucMcHitCol
Definition MucMcHit.h:89
Event::TofMcHitCol
ObjectVector< TofMcHit > TofMcHitCol
Definition TofMcHit.h:100
Event::EmcMcHitCol
ObjectVector< EmcMcHit > EmcMcHitCol
Definition EmcMcHit.h:132