/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Generator/GeneratorModule/GeneratorModule-00-01-05/src/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/MsgStream.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/DataSvc.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/Ranlux64Engine.h"
#include "CLHEP/Random/RandPoisson.h"

//#include "BesHepMC/GenVertex.h"
#include "HepMC/GenVertex.h"

#include "GeneratorObject/McGenEvent.h"


//---------------------------------------------------------------------------
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(messageService(), 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" << endreq;
    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 "
          << endreq;
    }
    else {
      log << MSG::INFO << "Fixed Number of Interactions per Event is: "
          << m_meanInteractions << endreq;
    }
  }
  else {
    m_pPoissonGenerator = new CLHEP::RandPoisson(*m_pRandomEngine, 
                                          m_meanInteractions);

    log << MSG::INFO << 
      "Poisson Distribution of Interactions per Event with Mean: "
        << m_meanInteractions << endreq;
  }
  // Initialize the generator itself
  StatusCode status = this->genInitialize();
  if (status.isFailure()) {
    log << MSG::ERROR << "Could not initialize Generator properly" << endreq;
    return status;
  }  
  StatusCode  status1 = this->genuserInitialize();
  if (status1.isFailure()) {
    log << MSG::ERROR << "Could not initialize user part properly" << endreq;
    return status1;
  } 
  return status;
}

//---------------------------------------------------------------------------
StatusCode GenModule::execute() {
  //---------------------------------------------------------------------------
    
  int numToGenerate;
  StatusCode status;

  MsgStream log(messageService(), name());

  log << MSG::DEBUG << "GenModule::execute()" << endreq;

  // 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(messageService(), name());
      log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
      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(messageService(), name());
        log << MSG::ERROR << "Could not register McGenEvent" << endreq;
        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(messageService(), 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 "
        << endreq;
    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 "
          << endreq;
      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 << endreq;
  
}

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;
      }
    }
  }
}

---