AthenaHepMCtoG4EventAction.cpp

Go to the documentation of this file.

//------------------------------------------------------------------
//
// ClassName:   AthenaHepMCtoG4EventAction
//
// Description: PrimaryGeneratorAction to convert HepMC event to G4Event
//
// Author:      Charles Leggett
//
// Date:        3-8-2001
//
// $Id: AthenaHepMCtoG4EventAction.cpp,v 1.1 2005/08/17 06:45:46 dengzy Exp $
// Simulation/G4Sim/G4Svc tag $Name: G4Svc-00-01-58 $
//
//------------------------------------------------------------------
 
#include "G4Svc/AthenaHepMCtoG4EventAction.h"
 
#include "HepMC/GenEvent.h"
 
#include "G4Event.hh"
#include "G4UImanager.hh"
#include "globals.hh"
 
#include "G4Geantino.hh"
#include "G4LorentzVector.hh"
#include "G4PrimaryParticle.hh"
#include "G4PrimaryVertex.hh"
 
#include "GaudiKernel/IMessageSvc.h"
 
using namespace HepMC;
 
AthenaHepMCtoG4EventAction::AthenaHepMCtoG4EventAction( int level ) : m_logLevel( level ) {}
 
AthenaHepMCtoG4EventAction::AthenaHepMCtoG4EventAction() : m_logLevel( int( MSG::INFO ) ) {}
 
AthenaHepMCtoG4EventAction::~AthenaHepMCtoG4EventAction() {}
 
//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 
// void AthenaHepMCtoG4EventAction::SetHepMCEvent(const HepMC::GenEvent *p) {
void AthenaHepMCtoG4EventAction::SetHepMCEvent( McEventCol::iterator p ) {
 
  if ( m_logLevel <= int( MSG::DEBUG ) )
  { std::cout << "AthenaHepMCtoG4EventAction: - Setting HepMC Event ptr" << std::endl; }
 
  //  p_evt = p;
  p_evtCollItr = p;
}
 
//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 
//  void AthenaHepMCtoG4EventAction::GeneratePrimaries(G4Event* anEvent)
//  {
//    if ( m_logLevel <= int(MSG::DEBUG) ) {
//      std::cout << "AthenaHepMCtoG4EventAction: - Generating Primaries"
//            << std::endl;
//    }
 
//    p_evt = (*p_evtCollItr)->pGenEvt();
 
//    if (p_evt == 0) {
 
//      if ( m_logLevel <= int(MSG::WARNING) ) {
//        std::cout << "AthenaHepMCtoG4EventAction: - HepMC ptr == null"
//          << std::endl;
//      }
 
//      return;
 
//    }
 
//    int nv = 1;
//    int mnv;
 
//    std::cout << "how many vertices: ";
//    std::cin >> mnv;
 
//    for ( HepMC::GenEvent::vertex_const_iterator v = p_evt->vertices_begin();
//          v != p_evt->vertices_end(); ++v ) {
//      bool IsGood=false;
//      for (HepMC::GenVertex::particle_iterator it=
//         (*v)->particles_begin(HepMC::children);
//       it!=(*v)->particles_end(HepMC::children);
//       it++) {
 
//        if (!(*it)->end_vertex() && (*it)->status()==1 ) {
//      IsGood=true;
//      break;
//        }
//      }
 
//      if (IsGood) {
//        std::cout << "found a good vertex, barcode: " << (*v)->barcode()
//          << std::endl;
 
//        G4LorentzVector lv=(*v)->position();
//        G4PrimaryVertex *vert= new G4PrimaryVertex(lv.x(),lv.y(),lv.z(),lv.t());
 
//        for (HepMC::GenVertex::particle_iterator it=
//           (*v)->particles_begin(HepMC::children);
//         it!=(*v)->particles_end(HepMC::children);
//         it++) {
 
//      //  if ( (*it)->status() != 1) continue;
 
//      int pdgcode=(*it)->pdg_id();
//      GenVertex *dcy = (*it)->end_vertex();
//      std::cout << "found a good particle, barcode: " << (*it)->barcode()
//            << "  pdg: " << pdgcode << "  status: " << (*it)->status()
//            << "  end: ";
//      if (dcy != 0) {
//        std::cout << dcy->barcode();
//      } else {
//        std::cout << "nil";
//      }
//      std::cout << std::endl;
 
