ManualGenerator.cxx

Go to the documentation of this file.

// *****************************************************************************
//
// ManualGenerator.cxx
//
//  A package to read four-vectors from an input file and send them through
//    the BES detector simulation. 
//
//  Notes:
//
//    *  To use, replace EvtGen in the job options file with:
//
//           #include "$MANUALGENERATORROOT/share/ManualGenerator.txt"
//           ManualGenerator.InputFile = "inputfile.txt";
//
//       where "inputfile.txt" is a text file containing four-vectors in
//       the following format for each event:
//
//           [number of final state particles in this event]
//           [PDG ID of particle 1]  [P_x] [P_y] [P_z] [E]
//           [PDG ID of particle 2]  [P_x] [P_y] [P_z] [E]
//           and so on.
//
//    *  This package deals only with final state particles.  The production
//       of the initial state particle (for example, the J/psi) must be handled
//       by a prior package (for example, KKMC).
//
//    *  Events should be generated with zero net 3-momentum.  Particles are
//       then internally boosted to account for the crossing angle, beam energy
//       spread, etc., as reflected by the four-momentum of the initial state
//       particle.
//
//    *  Events should be generated at the origin.  Particles are then offset
//       according to the production vertex of the initial state particle.
//
//    *  pi0 and Ks decays are handled by GEANT.  
//         (Note that B(pi0 --> gamma gamma) is not 100%.)
//
//
// September 2011 R. Mitchell
//
// --2012-8-10, zhao qingzhang and ping ronggang
//  Boost from CM system to Lab system has been implement out of generator, so 
//  remove the boost 
// *****************************************************************************
 
#include "ManualGenerator/ManualGenerator.h"
 
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/DataSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include "GeneratorObject/McGenEvent.h"
 
#include <stdlib.h>
 
 
// *****************************************************************************
//
//    CONSTRUCTOR
//
// *****************************************************************************
 
ManualGenerator::ManualGenerator(const string& name, ISvcLocator* pSvcLocator):Algorithm( name, pSvcLocator ){
 
 
    // get the name of the input file from the job options file
 
  declareProperty("InputFile", m_inputFileName = "");
 
}
 
 
 
// *****************************************************************************
//
//    INITIALIZE
//
// *****************************************************************************
 
 
StatusCode ManualGenerator::initialize(){
  //MsgStream log(messageService(), name());
  //log << MSG::INFO << "ManualGenerator initialize" << endreq;
 
 
    // open the input file
 
  m_inputFile.open(m_inputFileName.c_str());
  if (!m_inputFile){
    cout << "ManualGenerator: PROBLEMS OPENING FILE "
         << m_inputFileName << endl;
    exit(0);
  }
 
 
    // zero the event counter
 
  m_eventCounter = 0;
 
 
  return StatusCode::SUCCESS;
 
}
 
 
 
