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