In This Package:

GenHists.cc

Go to the documentation of this file.

#include "GenHists.h"

#include "Event/GenHeader.h"

#include "DetDesc/IGeometryInfo.h"
#include "DetDesc/DetectorElement.h"
#include "DetDesc/SolidBase.h"

#include "GaudiKernel/ITHistSvc.h"

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

#include "HepMC/GenEvent.h"

#include "TH1F.h"
#include "TH2F.h"
#include "TH1I.h"

#include <string>
using namespace std;


GenHists::GenHists(const std::string& name, ISvcLocator* pSvcLocator)
    : GaudiAlgorithm(name,pSvcLocator)
    , m_hsvc(0)
    , m_nKineParts(0)
    , m_kineEnergy(0)
    , m_xy(0), m_yz(0), m_xz(0)
{
    declareProperty("Location",m_location=DayaBay::GenHeaderLocation::Default,
                    "Location in the TES to get GenHeader");
    declareProperty("Volume",m_volume="",
                    "TDS path of volume used to generated vertices");
    declareProperty("FilePath",m_filepath="/file1/gen/",
                    "File path with with to register histograms.");
    declareProperty("MaxEnergy",m_maxEnergy = 10*CLHEP::MeV,
                    "Maximum energy for plot bounds.");
    declareProperty("EnergyUnits",m_energyUnits = CLHEP::MeV,
                    "Energy units for plot bounds.");
}

GenHists::~GenHists()
{
}

// These just save some typing below.  
#define REG1D(type,p,name,title,n,min,max)                              \
    do {                                                                \
        p = new type(name,title,n,min,max);                             \
        string path = m_filepath + name;                                \
        if (m_hsvc->regHist(path.c_str(),p).isFailure()) {              \
            error() << "Could not register \"" << path << endreq;       \
            delete p; p=0;                                              \
            return StatusCode::FAILURE;                                 \
        }                                                               \
    } while(false)
#define REG2D(type,p,name,title,nx,minx,maxx,ny,miny,maxy)              \
    do {                                                                \
        p = new type(name,title,nx,minx,maxx,ny,miny,maxy);             \
        string path = m_filepath + name;                                \
        if (m_hsvc->regHist(path.c_str(),p).isFailure()) {              \
            error() << "Could not register \"" << path << endreq;       \
            delete p; p=0;                                              \
            return StatusCode::FAILURE;                                 \
        }                                                               \
    } while(false)

StatusCode GenHists::initialize()
{
    this->GaudiAlgorithm::initialize();

    // Get the histogram service
    if ( service("THistSvc", m_hsvc).isFailure()) {
        error() << " No THistSvc available." << endreq;
        return StatusCode::FAILURE;
    } 

    // Book some histograms, using CPP macro to help reduce typing

    REG1D(TH1I,m_nKineParts,"nKineParts","Total number of primaries",10,0,10);
    REG1D(TH1F,m_kineEnergy,"kineEnergy","Particle Energy (MeV)",100,0,m_maxEnergy);
    REG1D(TH1F,m_vtxTime,"kineVtxTime","Times of vertices (seconds)",1500,0,15*60);
    return this->BookVertices();
}

StatusCode GenHists::BookVertices()
{
    // Get volume of generation to get bounds later
    DetectorElement* detelem=0;
    if ("" == m_volume) {
        warning() << "No Volume property set, guessing at volume bounds." << endreq;
    }
    if (!existDet<DetectorElement>(m_volume)) {
        warning() << "Failed to get detector element " << m_volume 
                  << ", will guess at volume bounds" << endreq;
    }
    else {
        detelem = getDet<DetectorElement>(m_volume);
    }

    double minx = -25*CLHEP::meter, maxx = +25*CLHEP::meter;
    double miny = -25*CLHEP::meter, maxy = +25*CLHEP::meter;
    double minz = -25*CLHEP::meter, maxz = +25*CLHEP::meter;
    if (detelem) {
        m_geo = detelem->geometry();
        const ILVolume* lvol = m_geo->lvolume();
        if (!lvol) {
            error() << "Could not get logical volume for: \"" << m_volume << "\"" << endreq;
            return StatusCode::FAILURE;
        }

        const SolidBase* solid = dynamic_cast<const SolidBase*>(lvol->solid());
        minx = solid->xMin(); maxx = solid->xMax();
        miny = solid->yMin(); maxy = solid->yMax();
        minz = solid->zMin(); maxz = solid->zMax();        
        info () << "Volume bounded by: "
                << "X: [" << minx/CLHEP::meter << " " << maxx/CLHEP::meter << "], "
                << "Y: [" << miny/CLHEP::meter << " " << maxy/CLHEP::meter << "], "
                << "Z: [" << minz/CLHEP::meter << " " << maxz/CLHEP::meter << "]"
                << endreq;
    }
    REG2D(TH2F,m_xy,"kineVertexXY","Primary Vertices, X-Y",
          100,minx/CLHEP::meter,maxx/CLHEP::meter,100,miny/CLHEP::meter,maxy/CLHEP::meter);
    REG2D(TH2F,m_yz,"kineVertexYZ","Primary Vertices, Y-Z",
          100,miny/CLHEP::meter,maxy/CLHEP::meter,100,minz/CLHEP::meter,maxz/CLHEP::meter);
    REG2D(TH2F,m_xz,"kineVertexXZ","Primary Vertices, Z-X",
          100,minx/CLHEP::meter,maxx/CLHEP::meter,100,minz/CLHEP::meter,maxz/CLHEP::meter);

    return StatusCode::SUCCESS;
}

StatusCode GenHists::execute()
{
    DayaBay::GenHeader* header = get<DayaBay::GenHeader>(m_location);
    
    HepMC::GenEvent* genevt = header->event();
    if (!genevt) { 
        warning() << "No event, giving up" << endreq;
        // we don't return failure as it is not fatal.
        return StatusCode::SUCCESS; 
    }

    HepMC::GenVertex* genvtx = genevt->signal_process_vertex();
    m_nKineParts->Fill(genvtx->particles_out_size());

    HepMC::GenVertex::particles_out_const_iterator 
        it = genvtx->particles_out_const_begin(), 
        done = genvtx->particles_out_const_end();

    double energy=0;
    for (; it != done; ++it) {
        HepMC::GenParticle* part = *it;
        HepMC::FourVector momentum = part->momentum();
        energy += momentum.e();
        m_kineEnergy->Fill(momentum.e()/m_energyUnits);
    }
    
    HepMC::FourVector r4 = genvtx->position();

    Gaudi::XYZPoint gp(r4.x(),r4.y(),r4.z());
    Gaudi::XYZPoint lp = m_geo->toLocal(gp);

    m_vtxTime->Fill(r4.t()/CLHEP::second);
    m_xy->Fill(lp.x()/CLHEP::meter,lp.y()/CLHEP::meter);
    m_yz->Fill(lp.y()/CLHEP::meter,lp.z()/CLHEP::meter);
    m_xz->Fill(lp.x()/CLHEP::meter,lp.z()/CLHEP::meter);

    return StatusCode::SUCCESS;
}

StatusCode GenHists::finalize()
{

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

---