GenModule.cxx

Go to the documentation of this file.

// -------------------------------------------------------------
// File: GeneratorModule/GenModule.cxx
// Description:
//    This class is the base class used to specify the behavior of all
//    event generator modules and is meant to capture the common behavior
//    of these modules.
//
// Header for this module:-
 
#include <fstream>
 
#include "GeneratorModule/GenModule.h"
 
// Framework Related Headers:-
 
#include "GaudiKernel/DataSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include "GaudiKernel/IIncidentSvc.h"
#include "GaudiKernel/Incident.h"
 
#include "GaudiKernel/IPartPropSvc.h"
 
// Other Packages used by this class:-
#include "CLHEP/Random/RandPoisson.h"
#include "CLHEP/Random/Ranlux64Engine.h"
 
// #include "BesHepMC/GenVertex.h"
#include "HepMC/GenVertex.h"
 
#include "GeneratorObject/McGenEvent.h"
 
DECLARE_COMPONENT( GenModule )
//---------------------------------------------------------------------------
GenModule::GenModule( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator )
    , m_pRandomEngine( 0 )
    , m_pPoissonGenerator( 0 )
//---------------------------------------------------------------------------
{
  declareProperty( "FixedMode", m_fixedMode = true );
  declareProperty( "MeanInt", m_meanInteractions = 1.0 );
  declareProperty( "RandomSeed", m_randomSeed = 1234567 );
  declareProperty( "StripPartons", m_StripVector );
  m_StripVector.push_back( 0 );
}
 
//---------------------------------------------------------------------------
GenModule::~GenModule() {
  // Delete random number objects if they were created
 
  if ( m_pPoissonGenerator ) delete m_pPoissonGenerator;
  if ( m_pRandomEngine ) delete m_pRandomEngine;
}
//---------------------------------------------------------------------------
 
//---------------------------------------------------------------------------
StatusCode GenModule::initialize() {
  //---------------------------------------------------------------------------
 
  // Inform the user what the mode and conditions are:
  MsgStream log( msgSvc(), name() );
 
  // Initialize StripPartons variables
  if ( m_StripVector[0] > 0 ) { StripPartonsInit(); }
  else
  {
    m_StripPartonSwitch = false;
    m_StripVector.clear();
  }
 
  // Get the Particle Properties Service
  IPartPropSvc* p_PartPropSvc;
  // static const bool CREATEIFNOTTHERE(true);
  // StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
  StatusCode PartPropStatus = service( "PartPropSvc", p_PartPropSvc );
  if ( !PartPropStatus.isSuccess() || 0 == p_PartPropSvc )
  {
    log << MSG::ERROR << " Could not initialize Particle Properties Service" << endmsg;
    return PartPropStatus;
  }
 
  m_particleTable = p_PartPropSvc->PDT();
 
  m_pRandomEngine = new CLHEP::Ranlux64Engine( m_randomSeed );
 
  if ( m_fixedMode )
  {
    if ( m_meanInteractions == 1.0 )
    { log << MSG::INFO << "Standard Initialization:  Single Interaction Mode " << endmsg; }
    else
    {
      log << MSG::INFO << "Fixed Number of Interactions per Event is: " << m_meanInteractions
          << endmsg;
    }
  }
  else
  {
    m_pPoissonGenerator = new CLHEP::RandPoisson( *m_pRandomEngine, m_meanInteractions );
 
    log << MSG::INFO
        << "Poisson Distribution of Interactions per Event with Mean: " << m_meanInteractions
        << endmsg;
  }
  // Initialize the generator itself
  StatusCode status = this->genInitialize();
  if ( status.isFailure() )
  {
    log << MSG::ERROR << "Could not initialize Generator properly" << endmsg;
    return status;
  }
  StatusCode status1 = this->genuserInitialize();
  if ( status1.isFailure() )
  {
    log << MSG::ERROR << "Could not initialize user part properly" << endmsg;
    return status1;
  }
  return status;
}
 
