ParticlePropertySvc.cpp

Go to the documentation of this file.

// $Id: ParticlePropertySvc.cpp,v 1.17 2008/10/27 16:41:33 marcocle Exp $
// ============================================================================
// CVS tag $Name:  $ , version $Revision: 1.17 $
// ============================================================================
//Include files
// ============================================================================
// STD&STL 
// ============================================================================
#include <cstdlib>
#include <cstring>
#include <fstream>
// ============================================================================
// GaudiKernel
// ============================================================================
#include "GaudiKernel/SvcFactory.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/ParticleProperty.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiKernel/IFileAccess.h"
// ============================================================================
//#include "GaudiKernel/ToStream.h"
// ============================================================================
// Local
// ============================================================================
#include "ParticlePropertySvc.h"

// ============================================================================
DECLARE_SERVICE_FACTORY(ParticlePropertySvc)
// ============================================================================
// ============================================================================
ParticlePropertySvc::ParticlePropertySvc 
( const std::string& name ,
  ISvcLocator*       svc  )
  : Service( name, svc ) 
  // the default name of input data file 
  , m_filename ( "ParticleTable.txt" ) 
  , m_other ()
  , m_particles ()
  // storages :
  , m_vectpp  () 
  , m_idmap   () 
  , m_namemap () 
  , m_stdhepidmap() 
  , m_pythiaidmap() 
  //
  , m_owned ()
  , m_replaced () 
  , m_fileAccess (0)
{
  // Redefine the default name: 
  if( NULL != getenv("PARAMFILESROOT") ) 
  {
    m_filename  = getenv( "PARAMFILESROOT" ) ;
    m_filename += "/data/ParticleTable.txt"  ;
  }
  //
  declareProperty ( "ParticlePropertiesFile" , m_filename  ) ;
  declareProperty ( "OtherFiles"             , m_other     ) ;
  declareProperty ( "Particles"              , m_particles ) ;
}
// ============================================================================
// ============================================================================
ParticlePropertySvc::~ParticlePropertySvc() 
{
  for ( Set::iterator i = m_owned.begin(); i != m_owned.end() ; ++i ) 
  { if ( 0 != *i ) { delete *i ; } }
} 
// ============================================================================
// ============================================================================
const InterfaceID& ParticlePropertySvc::type() const 
{ return IID_IParticlePropertySvc; }
// ============================================================================
// ============================================================================
StatusCode ParticlePropertySvc::initialize() 
{
  StatusCode sc = Service::initialize();
  if ( sc.isFailure() ) { return sc ; }
  
  MsgStream log( msgSvc() , name() ) ;
  
  sc = setProperties();  
  if ( sc.isFailure() ) 
  {
    log << MSG::ERROR << " Could not set the properties " << endreq ;
    return sc ;
  }
  
  
  sc = service("VFSSvc",m_fileAccess);
  if ( sc.isFailure() ) 
  {
    log << MSG::ERROR << " Cannot retrieve the VFS service " << endreq ;
    return sc ;
  }
  
  sc = parse();
  if ( sc.isFailure() ) 
  {
    log << MSG::ERROR << " Could not parse the file " << endreq ;
    return sc ;
  }
  if ( !m_particles.empty() ) 
  {
    sc = addParticles () ;
    if ( sc.isFailure() ) 
    {
      log << MSG::ERROR << " Could not treat particles! " << endreq ;
      return sc ;
    }
  }
  log << MSG::DEBUG << "ParticleProperties parsed successfully" << endreq;
  
  log << MSG::DEBUG << "Access properties" << endreq;
  // For debugging purposes print out the size of the internal maps  
  // particle name as key: all particles in database are present here
  log << MSG::DEBUG << "NameMap size =" << m_namemap.size() << endreq;
  // Geant3 ID as key: all particles in database are present here
  log << MSG::DEBUG << "GeantID Map size =" << m_idmap.size() << endreq;
  // StdHep ID as key: some particles have no valid StdHep ID 
  log << MSG::DEBUG << "StdHepID Map size =" << m_stdhepidmap.size()
      << endreq;
  // Pythia ID as key: some particles are not defined in Pythia
  log << MSG::DEBUG << "PythiaID Map size =" << m_pythiaidmap.size()
      << endreq;
  
  if ( !m_replaced.empty() ) 
  {
    log << MSG::INFO 
        << "Properties have beed redefined for "
        << m_replaced.size() << " particles : "
        << Gaudi::Utils::toString( m_replaced ) 
        << endreq ;
  }
  
  return StatusCode::SUCCESS ;
} 
// =============================================================================
// =============================================================================
StatusCode ParticlePropertySvc::finalize() 
{ 
  if ( !m_other.empty() ) 
  {
    MsgStream log( msgSvc() , name() ) ;
    log << MSG::INFO
        << "Additional Properties have beed read from files: "
        << Gaudi::Utils::toString ( m_other ) 
        << endreq ;
  }

  if ( !m_replaced.empty() ) 
  {
    MsgStream log( msgSvc() , name() ) ;
    log << MSG::ALWAYS 
        << "Properties have beed redefined for "
        << m_replaced.size() << " particles : "
        << Gaudi::Utils::toString( m_replaced ) 
        << endreq ;
  }

  if (m_fileAccess) {
    m_fileAccess->release();
    m_fileAccess = 0;
  }
  
  return Service::finalize () ; 
}
// =============================================================================
// =============================================================================
StatusCode ParticlePropertySvc::queryInterface
( const InterfaceID& riid, 
  void** ppvInterface ) 
{
  // chech the placeholder for result 
  if ( !ppvInterface ) { return StatusCode::FAILURE ; }
  //
  *ppvInterface = 0;
  //
  if ( riid == IID_IParticlePropertySvc ) 
  {
    *ppvInterface = static_cast<IParticlePropertySvc*>(this);
    addRef();                
    return StatusCode::SUCCESS ;    
  }
  //
  return  Service::queryInterface( riid, ppvInterface ) ; 
} 
// =============================================================================
// =============================================================================
StatusCode ParticlePropertySvc::push_back
( const std::string& particle ,
  int                geantId  , 
  int                jetsetId , 
  double             charge   , 
  double             mass     , 
  double             tlife    ,
  const std::string& evtName  , 
  int                pythiaId , 
  double             maxWidth ) 
{
  ParticleProperty* pp = new ParticleProperty 
    ( particle , geantId  , jetsetId , 
      charge   , mass     , tlife    , 
      evtName  , pythiaId , maxWidth ) ;
  // 
  m_owned.insert ( pp ) ;
  //
  return push_back( pp );
} 
// =============================================================================
// =============================================================================
StatusCode ParticlePropertySvc::push_back ( ParticleProperty* pp ) 
{
  if ( 0 == pp ) { return StatusCode::FAILURE ; }
  //
  { // try to add into Geant(3)ID map
    const int ID = pp->geantID() ;
    // is this already in the map?
    MapID::const_iterator ifind = m_idmap.find( ID ) ;
    if ( m_idmap.end() != ifind && 0 != m_idmap[ ID ]) 
    { 
      diff ( ifind->second , pp ) ;
      m_replaced.insert( m_idmap[ ID ]->particle() ) ; 
    }
    // put it into the map 
    m_idmap[ ID ] = pp ;
  }
  //
  { // try to add into Name map 
    const std::string& particle = pp->particle() ;
    // is this already in the map?
    MapName::const_iterator ifind = m_namemap.find( particle ) ;
    if ( m_namemap.end() != ifind && 0 != m_namemap[ particle ] ) 
    { 
      diff ( ifind->second , pp ) ;
      m_replaced.insert( m_namemap[ particle ]->particle() ) ; 
    }
    // put it into the map 
    m_namemap[ particle ] = pp ;    
  }
  //
  // add to StdHep map only if StdHep ID different from zero and if
  // not Cerenkov (StdHep ID = gamma)  
  if ( 0 != pp->jetsetID() && "Tcherenkov" != pp->particle() )
  { // try to add into StdHepID map
    const int ID = pp->jetsetID() ;
    // is this already in the map?
    MapStdHepID::const_iterator ifind = m_stdhepidmap.find( ID ) ;
    if ( m_stdhepidmap.end() != ifind && 0 != m_stdhepidmap[ ID ]) 
    { 
      diff ( ifind->second , pp ) ;
      m_replaced.insert( m_stdhepidmap[ ID ]->particle() ) ; 
    }
    // put it into the map 
    m_stdhepidmap[ ID ] = pp ;
  }
  //
  // add to Pythia map only if Pythia ID is different from 
  // zero ( StdHep id is always different from zero in this case )
  if ( 0 != pp->pythiaID() && 
       0 != pp->jetsetID() &&  
       "Tcherenkov" != pp->particle() )
  { // try to add into PythiaID map
    const int ID = pp->pythiaID() ;
    // is this already in the map?
    MapPythiaID::const_iterator ifind = m_pythiaidmap.find( ID ) ;
    if ( m_pythiaidmap.end() != ifind && 0 != m_pythiaidmap[ ID ]) 
    {
      diff ( ifind->second , pp ) ; 
      m_replaced.insert( m_pythiaidmap[ ID ]->particle() ) ; 
    }
    // put it into the map 
    m_pythiaidmap[ ID ] = pp ;
  }
  //
  return rebuild() ;
}
// =============================================================================
// =============================================================================
namespace 
{
  template <class MAP>
  void _remove_ ( MAP& m , const ParticleProperty* pp ) 
  {
    typename MAP::iterator i = m.begin() ;
    for ( ; m.end() != i ; ++i ) { if ( i->second == pp ) { break ; } }
    if  ( m.end() != i ) { m.erase ( i ) ; }
  } 
}
// ============================================================================
StatusCode ParticlePropertySvc::erase( const ParticleProperty* pp ) 
{
  if ( 0 == pp ) { return StatusCode::FAILURE ; }
  
  _remove_ ( m_idmap       , pp ) ;
  _remove_ ( m_namemap     , pp ) ;
  _remove_ ( m_stdhepidmap , pp ) ;
  _remove_ ( m_pythiaidmap , pp ) ;
  // 
  return rebuild() ;
} 
// ============================================================================
// ============================================================================
StatusCode ParticlePropertySvc::parse() 
{

  // parse "the main" file 
  StatusCode sc = parse ( m_filename ) ;
  if ( sc.isFailure() ) { return sc ; }
  
  // parse "other" files 
  for ( Files::const_iterator file = m_other.begin() ;
        m_other.end() != file ; ++file ) 
  {
    sc = parse ( *file ) ;
    if ( sc.isFailure() ) { return sc ; }
  }
  
  // Now check that the file format was consistent with what parser
  // expected
  if ( m_namemap.size() == 0 ) 
  {
    MsgStream log( msgSvc(), name() );
    log << MSG::ERROR 
        << "Format of input file inconsistent with what expected"
        << " - Check you are using ParticleData.txt" << endreq;
    return StatusCode::FAILURE;
  }
  
  return sc;
}
// ============================================================================
StatusCode ParticlePropertySvc::parse( const std::string& file ) 
{
  StatusCode sc = StatusCode::FAILURE;
  
  MsgStream log( msgSvc(), name() );
  char line[ 255 ];
  
  std::auto_ptr<std::istream> infileptr;
  if (m_fileAccess) infileptr = m_fileAccess->open(file);
  
  if ( infileptr.get() == 0 )
  {
    log << MSG::ERROR << "Unable to open properties file : " << file
        << endreq;
    return StatusCode::FAILURE ;
  }
  
  std::istream &infile = *infileptr;
  sc = StatusCode::SUCCESS;
  log << MSG::INFO 
      << "Opened particle properties file : " << file << endreq;
  
  while( infile ) 
  {
    // parse each line of the file (comment lines begin with # in the cdf
    // file,
    infile.getline( line, 255 );
    
    if ( line[0] == '#' ) continue;
    
    std::string par, gid, jid, chg, mas, lif, evt, pyt, mwi ; 
    char* token = strtok( line, " " );
    if ( token ) { par = token; token = strtok( NULL, " " );} else continue;
    if ( token ) { gid = token; token = strtok( NULL, " " );} else continue;
    if ( token ) { jid = token; token = strtok( NULL, " " );} else continue;
    if ( token ) { chg = token; token = strtok( NULL, " " );} else continue;
    if ( token ) { mas = token; token = strtok( NULL, " " );} else continue;
    if ( token ) { lif = token; token = strtok( NULL, " " );} else continue;
    if ( token ) { evt = token; token = strtok( NULL, " " );} else continue;
    if ( token ) { pyt = token; token = strtok( NULL, " " );} else continue;
    if ( token ) { mwi = token; token = strtok( NULL, " " );} else continue;
    if ( token != NULL ) continue;
    
    // In SICb cdf file mass and lifetime units are GeV and sec, specify it so
    // that they are converted to Gaudi units (MeV and ns)
    double mass = atof( mas.c_str() ) * Gaudi::Units::GeV;
    double tlife = atof( lif.c_str() ) * Gaudi::Units::s;
    long   ljid = atoi( jid.c_str() );
    long   lgid = atoi( gid.c_str() );
    long   lpyt = atoi( pyt.c_str() ) ;
    double mW = atof( mwi.c_str() ) * Gaudi::Units::GeV ;
    
    // Change the particles that do not corresopond to a pdg number
    if ( ljid == 0 ) {
      ljid = 10000000*lgid;
    }
    
    // add a particle property
    sc = push_back( par, lgid, ljid, 
                    atof( chg.c_str() ), mass, tlife, evt, lpyt, mW ) ;
    if ( sc.isFailure() )
    {
      log << MSG::ERROR 
          << "Error from ParticlePropertySvc::push_back for particle='" 
          << par << "'" << endreq ;
    }
    
  }
  
  //infile.close();
  
  return StatusCode::SUCCESS ;
}
// ============================================================================
// ============================================================================
const ParticleProperty* 
ParticlePropertySvc::anti ( const ParticleProperty* pp ) const 
{
  if ( 0 == pp ) { return 0 ; }
  const int     ID = pp->pdgID() ;
  const int antiID = -1 * ID     ;
  for ( const_iterator it = m_vectpp.begin() ; m_vectpp.end() != it ; ++it ) 
  {
    const ParticleProperty* ap = *it ;
    if ( 0 == ap               ) { continue  ; }                // CONTINUE 
    if ( antiID == ap->pdgID() ) { return ap ; }                // RETURN
  };
  //  
  return pp ;                                                   // RETURN 
}
// ============================================================================
// ============================================================================
StatusCode ParticlePropertySvc::setAntiParticles()
{ 
  // initialize particle<-->antiParticle relations 
  for ( iterator ip = m_vectpp.begin() ; m_vectpp.end() != ip ; ++ip ) 
  {
    ParticleProperty* pp = *ip ;
    if ( 0 == pp                    ) { continue ; }   // CONTINUE
    const ParticleProperty* ap = anti ( pp ) ;
    if ( 0 != ap                    ) { pp->setAntiParticle( ap ) ; }
  }
  return StatusCode::SUCCESS ;
}
// ============================================================================
// ============================================================================
namespace 
{
  template <class MAP>
  void _load_ ( MAP& m , ParticlePropertySvc::Set& result ) 
  {
    for ( typename MAP::iterator i = m.begin() ; m.end() != i ; ++i )
    { result.insert ( i->second ); }
  }
}
// ============================================================================
StatusCode ParticlePropertySvc::rebuild() 
{
  Set local ;
  m_vectpp.clear() ;
  m_vectpp.reserve ( m_idmap.size() + 100 ) ;
  // load information from maps into the set 
  _load_ ( m_idmap       , local ) ;
  _load_ ( m_namemap     , local ) ;
  _load_ ( m_stdhepidmap , local ) ;
  _load_ ( m_pythiaidmap , local ) ;
  // load information from set to the linear container vector 
  for ( Set::iterator i = local.begin() ; local.end() != i ; ++i )
  { m_vectpp.push_back ( *i ) ; }
  return setAntiParticles() ;
}
// ============================================================================
// treat additional particles 
// ============================================================================
StatusCode ParticlePropertySvc::addParticles() 
{

  MsgStream log ( msgSvc() , name() ) ;
  // loop over all "explicit" particles 
  for ( Particles::const_iterator ip = m_particles.begin() ; 
        m_particles.end() != ip ; ++ip ) 
  {
    const std::string item = *ip ;
    std::istringstream input( item ) ;
    // get the name
    std::string p_name   ;
    int         p_geant  ;
    int         p_jetset ;
    double      p_charge ;
    double      p_mass   ;
    double      p_ltime  ;
    std::string p_evtgen ;
    int         p_pythia ;
    double      p_maxwid ;
    if ( input 
         >> p_name
         >> p_geant 
         >> p_jetset 
         >> p_charge 
         >> p_mass 
         >> p_ltime 
         >> p_evtgen 
         >> p_pythia
         >> p_maxwid )
    {
      log << MSG::ALWAYS
          << " Add/Modify the particle: "
          << " name='"   << p_name   << "'"
          << " geant="   << p_geant   
          << " jetset="  << p_jetset  
          << " charge="  << p_charge  
          << " mass="    << p_mass     
          << " ltime="   << p_ltime   
          << " evtgen='" << p_evtgen << "'" 
          << " pythia="  << p_pythia  
          << " maxwid="  << p_maxwid << endreq ;
      //
      StatusCode sc = push_back
        ( p_name                        , 
          p_geant                       , 
          p_jetset                      ,
          p_charge                      , 
          p_mass    * Gaudi::Units::GeV , 
          p_ltime   * Gaudi::Units::s   , 
          p_evtgen                      ,
          p_pythia                      , 
          p_maxwid  * Gaudi::Units::GeV ) ;
      if ( sc.isFailure() ) { return sc ; }                        // RETURN 
    }
    else
    {
      log << MSG::ERROR
          << " could not parse '" << item << "'" << endreq ;
      return StatusCode::FAILURE ;                                 // RETURN 
    }
  }
  //
  return StatusCode::SUCCESS ;
}
// ============================================================================
bool ParticlePropertySvc::diff
( const ParticleProperty* o , 
  const ParticleProperty* n , 
  const MSG::Level        l ) const
{
  //
  if ( o == n ) { return false ; }
  //
  MsgStream log ( msgSvc() , name() ) ;
  log << l ;
  if ( 0 == o || 0 == n  ) 
  {
    log << MSG::WARNING << " ParticleProperty* point to NULL" << endreq ;
    return true ;                                                    // RETURN 
  }
  //
  bool result = false ;
  if ( o -> particle () != n -> particle () ) 
  { 
    result = true ;
    log << " Name:'"  << o -> particle () << "'/'" << n -> particle () << "'" ; 
  }
  if ( o -> geantID  () != n -> geantID  () ) 
  { 
    result = true ;
    log << " G3ID:"   << o -> geantID  () << "/"   << n -> geantID  () << "'" ; 
  }
  if ( o -> pdgID    () != n -> pdgID    () ) 
  { 
    result = true ;
    log << " PDGID:"  << o -> pdgID    () << "/"   << n -> pdgID    () << "'" ;
  }
  if ( o -> pythiaID () != n -> pythiaID () ) 
  { 
    result = true ;
    log << " PYID:"   << o -> pythiaID () << "/"   << n -> pythiaID () << "'" ; 
  }
  if ( o -> charge   () != n -> charge   () ) 
  { 
    result = true ;
    log << " Q:"      << o -> charge   () << "/"   << n -> charge   () << "'" ; 
  }
  if ( o -> mass     () != n -> mass     () ) 
  { 
    result = true ;
    log << " M:"      << o -> mass     () << "/"   << n -> mass     () << "'" ; 
  }
  if ( o -> lifetime () != n -> lifetime () ) 
  { 
    result = true ;
    log << " T:"      << o -> lifetime () << "/"   << n -> lifetime () << "'" ; 
  }
  if ( o -> evtGenName () != n -> evtGenName () ) 
  {
    result = true ;
    log << " EvtGen:" << o -> evtGenName () << "/"   << n -> evtGenName () << "'" ; 
  }
  if ( o -> maxWidth () != n -> maxWidth () ) 
  { 
    result = true ;
    log << " WMAX:"   << o -> maxWidth () << "/"   << n -> maxWidth () << "'" ; 
  }
  if ( result ) { log << endreq ; }
  //
  return result ;
}





// ============================================================================
// The END 
// ============================================================================