GtGenerator.cc

Go to the documentation of this file.

#include "GtGenerator.h"
#include "GenTools/IHepMCEventMutator.h"

#include "Event/GenHeader.h"

#include "Context/TimeStamp.h"

#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/IDataManagerSvc.h"
#include "GaudiKernel/IRndmEngine.h"

#include "GaudiKernel/IssueSeverity.h"

#include "CLHEP/Units/SystemOfUnits.h"
#include "HepMC/GenEvent.h"

#include <vector>
#include <limits>

GtGenerator::GtGenerator(const std::string& name, ISvcLocator* pSvcLocator)
    : DybAlgorithm<DayaBay::GenHeader>(name,pSvcLocator)
{
    declareProperty("TimeStamp",m_tsseconds = 0,
                    "Absolute start time in seconds from the unix epoch. "
                    "Caution: this is a bare integer, do not multiply by CLHEP seconds.");
    declareProperty("GenTools",m_genToolNames,
                    "Tools to generate HepMC::GenEvents");
    declareProperty("GenName", m_genName = "GtGenerator",
                    "Name of this generator for book keeping purposes.");
}

GtGenerator::~GtGenerator()
{
}

StatusCode GtGenerator::initialize()
{
    for (size_t ind=0; ind < m_genToolNames.size(); ++ind) {
        std::string gtn = m_genToolNames[ind];
        try {
            m_genTools.push_back(tool<IHepMCEventMutator>(gtn));
        }
        catch(const GaudiException& exg) {
            fatal() << "Failed to get generator: \"" << gtn << "\"" << endreq;
            return StatusCode::FAILURE;
        }
        info () << "Added gen tool " << gtn << endreq;
    }

    // Set the generator start time
    m_now = TimeStamp(m_tsseconds);

    return StatusCode::SUCCESS;
}

StatusCode GtGenerator::execute()
{
    // Look for pre-existing header object or make new one
    DayaBay::GenHeader* gen_header = MakeHeaderObject();


    HepMC::GenEvent* event = new HepMC::GenEvent;
    event->set_event_number(this->getExecNum());

    // Let each tool do its thing on the event
    for (size_t ind = 0; ind< m_genTools.size(); ++ind) {
        debug () << "Running gen tool #" << ind << " " << m_genToolNames[ind] << endreq;
        if (m_genTools[ind]->mutate(*event).isFailure()) {
            fatal() << "Generator Tool " << m_genToolNames[ind]
                    << " failed" << endreq;
            return StatusCode::FAILURE;
        }
    }

    // Find the smallest time
    double dt = std::numeric_limits<double>::max();
    {
        HepMC::GenEvent::vertex_const_iterator it, done = event->vertices_end();
        for (it = event->vertices_begin(); it != done; ++it) {
            double t = (*it)->position().t();
            if (t < dt) dt = t;
        }
    }
    debug () << "Smallest time of "
             << event->vertices_size() << " vertices = " << dt << endreq;

    // Slurp away seconds and ns and leave the rest
    dt = dt/CLHEP::second;      // bring into explicit units of seconds
    time_t seconds = (time_t)dt;
    double fractional = dt - seconds;
    int ns = (int)(1e9*fractional);
    double rest = fractional - (1e-9*ns);

    // how much we steal away from each vertex time
    double toffset = (dt - rest)*CLHEP::second; // put back into system of units

    // Increment our notion of now
    m_now.Add(TimeStamp(seconds,ns));
    gen_header->setTimeStamp(m_now);

    debug() << "Offset time is " << toffset
            << ", now is " << m_now << endreq;

    // Loop over vertices to set reduce time we stole and set time extent
    {
        double earliest = -1, latest = -1;
        HepMC::GenEvent::vertex_iterator it, done = event->vertices_end();
        for (it = event->vertices_begin(); it != done; ++it) {
            HepMC::FourVector vtx = (*it)->position();
            double t = vtx.t() - toffset;
            vtx.setT(t);
            (*it)->set_position(vtx);

            if (earliest < 0 || earliest > t) earliest = t;
            if (latest   < 0 || latest   < t)   latest = t;
        }

        TimeStamp ts_earliest(m_now);
        if (earliest > 0)
            ts_earliest.Add(earliest/CLHEP::second);
        gen_header->setEarliest(ts_earliest);

        TimeStamp ts_latest(m_now);
        if (latest > 0)
            ts_latest.Add(latest/CLHEP::second);
        gen_header->setLatest(ts_latest);
    }

    // GenHeader subclass
    gen_header->setEvent(event);
    gen_header->setGeneratorName(m_genName);

    return StatusCode::SUCCESS;
}

StatusCode GtGenerator::finalize()
{
    // Do NOT release the tools.  This is handled implicitly
    m_genTools.clear();

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