In This Package:

GtBeamerTool.cc

Go to the documentation of this file.

#include "GtBeamerTool.h"

#include "DetDesc/IDetectorElement.h"
#include "DetDesc/IGeometryInfo.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"

using namespace CLHEP;

GtBeamerTool::GtBeamerTool(const std::string& type,
                           const std::string& name,
                           const IInterface* parent)
    : GaudiTool(type,name,parent)
    , m_radius(0)
    , m_particlesPerEvent(1)
    , m_particleName("opticalphoton")
    , m_targetOffset(3,0)
    , m_sourceDirection(3,0)
    , m_distance(CLHEP::meter)
    , m_position(3,0)
    , m_direction(3,0)
    , m_pid(0)
{
    declareInterface<IHepMCEventMutator>(this);

    declareProperty("Radius",m_radius=CLHEP::meter,
                    "Radius of source disk.");

    declareProperty("ParticlesPerEvent",m_particlesPerEvent=1,
                    "Number of particles per event");

    declareProperty("ParticleName",m_particleName="opticalphoton",
                    "Name of particle.");

    declareProperty("Momentum",m_momentum=2.5*CLHEP::eV,
                    "Momentum of particles in beam.");

    declareProperty("TargetElement",m_targetElement="",
                    "DetecorElement to aim at.");
    
    declareProperty("TargetOffset",m_targetOffset,
                    "Offset from target origin to aim for.");

    declareProperty("SourceDirection",m_sourceDirection,
                    "Direction from target to source.");

    declareProperty("SourceDistance",m_distance,
                    "Distance from target to source.");

}

GtBeamerTool::~GtBeamerTool()
{
}

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

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

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

    // First try to find via PID
    double mass = 0;
    if (m_pid) {
        particle = ppSvc->findByStdHepID(m_pid);
        if (particle) {
            mass = particle->mass();
            m_particleName = particle->particle();
        }
    }
    // If failed, then try via name
    if (!particle) {
        particle = ppSvc->find(m_particleName);
        if (particle) {
            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;
    }


    if ("" == m_targetElement) {
        error() << "TargetElement not given" << endreq;
        return StatusCode::FAILURE;
    }
    IDetectorElement* targetDE =
        getDet<IDetectorElement>(m_targetElement);


    // Convert local direction to global coordinates
    Gaudi::XYZVector local_direction(m_sourceDirection[0],
                                     m_sourceDirection[1],
                                     m_sourceDirection[2]);
    Gaudi::XYZVector global_direction = 
        targetDE->geometry()->toGlobal(local_direction);
    m_direction.setX(global_direction.x());
    m_direction.setY(global_direction.y());
    m_direction.setZ(global_direction.z());
    m_direction = m_direction.unit();


    // Convert target point to global coordinates
    Gaudi::XYZPoint local_target(m_targetOffset[0],
                                 m_targetOffset[1],
                                 m_targetOffset[2]);
    Gaudi::XYZPoint global_target = 
        targetDE->geometry()->toGlobal(local_target);

    m_position.setX(global_target.x() - m_distance*global_direction.x());
    m_position.setY(global_target.y() - m_distance*global_direction.y());
    m_position.setZ(global_target.z() - m_distance*global_direction.z());


    m_fourmom = HepLorentzVector(m_momentum*m_direction,
                                 sqrt(m_momentum*m_momentum + mass*mass));


    return StatusCode::SUCCESS;
}

StatusCode GtBeamerTool::finalize()
{
    return StatusCode::SUCCESS;
}

StatusCode GtBeamerTool::mutate(HepMC::GenEvent& event)
{
    for (int ind=0; ind<m_particlesPerEvent; ++ind) {

        Hep3Vector ortho = m_direction.orthogonal();
        
        double angle = m_uni()*360*CLHEP::degree;
        ortho = ortho.rotate(angle,m_direction).unit();
        
        double ran_radius = sqrt(m_uni()*m_radius*m_radius);
        ortho = ran_radius*ortho + m_position;

        HepLorentzVector position(ortho,0);
        HepMC::GenVertex* vertex = new HepMC::GenVertex(position);

        HepMC::GenParticle* part = 
            new HepMC::GenParticle(m_fourmom,m_pid,1/*=status*/);
        double phi = m_uni()*360*CLHEP::degree;
        part->set_polarization(HepMC::Polarization(0.5*HepMC::HepMC_pi,phi));
        vertex->add_particle_out(part);

        if (!event.vertices_size()) 
            event.set_signal_process_vertex(vertex);
        event.add_vertex(vertex);
    }

    return StatusCode::SUCCESS;
}

---