//---------------------------------------------------------------------------
StatusCode GenModule::execute() {
  //---------------------------------------------------------------------------
 
  int        numToGenerate;
  StatusCode status;
 
  MsgStream log( msgSvc(), name() );
 
  log << MSG::DEBUG << "GenModule::execute()" << endmsg;
 
  // Decide how many interactions to generate
  if ( m_fixedMode ) { numToGenerate = (int)m_meanInteractions; }
  else { numToGenerate = m_pPoissonGenerator->fire(); }
 
  // Generate as many interactions as requested
  for ( int i = 0; i < numToGenerate; i++ )
  {
 
    // Call the code that generates an event
    status = this->callGenerator();
 
    // Create the MC event and send the GeneratorEvent stored in it to fillEvt
    GenEvent* evt = new GenEvent( 1, 1 );
    status        = this->fillEvt( evt );
 
    // Strip out the requested partons here.
    if ( m_StripPartonSwitch ) StripPartons( evt );
 
    // Check if the McCollection already exists
    SmartDataPtr<McGenEventCol> anMcCol( eventSvc(), "/Event/Gen" );
    if ( anMcCol != 0 )
    {
      // Add event to existing collection
      MsgStream log( msgSvc(), name() );
      log << MSG::INFO << "Add McGenEvent to existing collection" << endmsg;
      McGenEvent* mcEvent = new McGenEvent( evt );
      anMcCol->push_back( mcEvent );
    }
    else
    {
      // Create Collection and add  to the transient store
      McGenEventCol* mcColl  = new McGenEventCol;
      McGenEvent*    mcEvent = new McGenEvent( evt );
      mcColl->push_back( mcEvent );
 
      StatusCode sc = eventSvc()->registerObject( "/Event/Gen", mcColl );
 
      if ( sc != StatusCode::SUCCESS )
      {
        MsgStream log( msgSvc(), name() );
        log << MSG::ERROR << "Could not register McGenEvent" << endmsg;
        delete mcColl;
        delete evt;
        delete mcEvent;
        return StatusCode::FAILURE;
      }
      else
      {
        //        std::cout << " McGenEventCol made and stored" << std::endl;
      }
    }
  }
 
  // now call the incident service and set the signal.
  IIncidentSvc* incSvc;
  service( "IncidentSvc", incSvc );
  incSvc->fireIncident( Incident( name(), "McGenEvent Generated" ) );
 
  return status;
}
 
//---------------------------------------------------------------------------
StatusCode GenModule::finalize() {
  //---------------------------------------------------------------------------
 
  StatusCode status = this->genFinalize();
 
  return status;
}
 
//---------------------------------------------------------------------------
StatusCode GenModule::genInitialize() {
  //---------------------------------------------------------------------------
  return StatusCode::SUCCESS;
}
 
//---------------------------------------------------------------------------
StatusCode GenModule::genuserInitialize() {
  //---------------------------------------------------------------------------
  return StatusCode::SUCCESS;
}
 
//---------------------------------------------------------------------------
StatusCode GenModule::callGenerator() {
  //---------------------------------------------------------------------------
  return StatusCode::SUCCESS;
}
 
//---------------------------------------------------------------------------
StatusCode GenModule::genFinalize() {
  //---------------------------------------------------------------------------
  return StatusCode::SUCCESS;
}
 
//---------------------------------------------------------------------------
StatusCode GenModule::fillEvt( GenEvent* /* evt */ ) {
  //---------------------------------------------------------------------------
  return StatusCode::SUCCESS;
}
 
