In This Package:

SimHists.cc

Go to the documentation of this file.

#include "SimHists.h"

#include "Event/SimHeader.h"
#include "Event/SimHitHeader.h"
#include "Event/SimHitCollection.h"
#include "Event/SimHit.h"

#include "Conventions/Detectors.h"
#include "Conventions/Site.h"
#include "Conventions/DetectorId.h"

#include "GaudiKernel/ITHistSvc.h"

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

#include <map>
#include <string>
using namespace std;

typedef map<short int, pair<int,string> > detector_bins_t;
static detector_bins_t detector_bins;

SimHists::SimHists(const std::string& name, ISvcLocator* pSvcLocator)
    : GaudiAlgorithm(name,pSvcLocator)
    , m_hsvc(0)
    , m_nHitCollections(0)
    , m_nHitByDetector(0)
{
    declareProperty("Location",m_location=DayaBay::SimHeaderLocation::Default,
                    "Location in the TES to get HepMCEvents");
    declareProperty("FilePath",m_filepath="/file1/sim/",
                    "File path with with to register histograms.");

    if (!detector_bins.size()) {
        DayaBay::Detector dets[] = {
            DayaBay::Detector(Site::kDayaBay,DetectorId::kAD1),
            DayaBay::Detector(Site::kDayaBay,DetectorId::kAD2),
            DayaBay::Detector(Site::kDayaBay,DetectorId::kOWS),
            DayaBay::Detector(Site::kDayaBay,DetectorId::kIWS),
            DayaBay::Detector(Site::kDayaBay,DetectorId::kRPC),
            DayaBay::Detector(Site::kLingAo,DetectorId::kAD1),
            DayaBay::Detector(Site::kLingAo,DetectorId::kAD2),
            DayaBay::Detector(Site::kLingAo,DetectorId::kOWS),
            DayaBay::Detector(Site::kLingAo,DetectorId::kIWS),
            DayaBay::Detector(Site::kLingAo,DetectorId::kRPC),
            DayaBay::Detector(Site::kFar,DetectorId::kAD1),
            DayaBay::Detector(Site::kFar,DetectorId::kAD2),
            DayaBay::Detector(Site::kFar,DetectorId::kAD3),
            DayaBay::Detector(Site::kFar,DetectorId::kAD4),
            DayaBay::Detector(Site::kFar,DetectorId::kOWS),
            DayaBay::Detector(Site::kFar,DetectorId::kIWS),
            DayaBay::Detector(Site::kFar,DetectorId::kRPC),
            DayaBay::Detector(),
        };

        detector_bins[0] = pair<int,string>(0,"unknown");
        for (int ind=0; dets[ind].site(); ++ind) {
            detector_bins[dets[ind].siteDetPackedData()] 
                = pair<int,string>(ind+1,dets[ind].detName());
        }
    }
}

SimHists::~SimHists()
{
}

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

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

    m_nHitCollections = new TH1I("nHitCollections","Number of Hit Collections",15,0,15);
    if (m_hsvc->regHist(m_filepath+"nHitCollections",m_nHitCollections).isFailure()) {
        error() << "Could not register " << m_filepath+"nHitCollections" << endreq;
        delete m_nHitCollections; m_nHitCollections = 0;
        return StatusCode::FAILURE;
    }

    m_nHitByDetector = new TH2F("nHitByDetector","Number of hits in each detector",
                                15,0,15,1000,0,1000);
    if (m_hsvc->regHist(m_filepath+"nHitByDetector",m_nHitByDetector).isFailure()) {
        error() << "Could not register " << m_filepath+"nHitByDetector" << endreq;
        delete m_nHitByDetector; m_nHitByDetector = 0;
        return StatusCode::FAILURE;
    }
    // FIXME: label bins by detector name

    return StatusCode::SUCCESS;
}

StatusCode SimHists::execute()
{
    DayaBay::SimHeader* shead = get<DayaBay::SimHeader>(m_location);
    const DayaBay::SimHitHeader* shith = shead->hits();
    const DayaBay::SimHitHeader::hc_map& hcmap = shith->hitCollection();

    m_nHitCollections->Fill(hcmap.size());
    
    DayaBay::SimHitHeader::hc_map::const_iterator it, done = hcmap.end();
    for (it=hcmap.begin(); it != done; ++it) {
        DayaBay::Detector det(it->first);
        pair<int,string> bin = detector_bins[det.siteDetPackedData()];

        info () << "Got hit collection from " << det.detName() 
                << " in bin " << bin.first << endreq;

        size_t siz = it->second->collection().size();
        m_nHitByDetector->Fill(bin.first + 0.5,siz);
    }

    return StatusCode::SUCCESS;
}

StatusCode SimHists::finalize()
{

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

---