Gnrtr.cc

Go to the documentation of this file.

#include "Gnrtr.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 "HepMC/GenEvent.h"
 
#include <vector>
#include "CLHEP/Units/SystemOfUnits.h"
#include <limits>

Gnrtr::Gnrtr(const std::string& name, ISvcLocator* pSvcLocator)
    : StageProcessor<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("TimeStampNano",m_tsnanoseconds = 0,
                    "Time offset in nanoseconds to add to the absolute start time");

    m_pullMode=true;
    
    GT_constructor();
}

Gnrtr::~Gnrtr()
{
    GT_destructor();
}

StatusCode Gnrtr::initialize()
{
    StatusCode sc = this->StageProcessor<DayaBay::GenHeader>::initialize();
    if (sc.isFailure()) return sc;

    m_Start=true;

    // Although initialized to 0, however, no special meaning
    m_CurrentTime=0;

    m_now = TimeStamp(m_tsseconds,m_tsnanoseconds);

    return GT_initialize();
}

StatusCode Gnrtr::execute()
{
    if(m_Start || m_CurrentTime <= thisStage()->currentTime()) {
        m_Start=false;
        GnrtrData* pGD=0;
        //
        // Because this execute() is not triggered by sysExecute()
        // preExecute() and postExecute() are called here explicitly.
        preExecute();
        StatusCode sc = GT_execute(pGD);
        if(sc.isFailure()) {
            return sc;
        }
        postExecute();
        //
        thisStage()->pushElement(pGD);
        this->registerData(*pGD);
        m_CurrentTime=pGD->time();
        debug()<<"to grep: new data pushed out at time "<<pGD->time()<<endreq;
    }

    return StatusCode::SUCCESS;
}

StatusCode Gnrtr::finalize()
{
    return GT_finalize();
}


void Gnrtr::GT_constructor()
{
    declareProperty("GenTools",m_genToolNames,
                    "Tools to generate HepMC::GenEvents");
    declareProperty("GenName", m_genName = "GtGenerator",
                    "Name of this generator for book keeping purposes.");
}


void Gnrtr::GT_destructor()
{
}

StatusCode Gnrtr::GT_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;
    }

    return StatusCode::SUCCESS;
}

StatusCode Gnrtr::GT_execute(GnrtrData*& p_output)
{
    // 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();
            debug()<<"time of this vertex: "<<t<<endreq;
            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);

    // >> Fifteen
    p_output=new GnrtrData(*gen_header);
    // << Fifteen

    return StatusCode::SUCCESS;

}

StatusCode Gnrtr::GT_finalize()
{
    for (size_t ind=0; ind < m_genTools.size(); ++ind)
        m_genTools[ind]->release();
    m_genTools.clear();

    return this->StageProcessor<DayaBay::GenHeader>::finalize();
}

// Local Variables:
// c-basic-offset: 4
// End: