BDataQualityModule.cxx

Go to the documentation of this file.

#include "BDataQualityModule.h"

#include <BeamDataUtil/BDUniquify.h>
#include <BeamDataUtil/BDEarliest.h>
#include <BeamDataUtil/BDSpillAccessor.h>
#include <HistMan/HistMan.h>
#include <SpillTiming/SpillTimeFinder.h>
#include <DataUtil/GetRawBlock.h>
#include <DataUtil/DumpMom.h>
#include <RawData/RawBeamMonBlock.h>
#include <RawData/RawBeamMonHeaderBlock.h>
#include <RawData/RawBeamData.h>
#include <RawData/RawBeamSwicData.h>

#include <JobControl/JobCModuleRegistry.h>
#include <MessageService/MsgService.h>

#include <TSystem.h>

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


CVSID("$Id: BDataQualityModule.cxx,v 1.5 2007/11/11 08:16:11 rhatcher Exp $");
JOBMODULE(BDataQualityModule,"BDataQuality","BDP<-->STND data quality histogrammer");

const char* file_hist_dir = "BDP Data Quality/File Hists";
const char* db_hist_dir = "BDP Data Quality/DB Hists";

BDataQualityModule::BDataQualityModule()
    : fFirstSpillTime(0)
    , fLastSpillTime(0)
{
}
BDataQualityModule::~BDataQualityModule()
{
}

const char* tor_names[] = {
    "E:TORTGT","E:TRTGTD","E:TOR101","E:TR101D", 0
};

const char* important_devs[] = {
    "E:TORTGT","E:TRTGTD","E:TOR101","E:TR101D", 
    "E:HP121","E:VP121","E:HPTGT","E:VPTGT",
    "E:MTGTDS",
    "E:HADMDS","E:MMA1DS","E:MMA2DS","E:MMA3DS",
    "E:NSLINA","E:NSLINB","E:NSLINC","E:NSLIND",
    "I:NUTARZ",0
};

const char* some_swics[] = {
    "E:MTGTDS",
    "E:HADMDS","E:MMA1DS","E:MMA2DS","E:MMA3DS",
    0
};

JobCResult BDataQualityModule::Ana(const MomNavigator *mom)
{
    static BDUniquify uniquify;
    static int count=-1, mult=1;
    ++count;
    
    if (count>mult*10) {
        mult *= 10;
    }
    
    if (count%mult == 0)
        cerr << count << endl;

    
    if (!mom) return JobCResult::kFailed;

    HistMan hm(file_hist_dir);

    int ntotalblocks = DataUtil::GetRawBlocks<RawBeamMonBlock>(mom).size();
    hm.Fill1d("total_blocks",ntotalblocks);

    vector<BDUniquify::BlockPair> unique = uniquify.GetUniqueBlocks(*mom);
    size_t nblocks = unique.size();    
    
    if (nblocks && fFirstSpillTime == 0.0)
        fFirstSpillTime = unique[0].first->GetTimeStamp();
    if (nblocks)
        fLastSpillTime = unique[0].first->GetTimeStamp();

    if (nblocks != 1) {
        cerr << "Other than a single block in record " << count << ", " << nblocks << endl;
        DataUtil::dump_mom(mom,cerr);
    }

    hm.Fill1d("unique_blocks",nblocks);

    for (size_t ind=0; ind<nblocks; ++ind) {
        this->FillFile(*unique[ind].first,*unique[ind].second);
    }

    return JobCResult::kAOK;
}



void BDataQualityModule::BeginFile()
{
    HistMan hm(file_hist_dir);

    // Per MOM
    const char* hname = "unique_blocks";
    TH1I* hi = hm.Get<TH1I>(hname);
    if (hi) hi->Clear();
    else hm.Book<TH1I>(hname,"Unique BeamMonBlocks in RawRecord (should always be 1)",5,0,5);
    hname = "total_blocks";
    hi = hm.Get<TH1I>(hname);
    if (hi) hi->Clear();
    else hm.Book<TH1I>(hname,"Total BeamMonBlocks in RawRecord (should always be 1)",5,0,5);


    // Per block
    hname = "earliest_dae_stnd";
    TH1F* hf = hm.Get<TH1F>(hname);
    if (hf) hf->Clear();
    else hm.Book<TH1F>(hname,"Earliest DAE - STND (should be < 0.5)",1000,-1,1);

    hname = "earliest_vme_stnd";
    hf = hm.Get<TH1F>(hname);
    if (hf) hf->Clear();
    else hm.Book<TH1F>(hname,"Earliest VME - STND (should be < 0.5)",1000,-1,1);
    
    hname = "earliest_last_dae";
    hf = hm.Get<TH1F>(hname);
    if (hf) hf->Clear();
    else hm.Book<TH1F>(hname,"Rep period = Earliest - Last DAE (typically ~2sec)",
                       1000,0,10);

    hname = "dae_vme_dt";
    hf = hm.Get<TH1F>(hname);
    if (hf) hf->Clear();
    else hm.Book<TH1F>(hname,"DAE - VME time over TGT-Prof,Had,Mu monitors (should be ~0.5 sec)",
                       1000,-1,1);

    hname = "dae_vme_dt_min";
    TH2F *dae_vme_dt_min = hm.Get<TH2F>(hname);
    if (dae_vme_dt_min) dae_vme_dt_min->Clear();
    else dae_vme_dt_min =
             hm.Book<TH2F>(hname,"DAE - VME dt per min for TGT-PM,Had,Mu (>1.0 implies lost data)",
                           12*60,0,12*60, 100,-5,5);

    //

    hname = "spills_per_hour";
    hf = hm.Get<TH1F>(hname);
    if (hf) hf->Clear();
    else  hm.Book<TH1F>(hname,"Spills per hour",12,0,12);

    hname = "npot_per_hour";
    hf = hm.Get<TH1F>(hname);
    if (hf) hf->Clear();
    else hm.Book<TH1F>(hname,"Npot per hour",12,0,12);

    hname = "spills_per_min";
    hf = hm.Get<TH1F>(hname);
    if (hf) hf->Clear();
    else  hm.Book<TH1F>(hname,"Spills per min",12*60,0,12*60);

    hname = "npot_per_min";
    hf = hm.Get<TH1F>(hname);
    if (hf) hf->Clear();
    else hm.Book<TH1F>(hname,"Npot per min",12*60,0,12*60);
    
    hname = "toroids";
    hf = hm.Get<TH1F>(hname);
    if (hf) hf->Clear();
    else hf = hm.Book<TH1F>(hname,"Toroids",4,0,4);
    for (int ind=0; tor_names[ind]; ++ind) {
        hf->GetXaxis()->SetBinLabel(ind+1,tor_names[ind]);
    }

    hname = "tortgt_trtgtd";
    hf = hm.Get<TH1F>(hname);
    if (hf) hf->Clear();
    else hm.Book<TH1F>(hname,"(TORTGT - TRTGTD)/(TORTGT + TRTGTD)",100,-1,1);

    hname = "tor101_tr101d";
    hf = hm.Get<TH1F>(hname);
    if (hf) hf->Clear();
    else hm.Book<TH1F>(hname,"(TOR101 - TR101D)/(TOR101 + TR101D)",100,-1,1);

    hname = "tortgt_tor101";
    hf = hm.Get<TH1F>(hname);
    if (hf) hf->Clear();
    else hm.Book<TH1F>(hname,"(TORTGT - TOR101)/(TORTGT + TOR101)",100,-1,1);

    hname = "trtgtd_tr101d";
    hf = hm.Get<TH1F>(hname);
    if (hf) hf->Clear();
    else hm.Book<TH1F>(hname,"(TRTGTD - TR101D)/(TRTGTD + TR101D)",100,-1,1);
    

    int ndevs = 0;
    while (important_devs[ndevs]) ++ndevs;

    hname = "dev_count";
    hf = hm.Get<TH1F>(hname);
    if (hf) hf->Clear();
    else hf = hm.Book<TH1F>(hname,"Device counts",ndevs,0,ndevs);
    for (int ind=0; ind<ndevs; ++ind)
        hf->GetXaxis()->SetBinLabel(ind+1,important_devs[ind]);
}

