In This Package:

HepMCtoG4.cc

Go to the documentation of this file.

#include "HepMCtoG4.h"
#include "G4DataHelpers/G4DhPrimaryVertexInformation.h"
#include "G4DataHelpers/G4DhPrimaryParticleInformation.h"

#include "HepMC/GenParticle.h"
#include "HepMC/GenEvent.h"
#include "HepMC/GenVertex.h"

#include "GaudiAlg/GaudiCommon.h"  //djaffe
#include "GaudiKernel/ParticleProperty.h" // particle mass


#include "G4PrimaryParticle.hh"
#include "G4PrimaryVertex.hh"
#include "G4OpticalPhoton.hh"

#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Geometry/Normal3D.h"



HepMCtoG4::HepMCtoG4(const std::string& type,
                     const std::string& name,
                     const IInterface* parent)
    : GaudiTool(type,name,parent)
      , m_pps(0) // djaffe
{
    declareInterface<IHepMCtoG4>(this);
    
    declareProperty("ZeroTime",m_zeroTime = true,
                    "Set true and all times will be reset so the signal vertex is at t=0.");
}

HepMCtoG4::~HepMCtoG4()
{
}


StatusCode HepMCtoG4::initialize()
{
  m_pps = svc<IParticlePropertySvc>( "ParticlePropertySvc", true ) ; //djaffe
  return StatusCode::SUCCESS;
}

StatusCode HepMCtoG4::finalize()
{
  return StatusCode::SUCCESS;
}
// djaffe: Set optical photon mass to be zero and get the mass of all other particles
//         from the particle property service. If the mass isn't available from the
//         service or the calculated mass doesn't agree with the mass from the service
//         then issue a warning.
//         No warning is issued for nuclei that have been created by GenDecay with zero energy,
//         set Trac ticket #364 for details.
G4PrimaryParticle*
HepMCtoG4::makeG4particle(const HepMC::GenParticle& part, const HepMC::GenEvent* event, const HepMC::GenVertex* vertex)
{
    HepMC::FourVector p = part.momentum();

    G4PrimaryParticle* g4part = 0;
    
    double m, cm ;
    if (part.pdg_id() == 20022 /*opticalphoton*/) {
        g4part = new G4PrimaryParticle(G4OpticalPhoton::Definition(),
                                       p.px(),p.py(),p.pz());
        m = 0. ;
    }
    else {
        g4part = new G4PrimaryParticle(part.pdg_id(),p.px(),p.py(),p.pz());
        m = cm = sqrt(fabs( p.t()*p.t() - (p.x()*p.x() + p.y()*p.y() + p.z()*p.z()) ));

        ParticleProperty*  ppp = m_pps->findByStdHepID( part.pdg_id() );
        if ( ppp  ) 
          {
            m = ppp->mass();
            double r = 0.;
            if ( m > 0 ) { r = std::abs(cm-m)/m; } 
            if ( r > 0.001 || std::abs(cm-m) > 0.01*MeV )
              {
                warning() << "calculated mass = " << cm
                          << " differs from IParticlePropertySvc mass " << m
                          << " for pdg_id " << part.pdg_id()
                          << ". Using IParticlePropertySvc mass. " << endreq;
              }
          }
        else
          {
            if ( p.t() < 0.01*MeV && std::abs(part.pdg_id()) > 1000000000 )
              {
                debug() << "As expected, there is no IParticlePropertySvc mass for pdg_id "
                        << part.pdg_id() << " which is a zero-energy, recoil nucleus created by GenDecay. Use calculated mass " << m 
                        << " from particle energy " << p.t() << " momentum " << sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z())
                        << endreq;
              }
            else
              {
                warning() << "cannot find IParticlePropertySvc mass for pdg_id "
                          << part.pdg_id() << " and status " << part.status() << " . Using calculated mass " << m 
                          << " from particle energy " << p.t() << " momentum " << sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z())
                          << endreq;
              }
          }

    }
    g4part->SetMass(m);
    HepMC::ThreeVector pol = part.polarization().normal3d();
    g4part->SetPolarization(pol.x(),pol.y(),pol.z());

    g4part->SetUserInformation(new G4DhPrimaryParticleInformation(event,vertex,&part));
    return g4part;        
}

StatusCode HepMCtoG4::convert(const HepMC::GenEvent& input,
                              std::vector<G4PrimaryVertex*>& output)
{

    // Catch HepMC::GenEvent weighting
    bool hasEventWeight = false;
    double eventWeight = 1.0;
    // If HepMC event weights container has only one value, apply to
    // all vertices as described in HepMC v2 documentation.
    if( input.weights().size() == 1 ){ 
      hasEventWeight = true; 
      eventWeight *= input.weights().front();
    }

    // Loop over vertices in the event
    HepMC::GenEvent::vertex_const_iterator
        iVtx = input.vertices_begin(),
        doneVtx = input.vertices_end();
    for (/*nop*/; doneVtx != iVtx; ++iVtx) {
            
        G4PrimaryVertex* g4vtx = new G4PrimaryVertex;
        bool first_time = true;
        
        // Catch HepMC::GenVertex weighting
        bool hasVertexWeight = false;
        double vertexWeight = 1.0;
        // If HepMC::GenVertex weights container has only one value,
        // apply to current vertex as described in HepMC v2
        // documentation.
        if( (*iVtx)->weights().size() == 1 ){ 
          hasVertexWeight = true;
          vertexWeight *= (*iVtx)->weights().front();
        }
        if( hasEventWeight || hasVertexWeight ){
          g4vtx->SetWeight( eventWeight*vertexWeight );
        }

        // Loop over particles in the vertex
        HepMC::GenVertex::particles_out_const_iterator 
            iPart = (*iVtx)->particles_out_const_begin(),
            donePart = (*iVtx)->particles_out_const_end();
        for (/*nop*/; donePart != iPart; ++iPart) {
                
            // Only keep particles that are important for tracking
            if ((*iPart)->status() != 1) continue;

            if (first_time) {
                first_time = false;
                HepMC::FourVector v = (*iPart)->production_vertex()->position();
                g4vtx->SetPosition(v.x(),v.y(),v.z());
                if (m_zeroTime) 
                    g4vtx->SetT0(0.0);
                else
                    g4vtx->SetT0(v.t());
                g4vtx->SetUserInformation(new G4DhPrimaryVertexInformation(&input,*iVtx));
            }
            g4vtx->SetPrimary(makeG4particle(*(*iPart),&input,*iVtx));
        }

        if (!g4vtx->GetNumberOfParticle()) {
            delete g4vtx;
            continue;
        }

        output.push_back(g4vtx);
    }
    return StatusCode::SUCCESS;
}

---