void GenModule::StripPartonsInit( void ) {
  MsgStream log( msgSvc(), name() );
 
  for ( int i = 1; i < 9; ++i )
  {
    m_AllPartons.push_back( i );
    m_AllPartons.push_back( -i );
  }                             // Quarks
  m_AllPartons.push_back( 21 ); // gluon
 
  m_StripPartonSwitch = true;
  m_StripVector.erase( m_StripVector.begin(), m_StripVector.begin() + 1 );
  // When the user specifies only the StripPartonSwitch or gives a Particle Code
  // which is NOT a parton then strip ALL partons
  if ( m_StripVector.size() == 0 )
  {
    log << MSG::INFO << " !!!! NO SPECIFIC PARTON FOR STRIPOUT WAS REQUESTED => STRIPOUT ALL "
        << endmsg;
    m_StripVector = m_AllPartons;
  }
  else
  {
    bool                             value_ok = true;
    std::vector<int>::const_iterator i        = m_StripVector.begin();
    do {
      if ( std::find( m_AllPartons.begin(), m_AllPartons.end(), *i ) == m_AllPartons.end() )
        value_ok = false;
      ++i;
    } while ( i != m_StripVector.end() && value_ok );
    if ( !value_ok )
    {
      log << MSG::INFO << " !!!! ILEGAL PDG-ID FOR STRIPOUT WAS REQUESTED => STRIPOUT ALL "
          << endmsg;
      m_StripVector = m_AllPartons;
    }
  }
  log << MSG::INFO << " THE FOLLOWING PARTON(S) WILL BE STRIPED OUT ";
  for ( std::vector<int>::const_iterator ip = m_StripVector.begin(); ip != m_StripVector.end();
        ++ip )
    log << *ip << "  ";
  log << endmsg;
}
 
void GenModule::StripPartons( GenEvent* evt ) {
  //    for ( HepMC::GenEvent::particle_iterator p = evt->particles_begin();
  //          p != evt->particles_end(); ++p )
  //    {
  //  //        std::cout << " PDG, BAR " << (*p)->pdg_id() << " " << (*p)->barcode() <<
  //  std::endl;        if ( std::find(m_StripVector.begin(),              m_StripVector.end(),
  //               (*p)->pdg_id()) != m_StripVector.end() )
  //        {
  //  //        std::cout << " REMOVING " << (*p)->pdg_id() << " " << (*p)->barcode()
  //  << std::endl;         HepMC::GenVertex*   pvtx    =   (*p)->production_vertex();
  //        HepMC::GenVertex*   evtx    =   (*p)->end_vertex();
  //        if (pvtx) pvtx->remove_particle(*p);
  //        if (evtx) evtx->remove_particle(*p);
  //        delete  *p;
  //        }
 
  //    }
 
  // Loop on all vertices of the event and remove the particle from the vertex
  for ( HepMC::GenEvent::vertex_iterator vtx = evt->vertices_begin();
        vtx != evt->vertices_end(); ++vtx )
  {
    // Loop on all children particles and remove the ones that should be stripped out
    for ( HepMC::GenVertex::particle_iterator p = ( *vtx )->particles_begin( HepMC::children );
          p != ( *vtx )->particles_end( HepMC::children ); ++p )
    {
      if ( std::find( m_StripVector.begin(), m_StripVector.end(), ( *p )->pdg_id() ) !=
           m_StripVector.end() )
      {
        if ( ( *p )->end_vertex() ) ( *p )->end_vertex()->remove_particle( *p );
        if ( ( *p )->production_vertex() ) ( *p )->production_vertex()->remove_particle( *p );
        delete *p;
      }
    }
    // Loop on all parents particles and remove the ones that should be stripped out
    for ( HepMC::GenVertex::particle_iterator p = ( *vtx )->particles_begin( HepMC::parents );
          p != ( *vtx )->particles_end( HepMC::parents ); ++p )
    {
      if ( std::find( m_StripVector.begin(), m_StripVector.end(), ( *p )->pdg_id() ) !=
           m_StripVector.end() )
      {
        if ( ( *p )->end_vertex() ) ( *p )->end_vertex()->remove_particle( *p );
        if ( ( *p )->production_vertex() ) ( *p )->production_vertex()->remove_particle( *p );
        delete *p;
      }
    }
  }
}