CheckOefDbAlg.cxx

Go to the documentation of this file.

#include "EventModel/Event.h"
#include "EventModel/EventHeader.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
#include <iomanip>
#include <iostream>
 
#include "CheckOefDbAlg.h"
 
using namespace std;
DECLARE_COMPONENT( CheckOefDbAlg )
CheckOefDbAlg::CheckOefDbAlg(const std::string &name, ISvcLocator *pSvcLocator)
    : Algorithm(name, pSvcLocator) {
    declareProperty("Check_etsT1_hist", m_check_etsT1_time_only = 0); // if it is not 0, we will output a .root file named with this runID, and only check the time and etsT1.
}
 
StatusCode CheckOefDbAlg::initialize() {
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "in initialize()" << endmsg;
 
    StatusCode sc = service("InjSigIntervalSvc", m_InjSigIntervalSvc);
    if (sc != StatusCode::SUCCESS) {
        log << MSG::FATAL << "can not use InjSigIntervalSvc" << endmsg;
    }
 
    sc = service("InjSigTimeSvc", m_InjSigTimeSvc);
    if (sc != StatusCode::SUCCESS) {
        log << MSG::FATAL << "can not use InjSigTimeSvc" << endmsg;
    }
 
    sc = service("OfflineEvtFilterSvc", m_evtFilterSvc);
    if (sc != StatusCode::SUCCESS) {
        log << MSG::FATAL << "can not use OfflineEvtFilterSvc" << endmsg;
    }
    m_lastRun    = -1;
    m_hist_etsT1 = NULL;
    if (m_check_etsT1_time_only) {
 
        const int expected_range = 50000;
 
        m_hist_etsT1 = new TH1D("hist_etsT1", "hist_etsT1", expected_range * 2 + 1, -expected_range - .5, +expected_range + .5);
 
        cout << "creating new histograms" << endl;
    }
    return StatusCode::SUCCESS;
}
 
StatusCode CheckOefDbAlg::execute() {
    MsgStream log(msgSvc(), name());
    //Get EventHeader
    SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(), "/Event/EventHeader");
    if (!eventHeader) {
        log << MSG::FATAL << "Could not find Event Header" << endmsg;
        return StatusCode::FAILURE;
    }
    if (m_hist_etsT1){ // fill in the histogram
        unsigned long time_etsT1   = eventHeader->etsT1();
        long time_etsT1_timelike32 = ((time_etsT1 / 2000000) & 0x00000000ffffffffLU);
        // const int expected_range=50000;
        m_hist_etsT1->Fill(time_etsT1_timelike32);
    }
 
    if (m_lastRun == eventHeader->runNumber()) {
        return StatusCode::SUCCESS;
    }
    m_lastRun = eventHeader->runNumber();
    cout << "==============CheckOefDbAlg: Run " << m_lastRun << "====================" << endl;
 
    cout << "shielding window" << endl;
    cout << "RunFrom: " << setw(15) << m_evtFilterSvc->getRunFrom() << endl
         << "RunTo: " << setw(15) << m_evtFilterSvc->getRunTo() << endl
         << "EventFrom: " << setw(15) << m_evtFilterSvc->getEventFrom() << endl
         << "EventTo: " << setw(15) << m_evtFilterSvc->getEventTo() << endl
         << "shielding window Npar: " << setw(15) << m_evtFilterSvc->getNpar() << endl;
    int npar = m_evtFilterSvc->getNpar();
    for (int i = 0; i < npar; i++) {
        cout << "Flag:" << setw(15) << m_evtFilterSvc->getFlag(i) << endl
             << "TBegin: " << setw(15) << m_evtFilterSvc->getTBegin(i) << endl
             << "TEnd: " << setw(15) << m_evtFilterSvc->getTEnd(i) << endl;
    }
 
    cout << "injection interval" << endl;
    cout << "TInterval: " << setw(15) << m_InjSigIntervalSvc->getTInterval() << endl;
 
    cout << "IST" << endl;
    cout << "IST Npar: " << setw(15) << m_InjSigTimeSvc->getNpar() << endl;
    for (int i = 0; i < m_InjSigTimeSvc->getNpar(); i++) {
        cout << "IST Flag: " << setw(15) << m_InjSigTimeSvc->getFlag(i) << endl
             << "IST:" << setw(15) << m_InjSigTimeSvc->getIST(i) << endl;
    }
    cout << "CheckOefDbAlg: checking OEF, npar=" << npar << endl;
    for (int i = 0; i < npar; i++) {
        cout << "CheckOefDbAlg: OEF: flag:tBegin:tEnd = " << m_evtFilterSvc->getFlag(i)
             << ":" << setw(13) << m_evtFilterSvc->getTBegin(i)
             << ":" << setw(13) << m_evtFilterSvc->getTEnd(i) << endl;
    }
 
    return StatusCode::SUCCESS;
}
static std::pair<int, int> get_range_bin(TH1 *hist) {
    int nbins             = hist->GetNbinsX();
    int first_nonzero_bin = -1;
    int last_nonzero_bin  = nbins + 1;
    for (int i = 0; i <= nbins + 1; i++) {
        if (hist->GetBinContent(i) > 0) {
            first_nonzero_bin = i;
            break;
        }
    }
    for (int i = nbins + 1; i >= 0; i--) {
        if (hist->GetBinContent(i) > 0) {
            last_nonzero_bin = i;
            break;
        }
    }
    return std::make_pair(first_nonzero_bin, last_nonzero_bin);
}
static std::pair<double, double> get_new_range(TH1 *hist, std::pair<int, int> range_bin) {
    int nbins         = hist->GetNbinsX();
    int new_bin_first = range_bin.first - 5;
    if (new_bin_first < 0) new_bin_first = 0;
    int new_bin_last = range_bin.second + 5;
    if (new_bin_last > nbins + 1) new_bin_last = nbins + 1;
    return std::make_pair(hist->GetBinCenter(new_bin_first) - .6, hist->GetBinCenter(new_bin_last) + .6);
}
StatusCode CheckOefDbAlg::finalize() {
    if (m_hist_etsT1) {
        
        std::stringstream ss;
        ss << "hists_" << m_lastRun << ".root";
        cout << "Printing histogram to "<< ss.str().c_str() << endl;
        // next we try to find the where in the hist is not 0.
        // we will use this to find the range of the histogram
 
        std::pair<double, double> range_etsT1_value = get_new_range(m_hist_etsT1, get_range_bin(m_hist_etsT1));
        m_hist_etsT1->GetXaxis()->SetRangeUser(range_etsT1_value.first, range_etsT1_value.second);
 
        TFile *f = new TFile(ss.str().c_str(), "recreate");
        f->cd();
        m_hist_etsT1->Write();
        f->Close();
 
        delete f;
        delete m_hist_etsT1;
    }
    return StatusCode::SUCCESS;
}