//      const G4ThreeVector& p=(*it)->momentum();
//      G4PrimaryParticle *part;
//      if (pdgcode!=999) {
//        part=new G4PrimaryParticle (pdgcode,p.x()*GeV,p.y()*GeV,p.z()*GeV);
//      } else {
//        G4ParticleDefinition *pdef = G4Geantino::GeantinoDefinition();
//        part=new G4PrimaryParticle (pdef,p.x()*GeV,p.y()*GeV,p.z()*GeV);
//      }
//      vert->SetPrimary(part);
//        }
//        anEvent->AddPrimaryVertex(vert);
//        if (nv > mnv) { break; }
//        nv ++;
//      }
//    }
 
//    std::cout << "done generating event" << std::endl;
 
//    p_evtCollItr++;
 
//  }
 
//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 
void AthenaHepMCtoG4EventAction::GeneratePrimaries( G4Event* anEvent ) {
  if ( m_logLevel <= int( MSG::DEBUG ) )
  { std::cout << "AthenaHepMCtoG4EventAction: - Generating Primaries" << std::endl; }
 
  p_evt = (const HepMC::GenEvent*)( *p_evtCollItr );
 
  if ( p_evt == 0 )
  {
 
    if ( m_logLevel <= int( MSG::WARNING ) )
    { std::cout << "AthenaHepMCtoG4EventAction: - HepMC ptr == null" << std::endl; }
 
    return;
  }
 
  G4PrimaryVertex* vertex = new G4PrimaryVertex( 0., 0., 0., 0. );
 
  int n = 0;
  for ( HepMC::GenEvent::vertex_const_iterator v = p_evt->vertices_begin();
        v != p_evt->vertices_end(); ++v )
  {
    bool IsGood = false;
    for ( HepMC::GenVertex::particle_iterator it = ( *v )->particles_begin( HepMC::children );
          it != ( *v )->particles_end( HepMC::children ); it++ )
    {
 
      if ( !( *it )->end_vertex() && ( *it )->status() == 1 )
      {
        IsGood = true;
        break;
      }
    }
 
    if ( IsGood )
    {
      //        std::cout << "found a good vertex, barcode: " << (*v)->barcode()
      //        << std::endl;
 
      for ( HepMC::GenVertex::particle_iterator it =
                ( *v )->particles_begin( HepMC::children );
            it != ( *v )->particles_end( HepMC::children ); it++ )
      {
 
        if ( ( *it )->status() != 1 ) continue;
 
        int pdgcode = ( *it )->pdg_id();
        //      GenVertex *dcy = (*it)->end_vertex();
        //      std::cout << "found a good particle, barcode: " << (*it)->barcode()
        //            << "  pdg: " << pdgcode << "  status: " << (*it)->status()
        //            << "  end: ";
        //      if (dcy != 0) {
        //        std::cout << dcy->barcode();
        //      } else {
        //        std::cout << "nil";
        //      }
        //      std::cout << std::endl;
 
        const G4ThreeVector& p = ( *it )->momentum();
        G4PrimaryParticle*   part;
        if ( pdgcode != 999 )
        { part = new G4PrimaryParticle( pdgcode, p.x() * GeV, p.y() * GeV, p.z() * GeV ); }
        else
        {
          G4ParticleDefinition* pdef = G4Geantino::GeantinoDefinition();
          part = new G4PrimaryParticle( pdef, p.x() * GeV, p.y() * GeV, p.z() * GeV );
        }
 
        vertex->SetPosition( ( *v )->position().x(), ( *v )->position().y(),
                             ( *v )->position().z() );
        vertex->SetT0( ( *v )->position().t() );
 
        vertex->SetPrimary( part );
 
        n++;
      }
    }
  }
  anEvent->AddPrimaryVertex( vertex );
 
  if ( m_logLevel <= int( MSG::DEBUG ) )
  {
    std::cout << "AthenaHepMCtoG4EventAction: - done. G4Event has " << n << "  particles. "
              << std::endl;
  }
 
  p_evtCollItr++;
}