void BDataQualityModule::EndFile()
{
    HistMan hm(db_hist_dir);

    const char* hname = "";

    hname = "stnd_bdp_dt";
    TH1F *stnd_bdp_dt = hm.Get<TH1F>(hname);
    if (stnd_bdp_dt) stnd_bdp_dt->Clear();
    else stnd_bdp_dt =
             hm.Book<TH1F>(hname,"SpillTimeND - BDP Time (>1.0 implies lost data)",
                           1000,-1,1);

    hname = "stnd_bdp_dt_hour";
    TH2F *stnd_bdp_dt_hour = hm.Get<TH2F>(hname);
    if (stnd_bdp_dt_hour) stnd_bdp_dt_hour->Clear();
    else stnd_bdp_dt_hour =
             hm.Book<TH2F>(hname,"SpillTimeND - BDP Time per hour (>1.0 implies lost data)",
                           12,0,12, 100,-5,5);

    hname = "stnd_bdp_dt_min";
    TH2F *stnd_bdp_dt_min = hm.Get<TH2F>(hname);
    if (stnd_bdp_dt_min) stnd_bdp_dt_min->Clear();
    else stnd_bdp_dt_min =
             hm.Book<TH2F>(hname,"SpillTimeND - BDP Time per min (>1.0 implies lost data)",
                           12*60,0,12*60, 100,-5,5);
    
    // also need to check how many spills are in STND between the end
    // of the last file and the begin of this one.

    Detector::Detector_t det = Detector::kNear;
    SimFlag::SimFlag_t simflag = SimFlag::kData;
    VldTimeStamp vts(fFirstSpillTime);
    VldTimeStamp numi_start(2005,1,1,0,0,0);

    SpillTimeFinder& stf = SpillTimeFinder::Instance();
    BDSpillAccessor& sa = BDSpillAccessor::Get();

    while ((double)vts <= fLastSpillTime + 1.0) {
        VldContext vc(det,simflag,vts);
        VldTimeStamp stnd_vts = stf.GetTimeOfNextSpill(vc);
        
        if (stnd_vts < numi_start) {
            cerr << "Bailing on file due to getting early SpillTimeND: "
                 << stnd_vts
                 << " last vts = " << vts
                 << endl;
            break;
        }

        const BeamMonSpill* bms = sa.LoadSpill(stnd_vts);
        VldTimeStamp bms_vts;
        if (!bms) {
            cerr << "No BeamMonSpill for " << stnd_vts << endl;
            bms_vts = VldTimeStamp::GetEOT();
        }
        else {
            bms_vts = bms->SpillTime();
        }

        double dt = stnd_vts-bms_vts;

        stnd_bdp_dt->Fill(dt);
        stnd_bdp_dt_hour->Fill(((double)stnd_vts-fFirstSpillTime)/(60*60),dt);
        stnd_bdp_dt_min->Fill(((double)stnd_vts-fFirstSpillTime)/(60),dt);

        if (vts >= stnd_vts) {
            cerr << "Error, infinite loop, bailing with\n\tbms= " << vts
                 << " stnd= " << stnd_vts << endl;
            break;
        }
        vts = stnd_vts;

    }
    
    
}

