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