| # Gaudi Framework, version v25r0 |

Home |  | Generated: Mon Feb 17 2014

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
ParticlePropertySvc.cpp
Go to the documentation of this file.
1 // ============================================================================
2 // Include files
3 // ============================================================================
4 // STD&STL
5 // ============================================================================
6 #include <cstdlib>
7 #include <cstring>
8 #include <fstream>
9 // ============================================================================
10 // GaudiKernel
11 // ============================================================================
12 #include "GaudiKernel/ISvcLocator.h"
13 #include "GaudiKernel/MsgStream.h"
14 #include "GaudiKernel/ParticleProperty.h"
15 #include "GaudiKernel/PhysicalConstants.h"
16 #include "GaudiKernel/IFileAccess.h"
17 #include "GaudiKernel/System.h"
18 // ============================================================================
19 //#include "GaudiKernel/ToStream.h"
20 // ============================================================================
21 // Local
22 // ============================================================================
23 #include "ParticlePropertySvc.h"
24 // ============================================================================
25 namespace Gaudi {
29 DECLARE_COMPONENT(ParticlePropertySvc)
30 // ============================================================================
42 // ============================================================================
43 ParticlePropertySvc::ParticlePropertySvc
44 ( const std::string& name ,
45 ISvcLocator* svc )
46 : base_class( name, svc )
47 // the default name of input data file
48 , m_filename ( "ParticleTable.txt" )
49 , m_other ()
50 , m_particles ()
51 // storages :
52 , m_vectpp ()
53 , m_idmap ()
54 , m_namemap ()
55 , m_stdhepidmap()
56 , m_pythiaidmap()
57 //
58 , m_owned ()
59 , m_replaced ()
60 , m_fileAccess (0)
61 {
63 // Redefine the default name:
64 if( System::getEnv("PARAMFILESROOT", m_filename) )
65 {
66 m_filename += "/data/ParticleTable.txt";
67 }
68 //
69 declareProperty ( "ParticlePropertiesFile" , m_filename ) ;
70 declareProperty ( "OtherFiles" , m_other ) ;
71 declareProperty ( "Particles" , m_particles ) ;
72 }
73 // ============================================================================
75 // ============================================================================
76 ParticlePropertySvc::~ParticlePropertySvc()
77 {
78 for ( Set::iterator i = m_owned.begin(); i != m_owned.end() ; ++i )
79 { if ( 0 != *i ) { delete *i ; } }
80 }
81 // ============================================================================
83 // ============================================================================
84 StatusCode ParticlePropertySvc::initialize()
85 {
86 StatusCode sc = Service::initialize();
87 if ( sc.isFailure() ) { return sc ; }
88
89 MsgStream log( msgSvc() , name() ) ;
90
91 sc = setProperties();
92 if ( sc.isFailure() )
93 {
94 log << MSG::ERROR << " Could not set the properties " << endmsg ;
95 return sc ;
96 }
97
98
99 sc = service("VFSSvc",m_fileAccess);
100 if ( sc.isFailure() )
101 {
102 log << MSG::ERROR << " Cannot retrieve the VFS service " << endmsg ;
103 return sc ;
104 }
105
106 sc = parse();
107 if ( sc.isFailure() )
108 {
109 log << MSG::ERROR << " Could not parse the file " << endmsg ;
110 return sc ;
111 }
112 if ( !m_particles.empty() )
113 {
114 sc = addParticles () ;
115 if ( sc.isFailure() )
116 {
117 log << MSG::ERROR << " Could not treat particles! " << endmsg ;
118 return sc ;
119 }
120 }
121 log << MSG::DEBUG << "ParticleProperties parsed successfully" << endmsg;
122
123 log << MSG::DEBUG << "Access properties" << endmsg;
124 // For debugging purposes print out the size of the internal maps
125 // particle name as key: all particles in database are present here
126 log << MSG::DEBUG << "NameMap size =" << m_namemap.size() << endmsg;
127 // Geant3 ID as key: all particles in database are present here
128 log << MSG::DEBUG << "GeantID Map size =" << m_idmap.size() << endmsg;
129 // StdHep ID as key: some particles have no valid StdHep ID
130 log << MSG::DEBUG << "StdHepID Map size =" << m_stdhepidmap.size()
131 << endmsg;
132 // Pythia ID as key: some particles are not defined in Pythia
133 log << MSG::DEBUG << "PythiaID Map size =" << m_pythiaidmap.size()
134 << endmsg;
135
136 if ( !m_replaced.empty() )
137 {
138 log << MSG::INFO
139 << "Properties have been redefined for "
140 << m_replaced.size() << " particles : "
141 << Gaudi::Utils::toString( m_replaced )
142 << endmsg ;
143 }
144
145 return StatusCode::SUCCESS ;
146 }
147 // =============================================================================
149 // =============================================================================
150 StatusCode ParticlePropertySvc::finalize()
151 {
152 if ( !m_other.empty() )
153 {
154 MsgStream log( msgSvc() , name() ) ;
155 log << MSG::INFO
156 << "Additional Properties have been read from files: "
157 << Gaudi::Utils::toString ( m_other )
158 << endmsg ;
159 }
160
161 if ( !m_replaced.empty() )
162 {
163 MsgStream log( msgSvc() , name() ) ;
164 log << MSG::ALWAYS
165 << "Properties have been redefined for "
166 << m_replaced.size() << " particles : "
167 << Gaudi::Utils::toString( m_replaced )
168 << endmsg ;
169 }
170
171 if (m_fileAccess) {
172 m_fileAccess->release();
173 m_fileAccess = 0;
174 }
175
177 return Service::finalize () ;
178 }
179 // =============================================================================
181 // =============================================================================
182 StatusCode ParticlePropertySvc::push_back
183 ( const std::string& particle ,
184 int geantId ,
185 int jetsetId ,
186 double charge ,
187 double mass ,
188 double tlife ,
189 const std::string& evtName ,
190 int pythiaId ,
191 double maxWidth )
192 {
193 ParticleProperty* pp = new ParticleProperty
194 ( particle , geantId , jetsetId ,
195 charge , mass , tlife ,
196 evtName , pythiaId , maxWidth ) ;
197 //
198 m_owned.insert ( pp ) ;
199 //
200 return push_back( pp );
201 }
202 // =============================================================================
204 // =============================================================================
205 StatusCode ParticlePropertySvc::push_back ( ParticleProperty* pp )
206 {
207 if ( 0 == pp ) { return StatusCode::FAILURE ; }
208 //
209 { // try to add into Geant(3)ID map
210 const int ID = pp->geantID() ;
211 // is this already in the map?
212 MapID::const_iterator ifind = m_idmap.find( ID ) ;
213 if ( m_idmap.end() != ifind && 0 != m_idmap[ ID ])
214 {
215 diff ( ifind->second , pp ) ;
216 m_replaced.insert( m_idmap[ ID ]->particle() ) ;
217 }
218 // put it into the map
219 m_idmap[ ID ] = pp ;
220 }
221 //
222 { // try to add into Name map
223 const std::string& particle = pp->particle() ;
224 // is this already in the map?
225 MapName::const_iterator ifind = m_namemap.find( particle ) ;
226 if ( m_namemap.end() != ifind && 0 != m_namemap[ particle ] )
227 {
228 diff ( ifind->second , pp ) ;
229 m_replaced.insert( m_namemap[ particle ]->particle() ) ;
230 }
231 // put it into the map
232 m_namemap[ particle ] = pp ;
233 }
234 //
235 // add to StdHep map only if StdHep ID different from zero and if
236 // not Cerenkov (StdHep ID = gamma)
237 if ( 0 != pp->jetsetID() && "Tcherenkov" != pp->particle() )
238 { // try to add into StdHepID map
239 const int ID = pp->jetsetID() ;
240 // is this already in the map?
241 MapStdHepID::const_iterator ifind = m_stdhepidmap.find( ID ) ;
242 if ( m_stdhepidmap.end() != ifind && 0 != m_stdhepidmap[ ID ])
243 {
244 diff ( ifind->second , pp ) ;
245 m_replaced.insert( m_stdhepidmap[ ID ]->particle() ) ;
246 }
247 // put it into the map
248 m_stdhepidmap[ ID ] = pp ;
249 }
250 //
251 // add to Pythia map only if Pythia ID is different from
252 // zero ( StdHep id is always different from zero in this case )
253 if ( 0 != pp->pythiaID() &&
254 0 != pp->jetsetID() &&
255 "Tcherenkov" != pp->particle() )
256 { // try to add into PythiaID map
257 const int ID = pp->pythiaID() ;
258 // is this already in the map?
259 MapPythiaID::const_iterator ifind = m_pythiaidmap.find( ID ) ;
260 if ( m_pythiaidmap.end() != ifind && 0 != m_pythiaidmap[ ID ])
261 {
262 diff ( ifind->second , pp ) ;
263 m_replaced.insert( m_pythiaidmap[ ID ]->particle() ) ;
264 }
265 // put it into the map
266 m_pythiaidmap[ ID ] = pp ;
267 }
268 //
269 return rebuild() ;
270 }
271 // =============================================================================
273 // =============================================================================
274 namespace
275 {
276 template <class MAP>
277 void _remove_ ( MAP& m , const ParticleProperty* pp )
278 {
279 typename MAP::iterator i = m.begin() ;
280 for ( ; m.end() != i ; ++i ) { if ( i->second == pp ) { break ; } }
281 if ( m.end() != i ) { m.erase ( i ) ; }
282 }
283 }
284 // ============================================================================
285 StatusCode ParticlePropertySvc::erase( const ParticleProperty* pp )
286 {
287 if ( 0 == pp ) { return StatusCode::FAILURE ; }
288
289 _remove_ ( m_idmap , pp ) ;
290 _remove_ ( m_namemap , pp ) ;
291 _remove_ ( m_stdhepidmap , pp ) ;
292 _remove_ ( m_pythiaidmap , pp ) ;
293 //
294 return rebuild() ;
295 }
296 // ============================================================================
298 // ============================================================================
299 StatusCode ParticlePropertySvc::parse()
300 {
301
302 // parse "the main" file
303 StatusCode sc = parse ( m_filename ) ;
304 if ( sc.isFailure() ) { return sc ; }
305
306 // parse "other" files
307 for ( Files::const_iterator file = m_other.begin() ;
308 m_other.end() != file ; ++file )
309 {
310 sc = parse ( *file ) ;
311 if ( sc.isFailure() ) { return sc ; }
312 }
313
314 // Now check that the file format was consistent with what parser
315 // expected
316 if ( m_namemap.size() == 0 )
317 {
318 MsgStream log( msgSvc(), name() );
319 log << MSG::ERROR
320 << "Format of input file inconsistent with what expected"
321 << " - Check you are using ParticleData.txt" << endmsg;
322 return StatusCode::FAILURE;
323 }
324
325 return sc;
326 }
327 // ============================================================================
328 StatusCode ParticlePropertySvc::parse( const std::string& file )
329 {
330 StatusCode sc = StatusCode::FAILURE;
331
332 MsgStream log( msgSvc(), name() );
333 char line[ 255 ];
334
335 std::auto_ptr<std::istream> infileptr;
336 if (m_fileAccess) infileptr = m_fileAccess->open(file);
337
338 if ( infileptr.get() == 0 )
339 {
340 log << MSG::ERROR << "Unable to open properties file : " << file
341 << endmsg;
342 return StatusCode::FAILURE ;
343 }
344
345 std::istream &infile = *infileptr;
346 sc = StatusCode::SUCCESS;
347 log << MSG::INFO
348 << "Opened particle properties file : " << file << endmsg;
349
350 while( infile )
351 {
352 // parse each line of the file (comment lines begin with # in the cdf
353 // file,
354 infile.getline( line, 255 );
355
356 if ( line[0] == '#' ) continue;
357
359 #ifdef WIN32
360 // Disable warning
361 // C4996: 'strtok': This function or variable may be unsafe.
362 #pragma warning(disable:4996)
363 #endif
364 std::string par, gid, jid, chg, mas, lif, evt, pyt, mwi ;
365 char* token = strtok( line, " " );
366 if ( token ) { par = token; token = strtok( NULL, " " );} else continue;
367 if ( token ) { gid = token; token = strtok( NULL, " " );} else continue;
368 if ( token ) { jid = token; token = strtok( NULL, " " );} else continue;
369 if ( token ) { chg = token; token = strtok( NULL, " " );} else continue;
370 if ( token ) { mas = token; token = strtok( NULL, " " );} else continue;
371 if ( token ) { lif = token; token = strtok( NULL, " " );} else continue;
372 if ( token ) { evt = token; token = strtok( NULL, " " );} else continue;
373 if ( token ) { pyt = token; token = strtok( NULL, " " );} else continue;
374 if ( token ) { mwi = token; token = strtok( NULL, " " );} else continue;
375 if ( token != NULL ) continue;
376
377 // In SICb cdf file mass and lifetime units are GeV and sec, specify it so
378 // that they are converted to Gaudi units (MeV and ns)
379 double mass = atof( mas.c_str() ) * Gaudi::Units::GeV;
380 double tlife = atof( lif.c_str() ) * Gaudi::Units::s;
381 long ljid = atoi( jid.c_str() );
382 long lgid = atoi( gid.c_str() );
383 long lpyt = atoi( pyt.c_str() ) ;
384 double mW = atof( mwi.c_str() ) * Gaudi::Units::GeV ;
385
386 // Change the particles that do not correspond to a pdg number
387 if ( ljid == 0 ) {
388 ljid = 10000000*lgid;
389 }
390
391 // add a particle property
392 sc = push_back( par, lgid, ljid,
393 atof( chg.c_str() ), mass, tlife, evt, lpyt, mW ) ;
394 if ( sc.isFailure() )
395 {
396 log << MSG::ERROR
397 << "Error from ParticlePropertySvc::push_back for particle='"
398 << par << "'" << endmsg ;
399 }
400
401 }
402
403 //infile.close();
404
405 return StatusCode::SUCCESS ;
406 }
407 // ============================================================================
413 // ============================================================================
414 const ParticleProperty*
415 ParticlePropertySvc::anti ( const ParticleProperty* pp ) const
416 {
417 if ( 0 == pp ) { return 0 ; }
418 const int ID = pp->pdgID() ;
419 const int antiID = -1 * ID ;
420 for ( const_iterator it = m_vectpp.begin() ; m_vectpp.end() != it ; ++it )
421 {
422 const ParticleProperty* ap = *it ;
423 if ( 0 == ap ) { continue ; } // CONTINUE
424 if ( antiID == ap->pdgID() ) { return ap ; } // RETURN
425 };
426 //
427 return pp ; // RETURN
428 }
429 // ============================================================================
434 // ============================================================================
435 StatusCode ParticlePropertySvc::setAntiParticles()
436 {
437 // initialize particle<-->antiParticle relations
438 for ( iterator ip = m_vectpp.begin() ; m_vectpp.end() != ip ; ++ip )
439 {
440 ParticleProperty* pp = *ip ;
441 if ( 0 == pp ) { continue ; } // CONTINUE
442 const ParticleProperty* ap = anti ( pp ) ;
443 if ( 0 != ap ) { pp->setAntiParticle( ap ) ; }
444 }
445 return StatusCode::SUCCESS ;
446 }
447 // ============================================================================
449 // ============================================================================
450 namespace
451 {
453 template <class MAP>
454 void _load_ ( MAP& m , ParticlePropertySvc::Set& result )
455 {
456 for ( typename MAP::iterator i = m.begin() ; m.end() != i ; ++i )
457 { result.insert ( i->second ); }
458 }
459 }
460 // ============================================================================
461 StatusCode ParticlePropertySvc::rebuild()
462 {
463 Set local ;
464 m_vectpp.clear() ;
465 m_vectpp.reserve ( m_idmap.size() + 100 ) ;
466 // load information from maps into the set
467 _load_ ( m_idmap , local ) ;
468 _load_ ( m_namemap , local ) ;
469 _load_ ( m_stdhepidmap , local ) ;
470 _load_ ( m_pythiaidmap , local ) ;
471 // load information from set to the linear container vector
472 for ( Set::iterator i = local.begin() ; local.end() != i ; ++i )
473 { m_vectpp.push_back ( *i ) ; }
474 return setAntiParticles() ;
475 }
476 // ============================================================================
477 // treat additional particles
478 // ============================================================================
479 StatusCode ParticlePropertySvc::addParticles()
480 {
481
482 MsgStream log ( msgSvc() , name() ) ;
483 // loop over all "explicit" particles
484 for ( Particles::const_iterator ip = m_particles.begin() ;
485 m_particles.end() != ip ; ++ip )
486 {
487 const std::string item = *ip ;
488 std::istringstream input( item ) ;
489 // get the name
490 std::string p_name ;
491 int p_geant ;
492 int p_jetset ;
493 double p_charge ;
494 double p_mass ;
495 double p_ltime ;
496 std::string p_evtgen ;
497 int p_pythia ;
498 double p_maxwid ;
499 if ( input
500 >> p_name
501 >> p_geant
502 >> p_jetset
503 >> p_charge
504 >> p_mass
505 >> p_ltime
506 >> p_evtgen
507 >> p_pythia
508 >> p_maxwid )
509 {
510 log << MSG::ALWAYS
511 << " Add/Modify the particle: "
512 << " name='" << p_name << "'"
513 << " geant=" << p_geant
514 << " jetset=" << p_jetset
515 << " charge=" << p_charge
516 << " mass=" << p_mass
517 << " ltime=" << p_ltime
518 << " evtgen='" << p_evtgen << "'"
519 << " pythia=" << p_pythia
520 << " maxwid=" << p_maxwid << endmsg ;
521 //
522 StatusCode sc = push_back
523 ( p_name ,
524 p_geant ,
525 p_jetset ,
526 p_charge ,
527 p_mass * Gaudi::Units::GeV ,
528 p_ltime * Gaudi::Units::s ,
529 p_evtgen ,
530 p_pythia ,
531 p_maxwid * Gaudi::Units::GeV ) ;
532 if ( sc.isFailure() ) { return sc ; } // RETURN
533 }
534 else
535 {
536 log << MSG::ERROR
537 << " could not parse '" << item << "'" << endmsg ;
538 return StatusCode::FAILURE ; // RETURN
539 }
540 }
541 //
542 return StatusCode::SUCCESS ;
543 }
544 // ============================================================================
545 #ifdef __ICC
546 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
547 // The comparison are meant
548 #pragma warning(push)
549 #pragma warning(disable:1572)
550 #endif
551 bool ParticlePropertySvc::diff
552 ( const ParticleProperty* o ,
553 const ParticleProperty* n ,
554 const MSG::Level l ) const
555 {
556 //
557 if ( o == n ) { return false ; }
558 //
559 MsgStream log ( msgSvc() , name() ) ;
560 log << l ;
561 if ( 0 == o || 0 == n )
562 {
563 log << MSG::WARNING << " ParticleProperty* point to NULL" << endmsg ;
564 return true ; // RETURN
565 }
566 //
567 bool result = false ;
568 if ( o -> particle () != n -> particle () )
569 {
570 result = true ;
571 log << " Name:'" << o -> particle () << "'/'" << n -> particle () << "'" ;
572 }
573 if ( o -> geantID () != n -> geantID () )
574 {
575 result = true ;
576 log << " G3ID:" << o -> geantID () << "/" << n -> geantID () << "'" ;
577 }
578 if ( o -> pdgID () != n -> pdgID () )
579 {
580 result = true ;
581 log << " PDGID:" << o -> pdgID () << "/" << n -> pdgID () << "'" ;
582 }
583 if ( o -> pythiaID () != n -> pythiaID () )
584 {
585 result = true ;
586 log << " PYID:" << o -> pythiaID () << "/" << n -> pythiaID () << "'" ;
587 }
588 if ( o -> charge () != n -> charge () )
589 {
590 result = true ;
591 log << " Q:" << o -> charge () << "/" << n -> charge () << "'" ;
592 }
593 if ( o -> mass () != n -> mass () )
594 {
595 result = true ;
596 log << " M:" << o -> mass () << "/" << n -> mass () << "'" ;
597 }
598 if ( o -> lifetime () != n -> lifetime () )
599 {
600 result = true ;
601 log << " T:" << o -> lifetime () << "/" << n -> lifetime () << "'" ;
602 }
603 if ( o -> evtGenName () != n -> evtGenName () )
604 {
605 result = true ;
606 log << " EvtGen:" << o -> evtGenName () << "/" << n -> evtGenName () << "'" ;
607 }
608 if ( o -> maxWidth () != n -> maxWidth () )
609 {
610 result = true ;
611 log << " WMAX:" << o -> maxWidth () << "/" << n -> maxWidth () << "'" ;
612 }
613 if ( result ) { log << endmsg ; }
614 //
615 return result ;
616 }
617
618 }
619 #ifdef __ICC
620 // re-enable icc remark #1572
621 #pragma warning(pop)
622 #endif
623 // ============================================================================
624 // The END
625 // ============================================================================
626