void BDataQualityModule::FillFile(const RawBeamMonHeaderBlock& head,
                                  const RawBeamMonBlock& block)
{
    static BDEarliest bde;
    bde.SetSpill(head,block);
    double dae=0,vme=0,last_dae=0,last_vme=0;
    bde.GetTimestamps(dae,vme);
    bde.GetLastTimestamps(last_dae,last_vme);


    Detector::Detector_t det = Detector::kNear;
    SimFlag::SimFlag_t simflag = SimFlag::kData;
    VldTimeStamp dae_vts(dae);
    VldContext vc(det,simflag,dae_vts);

    SpillTimeFinder& stf = SpillTimeFinder::Instance();
    VldTimeStamp stnd_vts = stf.GetTimeOfNearestSpill(vc);
    double stnd = stnd_vts;

    HistMan hm(file_hist_dir);

    const char* hname = "earliest_dae_stnd";
    hm.Fill1d(hname,dae-stnd);

    hname = "earliest_vme_stnd";
    hm.Fill1d(hname,vme-stnd);
    
    hname = "earliest_last_dae";
    hm.Fill1d(hname,dae-last_dae);
    

    hname = "dae_vme_dt";
    TH1F* dae_vme_dt =hm.Get<TH1F>(hname);
    hname = "dae_vme_dt_min";
    TH2F* dae_vme_dt_min = hm.Get<TH2F>(hname);

    RawBeamSwicData swic;
    for (int ind=0; some_swics[ind]; ++ind) {
        const RawBeamData* rbd = block[some_swics[ind]];
        if (!rbd) continue;
        VldTimeStamp dae_vts(rbd->GetSeconds(), 1000000*rbd->GetMsecs());

        swic.SetData(*rbd);
        VldTimeStamp vme_vts(swic.VmeSeconds(), swic.VmeNanoseconds());

        double dt = (double)(dae_vts - vme_vts);

        if (dae_vme_dt)
            dae_vme_dt->Fill(dt);

        if (dae_vme_dt_min)
            dae_vme_dt_min->Fill((dae-fFirstSpillTime)/60.,dt);
    }

    hname = "spills_per_hour";
    hm.Fill1d(hname,(dae-fFirstSpillTime)/(60*60));

    hname = "spills_per_min";
    hm.Fill1d(hname,(dae-fFirstSpillTime)/(60));

    double pot=0, pots[4] = {0};

    for (int ind=0; tor_names[ind]; ++ind) {
        const RawBeamData* tor_dev = block[tor_names[ind]];
        if (tor_dev && tor_dev->GetDataLength())
            pots[ind] = tor_dev->GetData()[0];
        if (pot == 0) pot = pots[ind];
        hm.Fill1d("toroids",ind,pots[ind]);
    }
    

    hname = "npot_per_hour";
    hm.Fill1d(hname,(dae-fFirstSpillTime)/(60*60),pot);

    hname = "npot_per_min";
    hm.Fill1d(hname,(dae-fFirstSpillTime)/(60),pot);

    const char* hnames[] = {
        "tortgt_trtgtd",//0,1
        "tor101_tr101d",//2,3
        "tortgt_tor101",//0,2
        "trtgtd_tr101d",//1,3
        0
    };
    for (int ind=0; ind<2; ++ind) {
        double sum = pots[ind*2]+pots[ind*2+1];
        TH1F* h = hm.Get<TH1F>(hnames[ind]);
        if (sum == 0.0)
            h->Fill(-1);
        else
            h->Fill((pots[ind*2]-pots[ind*2+1])/sum);

        sum = pots[ind]+pots[ind+2];
        h = hm.Get<TH1F>(hnames[ind+2]);
        if (sum == 0.0)
            h->Fill(-1);
        else
            h->Fill((pots[ind]-pots[ind+2])/sum);

    }

    hname = "dev_count";
    TH1F* id = hm.Get<TH1F>(hname);
    if (id) {
        for (int ind=0; important_devs[ind]; ++ind) {
            const RawBeamData* dev = block[important_devs[ind]];
            if (dev && dev->GetDataLength())
                id->Fill(ind);
            else 
                cerr << "Failed to get " << important_devs[ind] << endl;
        }
    }
    else cerr << "Failed to get dev_count\n";

}

    

---