// *****************************************************************************
//
//    EXECUTE
//
// *****************************************************************************
 
 
StatusCode ManualGenerator::execute(){ 
 
 
    // read an event from the input file
 
  int nParticles;
  int idParticles[100];
  double pxParticle;
  double pyParticle;
  double pzParticle;
  double eParticle;
  CLHEP::HepLorentzVector pParticles[100];
 
  if (m_inputFile >> nParticles){
    for (int i = 0; i < nParticles; i++){
      m_inputFile >> idParticles[i];
      m_inputFile >> pxParticle;  pParticles[i].setPx(pxParticle);
      m_inputFile >> pyParticle;  pParticles[i].setPy(pyParticle);
      m_inputFile >> pzParticle;  pParticles[i].setPz(pzParticle);
      m_inputFile >> eParticle;   pParticles[i].setE(eParticle);
    }
  }
 
  else{
    cout << "ManualGenerator: RAN OUT OF INPUT EVENTS FROM FILE " 
         << m_inputFileName << endl;
    return StatusCode::SUCCESS;
  }
  
 
 
    // print out a message for this event to make sure all is well
 
  cout << "ManualGenerator: PROCESSING EVENT " << ++m_eventCounter << endl;
  cout << nParticles << endl;
  for (int i = 0; i < nParticles; i++){
    cout << idParticles[i] << "\t";
    cout << pParticles[i].px() << "\t";
    cout << pParticles[i].py() << "\t";
    cout << pParticles[i].pz() << "\t";
    cout << pParticles[i].e() << endl;
  }
 
 
    // retrieve the original MC event (HepMC::GenEvent) from the framework
 
  HepMC::GenEvent* genEvent = 0;
  SmartDataPtr<McGenEventCol> initialMcCol(eventSvc(), "/Event/Gen");
  if (initialMcCol == 0){
//    HepMC::GenEvent* genEvent = new GenEvent();  
//    genEvent->set_event_number(0);
 
    //create a Transient store 
//    McGenEventCol *mcColl = new McGenEventCol;
//    McGenEvent* mcEvent = new McGenEvent(genEvent);
//    mcColl->push_back(mcEvent);
//    StatusCode sc = eventSvc()->registerObject("/Event/Gen", mcColl);
//    if (sc != StatusCode::SUCCESS) return StatusCode::FAILURE;    
 
    cout << "ManualGenerator: COULD NOT FIND THE ORIGINAL MC EVENT COLLECTION" << endl;
    exit(0);
  }
//  else
//  {
    genEvent = (*(initialMcCol->begin()))->getGenEvt();
//  }
  //genEvent->print();
 
 
 
    // loop over MC particles and look for the particle we want to decay
 
  HepMC::GenParticle* decayParticle = NULL;
  for(HepMC::GenEvent::particle_const_iterator igenParticle = genEvent->particles_begin();
      igenParticle != genEvent->particles_end(); ++igenParticle ){
    decayParticle = *igenParticle;
    if (decayParticle->status() == 1 && abs(decayParticle->pdg_id()) != 22) break;
  }
  if (!decayParticle){
    cout << "ManualGenerator: COULD NOT FIND A PARTICLE TO DECAY" << endl;
    exit(0);
  }
 
 
    // create a new vertex from the decaying particle's production vertex
    //  and add it to the event
 
  HepMC::GenVertex* initialVertex = decayParticle->production_vertex();
  HepMC::GenVertex* newVertex = new HepMC::GenVertex(initialVertex->position());
  genEvent->add_vertex(newVertex);
  newVertex->add_particle_in(decayParticle);
 
 
 
    // boost the input four-vectors into the decaying particle's rest frame
 
  CLHEP::HepLorentzVector initialFourMomentum(decayParticle->momentum().px(),
                                              decayParticle->momentum().py(),
                                              decayParticle->momentum().pz(),
                                              decayParticle->momentum().e());
  CLHEP::HepLorentzVector finalFourMomentum(0,0,0,0);
  for (int i = 0; i < nParticles; i++){
  //  pParticles[i].boost(initialFourMomentum.boostVector());
    finalFourMomentum += pParticles[i];
  }
  if ((fabs(initialFourMomentum.px() - finalFourMomentum.px()) > 0.0001) || 
      (fabs(initialFourMomentum.py() - finalFourMomentum.py()) > 0.0001) ||
      (fabs(initialFourMomentum.pz() - finalFourMomentum.pz()) > 0.0001) ||
      (fabs(initialFourMomentum.e()  - finalFourMomentum.e())  > 0.01)){
    cout << "ManualGenerator: INITIAL AND FINAL STATES HAVE INCONSISTENT 4-MOMENTA" << endl;
    cout << "delta px = " << initialFourMomentum.px() - finalFourMomentum.px() << endl; 
    cout << "delta py = " << initialFourMomentum.py() - finalFourMomentum.py() << endl; 
    cout << "delta pz = " << initialFourMomentum.pz() - finalFourMomentum.pz() << endl; 
    cout << "delta e = " << initialFourMomentum.e() - finalFourMomentum.e() << endl; 
   // exit(0);
  }
 
 
    // add the input four-vectors to the new vertex
 
  for (int i = 0; i < nParticles; i++){
    newVertex->add_particle_out(new GenParticle(pParticles[i], idParticles[i], 1));
  }
  //genEvent->print();
 
 
 
  return StatusCode::SUCCESS; 
}
 
 
 
 
// *****************************************************************************
//
//    FINALIZE
//
// *****************************************************************************
 
 
StatusCode ManualGenerator::finalize() 
{
  m_inputFile.close();
  return StatusCode::SUCCESS;
}