| # Gaudi Framework, version v23r0 |

Home |  | Generated: Mon Jan 30 2012

# ParticlePropertySvc.cpp

Go to the documentation of this file.

    00001 // ============================================================================
    00002 // Include files
    00003 // ============================================================================
    00004 // STD&STL
    00005 // ============================================================================
    00006 #include <cstdlib>
    00007 #include <cstring>
    00008 #include <fstream>
    00009 // ============================================================================
    00010 // GaudiKernel
    00011 // ============================================================================
    00012 #include "GaudiKernel/SvcFactory.h"
    00013 #include "GaudiKernel/ISvcLocator.h"
    00014 #include "GaudiKernel/MsgStream.h"
    00015 #include "GaudiKernel/ParticleProperty.h"
    00016 #include "GaudiKernel/PhysicalConstants.h"
    00017 #include "GaudiKernel/IFileAccess.h"
    00018 #include "GaudiKernel/System.h"
    00019 // ============================================================================
    00020 //#include "GaudiKernel/ToStream.h"
    00021 // ============================================================================
    00022 // Local
    00023 // ============================================================================
    00024 #include "ParticlePropertySvc.h"
    00025 // ============================================================================
    00026 namespace Gaudi {
    00030 DECLARE_SERVICE_FACTORY(ParticlePropertySvc)
    00031 // ============================================================================
    00043 // ============================================================================
    00044 ParticlePropertySvc::ParticlePropertySvc
    00045 ( const std::string& name ,
    00046   ISvcLocator*       svc  )
    00047   : base_class( name, svc )
    00048   // the default name of input data file
    00049   , m_filename ( "ParticleTable.txt" )
    00050   , m_other ()
    00051   , m_particles ()
    00052   // storages :
    00053   , m_vectpp  ()
    00054   , m_idmap   ()
    00055   , m_namemap ()
    00056   , m_stdhepidmap()
    00057   , m_pythiaidmap()
    00058   //
    00059   , m_owned ()
    00060   , m_replaced ()
    00061   , m_fileAccess (0)
    00062 {
    00064   // Redefine the default name:
    00065   if( System::getEnv("PARAMFILESROOT", m_filename) )
    00066   {
    00067     m_filename += "/data/ParticleTable.txt";
    00068   }
    00069   //
    00070   declareProperty ( "ParticlePropertiesFile" , m_filename  ) ;
    00071   declareProperty ( "OtherFiles"             , m_other     ) ;
    00072   declareProperty ( "Particles"              , m_particles ) ;
    00073 }
    00074 // ============================================================================
    00076 // ============================================================================
    00077 ParticlePropertySvc::~ParticlePropertySvc()
    00078 {
    00079   for ( Set::iterator i = m_owned.begin(); i != m_owned.end() ; ++i )
    00080   { if ( 0 != *i ) { delete *i ; } }
    00081 }
    00082 // ============================================================================
    00084 // ============================================================================
    00085 StatusCode ParticlePropertySvc::initialize()
    00086 {
    00087   StatusCode sc = Service::initialize();
    00088   if ( sc.isFailure() ) { return sc ; }
    00089
    00090   MsgStream log( msgSvc() , name() ) ;
    00091
    00092   sc = setProperties();
    00093   if ( sc.isFailure() )
    00094   {
    00095     log << MSG::ERROR << " Could not set the properties " << endmsg ;
    00096     return sc ;
    00097   }
    00098
    00099
    00100   sc = service("VFSSvc",m_fileAccess);
    00101   if ( sc.isFailure() )
    00102   {
    00103     log << MSG::ERROR << " Cannot retrieve the VFS service " << endmsg ;
    00104     return sc ;
    00105   }
    00106
    00107   sc = parse();
    00108   if ( sc.isFailure() )
    00109   {
    00110     log << MSG::ERROR << " Could not parse the file " << endmsg ;
    00111     return sc ;
    00112   }
    00113   if ( !m_particles.empty() )
    00114   {
    00115     sc = addParticles () ;
    00116     if ( sc.isFailure() )
    00117     {
    00118       log << MSG::ERROR << " Could not treat particles! " << endmsg ;
    00119       return sc ;
    00120     }
    00121   }
    00122   log << MSG::DEBUG << "ParticleProperties parsed successfully" << endmsg;
    00123
    00124   log << MSG::DEBUG << "Access properties" << endmsg;
    00125   // For debugging purposes print out the size of the internal maps
    00126   // particle name as key: all particles in database are present here
    00127   log << MSG::DEBUG << "NameMap size =" << m_namemap.size() << endmsg;
    00128   // Geant3 ID as key: all particles in database are present here
    00129   log << MSG::DEBUG << "GeantID Map size =" << m_idmap.size() << endmsg;
    00130   // StdHep ID as key: some particles have no valid StdHep ID
    00131   log << MSG::DEBUG << "StdHepID Map size =" << m_stdhepidmap.size()
    00132       << endmsg;
    00133   // Pythia ID as key: some particles are not defined in Pythia
    00134   log << MSG::DEBUG << "PythiaID Map size =" << m_pythiaidmap.size()
    00135       << endmsg;
    00136
    00137   if ( !m_replaced.empty() )
    00138   {
    00139     log << MSG::INFO
    00140         << "Properties have been redefined for "
    00141         << m_replaced.size() << " particles : "
    00142         << Gaudi::Utils::toString( m_replaced )
    00143         << endmsg ;
    00144   }
    00145
    00146   return StatusCode::SUCCESS ;
    00147 }
    00148 // =============================================================================
    00150 // =============================================================================
    00151 StatusCode ParticlePropertySvc::finalize()
    00152 {
    00153   if ( !m_other.empty() )
    00154   {
    00155     MsgStream log( msgSvc() , name() ) ;
    00156     log << MSG::INFO
    00157         << "Additional Properties have been read from files: "
    00158         << Gaudi::Utils::toString ( m_other )
    00159         << endmsg ;
    00160   }
    00161
    00162   if ( !m_replaced.empty() )
    00163   {
    00164     MsgStream log( msgSvc() , name() ) ;
    00165     log << MSG::ALWAYS
    00166         << "Properties have been redefined for "
    00167         << m_replaced.size() << " particles : "
    00168         << Gaudi::Utils::toString( m_replaced )
    00169         << endmsg ;
    00170   }
    00171
    00172   if (m_fileAccess) {
    00173     m_fileAccess->release();
    00174     m_fileAccess = 0;
    00175   }
    00176
    00178   return Service::finalize () ;
    00179 }
    00180 // =============================================================================
    00182 // =============================================================================
    00183 StatusCode ParticlePropertySvc::push_back
    00184 ( const std::string& particle ,
    00185   int                geantId  ,
    00186   int                jetsetId ,
    00187   double             charge   ,
    00188   double             mass     ,
    00189   double             tlife    ,
    00190   const std::string& evtName  ,
    00191   int                pythiaId ,
    00192   double             maxWidth )
    00193 {
    00194   ParticleProperty* pp = new ParticleProperty
    00195     ( particle , geantId  , jetsetId ,
    00196       charge   , mass     , tlife    ,
    00197       evtName  , pythiaId , maxWidth ) ;
    00198   //
    00199   m_owned.insert ( pp ) ;
    00200   //
    00201   return push_back( pp );
    00202 }
    00203 // =============================================================================
    00205 // =============================================================================
    00206 StatusCode ParticlePropertySvc::push_back ( ParticleProperty* pp )
    00207 {
    00208   if ( 0 == pp ) { return StatusCode::FAILURE ; }
    00209   //
    00210   { // try to add into Geant(3)ID map
    00211     const int ID = pp->geantID() ;
    00212     // is this already in the map?
    00213     MapID::const_iterator ifind = m_idmap.find( ID ) ;
    00214     if ( m_idmap.end() != ifind && 0 != m_idmap[ ID ])
    00215     {
    00216       diff ( ifind->second , pp ) ;
    00217       m_replaced.insert( m_idmap[ ID ]->particle() ) ;
    00218     }
    00219     // put it into the map
    00220     m_idmap[ ID ] = pp ;
    00221   }
    00222   //
    00223   { // try to add into Name map
    00224     const std::string& particle = pp->particle() ;
    00225     // is this already in the map?
    00226     MapName::const_iterator ifind = m_namemap.find( particle ) ;
    00227     if ( m_namemap.end() != ifind && 0 != m_namemap[ particle ] )
    00228     {
    00229       diff ( ifind->second , pp ) ;
    00230       m_replaced.insert( m_namemap[ particle ]->particle() ) ;
    00231     }
    00232     // put it into the map
    00233     m_namemap[ particle ] = pp ;
    00234   }
    00235   //
    00236   // add to StdHep map only if StdHep ID different from zero and if
    00237   // not Cerenkov (StdHep ID = gamma)
    00238   if ( 0 != pp->jetsetID() && "Tcherenkov" != pp->particle() )
    00239   { // try to add into StdHepID map
    00240     const int ID = pp->jetsetID() ;
    00241     // is this already in the map?
    00242     MapStdHepID::const_iterator ifind = m_stdhepidmap.find( ID ) ;
    00243     if ( m_stdhepidmap.end() != ifind && 0 != m_stdhepidmap[ ID ])
    00244     {
    00245       diff ( ifind->second , pp ) ;
    00246       m_replaced.insert( m_stdhepidmap[ ID ]->particle() ) ;
    00247     }
    00248     // put it into the map
    00249     m_stdhepidmap[ ID ] = pp ;
    00250   }
    00251   //
    00252   // add to Pythia map only if Pythia ID is different from
    00253   // zero ( StdHep id is always different from zero in this case )
    00254   if ( 0 != pp->pythiaID() &&
    00255        0 != pp->jetsetID() &&
    00256        "Tcherenkov" != pp->particle() )
    00257   { // try to add into PythiaID map
    00258     const int ID = pp->pythiaID() ;
    00259     // is this already in the map?
    00260     MapPythiaID::const_iterator ifind = m_pythiaidmap.find( ID ) ;
    00261     if ( m_pythiaidmap.end() != ifind && 0 != m_pythiaidmap[ ID ])
    00262     {
    00263       diff ( ifind->second , pp ) ;
    00264       m_replaced.insert( m_pythiaidmap[ ID ]->particle() ) ;
    00265     }
    00266     // put it into the map
    00267     m_pythiaidmap[ ID ] = pp ;
    00268   }
    00269   //
    00270   return rebuild() ;
    00271 }
    00272 // =============================================================================
    00274 // =============================================================================
    00275 namespace
    00276 {
    00277   template <class MAP>
    00278   void _remove_ ( MAP& m , const ParticleProperty* pp )
    00279   {
    00280     typename MAP::iterator i = m.begin() ;
    00281     for ( ; m.end() != i ; ++i ) { if ( i->second == pp ) { break ; } }
    00282     if  ( m.end() != i ) { m.erase ( i ) ; }
    00283   }
    00284 }
    00285 // ============================================================================
    00286 StatusCode ParticlePropertySvc::erase( const ParticleProperty* pp )
    00287 {
    00288   if ( 0 == pp ) { return StatusCode::FAILURE ; }
    00289
    00290   _remove_ ( m_idmap       , pp ) ;
    00291   _remove_ ( m_namemap     , pp ) ;
    00292   _remove_ ( m_stdhepidmap , pp ) ;
    00293   _remove_ ( m_pythiaidmap , pp ) ;
    00294   //
    00295   return rebuild() ;
    00296 }
    00297 // ============================================================================
    00299 // ============================================================================
    00300 StatusCode ParticlePropertySvc::parse()
    00301 {
    00302
    00303   // parse "the main" file
    00304   StatusCode sc = parse ( m_filename ) ;
    00305   if ( sc.isFailure() ) { return sc ; }
    00306
    00307   // parse "other" files
    00308   for ( Files::const_iterator file = m_other.begin() ;
    00309         m_other.end() != file ; ++file )
    00310   {
    00311     sc = parse ( *file ) ;
    00312     if ( sc.isFailure() ) { return sc ; }
    00313   }
    00314
    00315   // Now check that the file format was consistent with what parser
    00316   // expected
    00317   if ( m_namemap.size() == 0 )
    00318   {
    00319     MsgStream log( msgSvc(), name() );
    00320     log << MSG::ERROR
    00321         << "Format of input file inconsistent with what expected"
    00322         << " - Check you are using ParticleData.txt" << endmsg;
    00323     return StatusCode::FAILURE;
    00324   }
    00325
    00326   return sc;
    00327 }
    00328 // ============================================================================
    00329 StatusCode ParticlePropertySvc::parse( const std::string& file )
    00330 {
    00331   StatusCode sc = StatusCode::FAILURE;
    00332
    00333   MsgStream log( msgSvc(), name() );
    00334   char line[ 255 ];
    00335
    00336   std::auto_ptr<std::istream> infileptr;
    00337   if (m_fileAccess) infileptr = m_fileAccess->open(file);
    00338
    00339   if ( infileptr.get() == 0 )
    00340   {
    00341     log << MSG::ERROR << "Unable to open properties file : " << file
    00342         << endmsg;
    00343     return StatusCode::FAILURE ;
    00344   }
    00345
    00346   std::istream &infile = *infileptr;
    00347   sc = StatusCode::SUCCESS;
    00348   log << MSG::INFO
    00349       << "Opened particle properties file : " << file << endmsg;
    00350
    00351   while( infile )
    00352   {
    00353     // parse each line of the file (comment lines begin with # in the cdf
    00354     // file,
    00355     infile.getline( line, 255 );
    00356
    00357     if ( line[0] == '#' ) continue;
    00358
    00360 #ifdef WIN32
    00361 // Disable warning
    00362 //   C4996: 'strtok': This function or variable may be unsafe.
    00363 #pragma warning(disable:4996)
    00364 #endif
    00365     std::string par, gid, jid, chg, mas, lif, evt, pyt, mwi ;
    00366     char* token = strtok( line, " " );
    00367     if ( token ) { par = token; token = strtok( NULL, " " );} else continue;
    00368     if ( token ) { gid = token; token = strtok( NULL, " " );} else continue;
    00369     if ( token ) { jid = token; token = strtok( NULL, " " );} else continue;
    00370     if ( token ) { chg = token; token = strtok( NULL, " " );} else continue;
    00371     if ( token ) { mas = token; token = strtok( NULL, " " );} else continue;
    00372     if ( token ) { lif = token; token = strtok( NULL, " " );} else continue;
    00373     if ( token ) { evt = token; token = strtok( NULL, " " );} else continue;
    00374     if ( token ) { pyt = token; token = strtok( NULL, " " );} else continue;
    00375     if ( token ) { mwi = token; token = strtok( NULL, " " );} else continue;
    00376     if ( token != NULL ) continue;
    00377
    00378     // In SICb cdf file mass and lifetime units are GeV and sec, specify it so
    00379     // that they are converted to Gaudi units (MeV and ns)
    00380     double mass = atof( mas.c_str() ) * Gaudi::Units::GeV;
    00381     double tlife = atof( lif.c_str() ) * Gaudi::Units::s;
    00382     long   ljid = atoi( jid.c_str() );
    00383     long   lgid = atoi( gid.c_str() );
    00384     long   lpyt = atoi( pyt.c_str() ) ;
    00385     double mW = atof( mwi.c_str() ) * Gaudi::Units::GeV ;
    00386
    00387     // Change the particles that do not correspond to a pdg number
    00388     if ( ljid == 0 ) {
    00389       ljid = 10000000*lgid;
    00390     }
    00391
    00392     // add a particle property
    00393     sc = push_back( par, lgid, ljid,
    00394                     atof( chg.c_str() ), mass, tlife, evt, lpyt, mW ) ;
    00395     if ( sc.isFailure() )
    00396     {
    00397       log << MSG::ERROR
    00398           << "Error from ParticlePropertySvc::push_back for particle='"
    00399           << par << "'" << endmsg ;
    00400     }
    00401
    00402   }
    00403
    00404   //infile.close();
    00405
    00406   return StatusCode::SUCCESS ;
    00407 }
    00408 // ============================================================================
    00414 // ============================================================================
    00415 const ParticleProperty*
    00416 ParticlePropertySvc::anti ( const ParticleProperty* pp ) const
    00417 {
    00418   if ( 0 == pp ) { return 0 ; }
    00419   const int     ID = pp->pdgID() ;
    00420   const int antiID = -1 * ID     ;
    00421   for ( const_iterator it = m_vectpp.begin() ; m_vectpp.end() != it ; ++it )
    00422   {
    00423     const ParticleProperty* ap = *it ;
    00424     if ( 0 == ap               ) { continue  ; }                // CONTINUE
    00425     if ( antiID == ap->pdgID() ) { return ap ; }                // RETURN
    00426   };
    00427   //
    00428   return pp ;                                                   // RETURN
    00429 }
    00430 // ============================================================================
    00435 // ============================================================================
    00436 StatusCode ParticlePropertySvc::setAntiParticles()
    00437 {
    00438   // initialize particle<-->antiParticle relations
    00439   for ( iterator ip = m_vectpp.begin() ; m_vectpp.end() != ip ; ++ip )
    00440   {
    00441     ParticleProperty* pp = *ip ;
    00442     if ( 0 == pp                    ) { continue ; }   // CONTINUE
    00443     const ParticleProperty* ap = anti ( pp ) ;
    00444     if ( 0 != ap                    ) { pp->setAntiParticle( ap ) ; }
    00445   }
    00446   return StatusCode::SUCCESS ;
    00447 }
    00448 // ============================================================================
    00450 // ============================================================================
    00451 namespace
    00452 {
    00454   template <class MAP>
    00455   void _load_ ( MAP& m , ParticlePropertySvc::Set& result )
    00456   {
    00457     for ( typename MAP::iterator i = m.begin() ; m.end() != i ; ++i )
    00458     { result.insert ( i->second ); }
    00459   }
    00460 }
    00461 // ============================================================================
    00462 StatusCode ParticlePropertySvc::rebuild()
    00463 {
    00464   Set local ;
    00465   m_vectpp.clear() ;
    00466   m_vectpp.reserve ( m_idmap.size() + 100 ) ;
    00467   // load information from maps into the set
    00468   _load_ ( m_idmap       , local ) ;
    00469   _load_ ( m_namemap     , local ) ;
    00470   _load_ ( m_stdhepidmap , local ) ;
    00471   _load_ ( m_pythiaidmap , local ) ;
    00472   // load information from set to the linear container vector
    00473   for ( Set::iterator i = local.begin() ; local.end() != i ; ++i )
    00474   { m_vectpp.push_back ( *i ) ; }
    00475   return setAntiParticles() ;
    00476 }
    00477 // ============================================================================
    00478 // treat additional particles
    00479 // ============================================================================
    00480 StatusCode ParticlePropertySvc::addParticles()
    00481 {
    00482
    00483   MsgStream log ( msgSvc() , name() ) ;
    00484   // loop over all "explicit" particles
    00485   for ( Particles::const_iterator ip = m_particles.begin() ;
    00486         m_particles.end() != ip ; ++ip )
    00487   {
    00488     const std::string item = *ip ;
    00489     std::istringstream input( item ) ;
    00490     // get the name
    00491     std::string p_name   ;
    00492     int         p_geant  ;
    00493     int         p_jetset ;
    00494     double      p_charge ;
    00495     double      p_mass   ;
    00496     double      p_ltime  ;
    00497     std::string p_evtgen ;
    00498     int         p_pythia ;
    00499     double      p_maxwid ;
    00500     if ( input
    00501          >> p_name
    00502          >> p_geant
    00503          >> p_jetset
    00504          >> p_charge
    00505          >> p_mass
    00506          >> p_ltime
    00507          >> p_evtgen
    00508          >> p_pythia
    00509          >> p_maxwid )
    00510     {
    00511       log << MSG::ALWAYS
    00512           << " Add/Modify the particle: "
    00513           << " name='"   << p_name   << "'"
    00514           << " geant="   << p_geant
    00515           << " jetset="  << p_jetset
    00516           << " charge="  << p_charge
    00517           << " mass="    << p_mass
    00518           << " ltime="   << p_ltime
    00519           << " evtgen='" << p_evtgen << "'"
    00520           << " pythia="  << p_pythia
    00521           << " maxwid="  << p_maxwid << endmsg ;
    00522       //
    00523       StatusCode sc = push_back
    00524         ( p_name                        ,
    00525           p_geant                       ,
    00526           p_jetset                      ,
    00527           p_charge                      ,
    00528           p_mass    * Gaudi::Units::GeV ,
    00529           p_ltime   * Gaudi::Units::s   ,
    00530           p_evtgen                      ,
    00531           p_pythia                      ,
    00532           p_maxwid  * Gaudi::Units::GeV ) ;
    00533       if ( sc.isFailure() ) { return sc ; }                        // RETURN
    00534     }
    00535     else
    00536     {
    00537       log << MSG::ERROR
    00538           << " could not parse '" << item << "'" << endmsg ;
    00539       return StatusCode::FAILURE ;                                 // RETURN
    00540     }
    00541   }
    00542   //
    00543   return StatusCode::SUCCESS ;
    00544 }
    00545 // ============================================================================
    00546 #ifdef __ICC
    00547 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
    00548 //   The comparison are meant
    00549 #pragma warning(push)
    00550 #pragma warning(disable:1572)
    00551 #endif
    00552 bool ParticlePropertySvc::diff
    00553 ( const ParticleProperty* o ,
    00554   const ParticleProperty* n ,
    00555   const MSG::Level        l ) const
    00556 {
    00557   //
    00558   if ( o == n ) { return false ; }
    00559   //
    00560   MsgStream log ( msgSvc() , name() ) ;
    00561   log << l ;
    00562   if ( 0 == o || 0 == n  )
    00563   {
    00564     log << MSG::WARNING << " ParticleProperty* point to NULL" << endmsg ;
    00565     return true ;                                                    // RETURN
    00566   }
    00567   //
    00568   bool result = false ;
    00569   if ( o -> particle () != n -> particle () )
    00570   {
    00571     result = true ;
    00572     log << " Name:'"  << o -> particle () << "'/'" << n -> particle () << "'" ;
    00573   }
    00574   if ( o -> geantID  () != n -> geantID  () )
    00575   {
    00576     result = true ;
    00577     log << " G3ID:"   << o -> geantID  () << "/"   << n -> geantID  () << "'" ;
    00578   }
    00579   if ( o -> pdgID    () != n -> pdgID    () )
    00580   {
    00581     result = true ;
    00582     log << " PDGID:"  << o -> pdgID    () << "/"   << n -> pdgID    () << "'" ;
    00583   }
    00584   if ( o -> pythiaID () != n -> pythiaID () )
    00585   {
    00586     result = true ;
    00587     log << " PYID:"   << o -> pythiaID () << "/"   << n -> pythiaID () << "'" ;
    00588   }
    00589   if ( o -> charge   () != n -> charge   () )
    00590   {
    00591     result = true ;
    00592     log << " Q:"      << o -> charge   () << "/"   << n -> charge   () << "'" ;
    00593   }
    00594   if ( o -> mass     () != n -> mass     () )
    00595   {
    00596     result = true ;
    00597     log << " M:"      << o -> mass     () << "/"   << n -> mass     () << "'" ;
    00598   }
    00599   if ( o -> lifetime () != n -> lifetime () )
    00600   {
    00601     result = true ;
    00602     log << " T:"      << o -> lifetime () << "/"   << n -> lifetime () << "'" ;
    00603   }
    00604   if ( o -> evtGenName () != n -> evtGenName () )
    00605   {
    00606     result = true ;
    00607     log << " EvtGen:" << o -> evtGenName () << "/"   << n -> evtGenName () << "'" ;
    00608   }
    00609   if ( o -> maxWidth () != n -> maxWidth () )
    00610   {
    00611     result = true ;
    00612     log << " WMAX:"   << o -> maxWidth () << "/"   << n -> maxWidth () << "'" ;
    00613   }
    00614   if ( result ) { log << endmsg ; }
    00615   //
    00616   return result ;
    00617 }
    00618
    00619 }
    00620 #ifdef __ICC
    00621 // re-enable icc remark #1572
    00622 #pragma warning(pop)
    00623 #endif
    00624 // ============================================================================
    00625 // The END
    00626 // ============================================================================
    00627

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines