GtGunGenTool.cc

Go to the documentation of this file.

#include "GtGunGenTool.h"

#include "GaudiKernel/IRndmGenSvc.h"
#include "GaudiKernel/IParticlePropertySvc.h"
#include "GaudiKernel/ParticleProperty.h"
#include "GaudiKernel/Vector3DTypes.h"

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

#include "CLHEP/Units/SystemOfUnits.h"
#include "CLHEP/Vector/ThreeVector.h"

using namespace CLHEP;

GtGunGenTool::GtGunGenTool(const std::string& type,
                           const std::string& name, 
                           const IInterface* parent)
    : GaudiTool(type,name,parent)
    , m_direction(3,0)
    , m_pid(0)
    , m_mass(0)
{
    declareInterface<IHepMCEventMutator>(this);

    declareProperty("ParticleName", m_particleName = "", 
                    "PDG particle name");
    declareProperty("ParticlesPerEvent",m_particlesPerEvent = 1,
                    "Number of particles per event to generate.");
    declareProperty("MomentumMode", m_momentumMode = "Fixed","Momentum mode");
    declareProperty("Momentum", m_momentum = 1.0*MeV,"Particle momentum");
    m_momentum.verifier().setLower(0.*eV);
    declareProperty("MomentumSpread", m_momentumSpread = 0,"Momentum spread");
    declareProperty("MomentumInterpretation",m_momentumInterpretation = "Momentum","Valid interpretation values are Momentum(default), KineticEnergy, TotalEnergy") ;
    declareProperty("DirectionMode", m_directionMode = "Fixed","Direction mode");
    m_direction[2] = 1.0;
    declareProperty("Direction", m_direction,"Direction");
    declareProperty("DirectionSpread",  m_directionSpread = 0,"Direction spread for cos(zenith angle)");
    declareProperty("PolarizeMode", m_polarizeMode = "None",
                    "Type of polarization of the particle.");
}

GtGunGenTool::~GtGunGenTool()
{
}

StatusCode GtGunGenTool::initialize()
{
    this->GaudiTool::initialize();

    // Get mass of particle for later calculation of energy
    IParticlePropertySvc * ppSvc = 
        svc< IParticlePropertySvc >( "ParticlePropertySvc" , true ) ;
    ParticleProperty* particle = 0;

    // First try to find via PID
    if (m_pid) {
        particle = ppSvc->findByStdHepID(m_pid);
        if (particle) {
            m_mass = particle->mass();
            m_particleName = particle->particle();
        }
    }
    // If failed, then try via name
    if (!particle) {
        particle = ppSvc->find(m_particleName);
        if (particle) {
            m_mass = particle->mass();
            m_pid = particle->pdgID();
        }
    }
    // If still fail, complain
    if (!particle) {
        fatal() << "Failed to find particle named \"" 
                << m_particleName << "\" and with ID " 
                << m_pid << endreq;
        return StatusCode::FAILURE;
    }

    IRndmGenSvc *rgs = 0;
    if (service("RndmGenSvc",rgs,true).isFailure()) {
        fatal() << "Failed to get random service" << endreq;
        return StatusCode::FAILURE;        
    }

    StatusCode sc;

    sc = m_gauss.initialize(rgs, Rndm::Gauss(0,1));
    if (sc.isFailure()) {
        fatal() << "Failed to initialize gaussian random numbers" << endreq;
        return StatusCode::FAILURE;
    }
    sc = m_uni.initialize(rgs, Rndm::Flat(0,1));
    if (sc.isFailure()) {
        fatal() << "Failed to initialize uniform random numbers" << endreq;
        return StatusCode::FAILURE;
    }

    // force normalization
    Hep3Vector unit(m_direction[0],m_direction[1],m_direction[2]);
    unit = unit.unit();
    for (int ind=0; ind<3; ++ind) m_direction[ind] = unit[ind];
 
    return StatusCode::SUCCESS;
}

StatusCode GtGunGenTool::finalize()
{

    return this->GaudiTool::finalize();
}


StatusCode GtGunGenTool::mutate(HepMC::GenEvent& event)
{
    for (int ind=0; ind<m_particlesPerEvent; ++ind) {
        if (this->oneVertex(event).isFailure()) 
            return StatusCode::FAILURE;
    }
    return StatusCode::SUCCESS;
}
StatusCode GtGunGenTool::oneVertex(HepMC::GenEvent& event)
{
    CLHEP::Hep3Vector dir;

    if (m_directionMode == "Fixed") {
        dir.set(m_direction[0],m_direction[1],m_direction[2]);
    }
    else {
        // cos(angle) from mean direction
        double costh=0;
        do {
            costh = generateNumber(m_directionMode,1,m_directionSpread);
        } while (costh < -1.0||costh > 1.0);
        
        double sinth = sqrt(1-costh*costh);

        // azimuth angle around mean direction
        double phi = generateNumber("Uniform",0.0,360*CLHEP::degree);
        double cosphi = cos(phi);
        double sinphi = sin(phi);
        
        dir.set(sinth*cosphi,sinth*sinphi,costh);

        // rotate along Y-axis by the mean direction's theta
        dir.rotateY(acos(m_direction[2]));
        // rotate along Z-axis by the mean direction's phi
        dir.rotateZ(atan2(m_direction[1],m_direction[0]));

    }

    dir = dir.unit();

    double momentum;
    do{
      momentum = generateNumber(m_momentumMode,m_momentum,m_momentumSpread);
    } while (momentum < 0.0);

    // What is the interpretation of momentum?
    double energy = sqrt(momentum*momentum + m_mass*m_mass); 
    if (m_momentumInterpretation == "TotalEnergy") {
      energy = momentum ; 
      if ( energy < m_mass ) {
        fatal() << "Cannot set total energy(MeV) = " << energy/CLHEP::MeV 
                << " of " << m_particleName 
                << " to be less than particle mass(MeV) = " << m_mass/CLHEP::MeV << endreq;
        return StatusCode::FAILURE;
      }
      momentum = sqrt( (energy+m_mass) * (energy-m_mass) );
    }

    if (m_momentumInterpretation == "KineticEnergy") {
      energy = momentum + m_mass;
      momentum = sqrt( (energy+m_mass) * (energy-m_mass) );
    }
    
    HepMC::ThreeVector p(momentum*dir);
    HepMC::FourVector fourmom(p.x(),p.y(),p.z(),energy);

    debug() << m_particleName << " (" << m_pid << ") generated with 4momentum = " 
            << momentum
            << " dir = " << dir
            << " momentum(MeV) = " << momentum/CLHEP::MeV
            << " energy(MeV) = " << energy/CLHEP::MeV
            << " kinetic energy(MeV) = " << (energy-m_mass)/CLHEP::MeV
            << endreq;

    HepMC::GenParticle* particle = new HepMC::GenParticle(fourmom,m_pid,1/*=status*/);

    if (m_polarizeMode == "Random") {
        double phi = generateNumber("Uniform",0.0,360*CLHEP::degree);
        particle->set_polarization(HepMC::Polarization(0.5*HepMC::HepMC_pi,phi));
    }

    HepMC::GenVertex* vertex = event.signal_process_vertex();
    if (!vertex) {
        vertex = new HepMC::GenVertex(HepMC::FourVector(0,0,0,0));
        event.set_signal_process_vertex(vertex);
    }

    vertex->add_particle_out(particle);

    return StatusCode::SUCCESS;
}

double GtGunGenTool::generateNumber(const std::string& mode, double mean, double spread)
{
    if (mode == "Smeared") {
        return m_gauss() * spread + mean;
    }
    if (mode == "Uniform") {
        return (m_uni()*2-1) * spread +mean ;
    }
    return mean;
}