GtTimeratorTool.cc

Go to the documentation of this file.

#include "GtTimeratorTool.h"

#include "GaudiKernel/IssueSeverity.h"
#include "GaudiKernel/IRndmGenSvc.h"
#include "HepMC/GenEvent.h"
#include "HepMC/GenVertex.h"
#include "CLHEP/Units/SystemOfUnits.h"

#include <cmath>

GtTimeratorTool::GtTimeratorTool(const std::string& type,
                                 const std::string& name,
                                 const IInterface* parent)
    : GaudiTool(type,name,parent)
{
    declareInterface<IHepMCEventMutator>(this);

    declareProperty("Mode",m_mode="Absolut","Mode for time calculation method.");
    declareProperty("Distribution",m_distrib="Exponential","Event timing distribution");
    declareProperty("LifeTime",m_lifetime=0,"Life time for distribution.");
}

GtTimeratorTool::~GtTimeratorTool()
{
}

StatusCode GtTimeratorTool::initialize()
{
    info () << "initialize timer tool with mean delta time of " 
            << m_lifetime/CLHEP::second << endreq;

    // Non-positive lifetimes are meaningless and cause problems, particularly for 
    // the Sim15 method
    if ( m_lifetime <= 0 )
      {
        fatal() << "Non-positive lifetime specified. This is meaningless." << endreq;
        return StatusCode::FAILURE;
      }

    // Unfortunately GaudiTool doesn't provide a short cut for us
    IRndmGenSvc *rgs = 0;
    if (service("RndmGenSvc",rgs,true).isFailure()) {
        fatal() << "Failed to get random service" << endreq;
        return StatusCode::FAILURE;        
    }

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

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

StatusCode GtTimeratorTool::mutate(HepMC::GenEvent& event)
{
    // Exponential Distribution: t(n) = t(n-1) * exp(-t/tL)
    // random t(n) = -ln(U)*tL + t(n-1)

    double u = m_rand();
    double dt = 0;

    if (m_distrib == "Exponential")  dt = ((-1.0 * log(u)) * m_lifetime);
    else if (m_distrib == "Uniform") dt = m_lifetime;
    else {
        fatal() << "Not a recognized event timing distribution option" << endreq;
        return StatusCode::FAILURE;
    }
    
    debug() << "Incrementing current time of by " 
            << dt/CLHEP::second << " seconds"
            << endreq;
    
    double event_time = dt;

    HepMC::GenEvent::vertex_iterator vtx, done = event.vertices_end();
    for (vtx = event.vertices_begin(); vtx != done; ++vtx) {
        HepMC::FourVector position = (*vtx)->position();
        double vertex_time = event_time;
        if (m_mode == "Relative")
              vertex_time += position.t();
        position.setT(vertex_time);
        (*vtx)->set_position(position);
    }
    return StatusCode::SUCCESS;
}