BDTestData.cxx

Go to the documentation of this file.

#include "BDTestData.h"

#include <BeamDataUtil/BDEarliest.h>
#include <BeamDataUtil/BDTarget.h>

#include <Validity/VldContext.h>
#include <RawData/RawRecord.h>
#include <RawData/RawBeamMonHeaderBlock.h>
#include <RawData/RawBeamMonBlock.h>
#include <RawData/RawDataBlock.h>
#include <RawData/RawBeamData.h>

#include <SpillTiming/SpillTimeFinder.h>

#include <TFile.h>
#include <TTree.h>
#include <TSystem.h>

#include <iostream>
#include <vector>
#include <list>
using namespace std;


BDTestData::BDTestData()
{
    this->Reset();
}
BDTestData::~BDTestData() {}
void BDTestData::Reset() 
{
    earliest = 0.0;
    tortgt = trtgtd = tor101 = tr101d = 0.0;
    nbunch = 0;
    tgtdist = -99999;
    horni = 0;
    for (int ind=0; ind<8; ++ind) {
        hitgt[ind] = vitgt[ind] = 0.0;
        xtgt[ind] = ytgt[ind] = -99999;
    }
    xsigma = 0.0;
    ysigma = 0.0;
}

static void set_toroid(const char* n, float& t, const RawBeamMonBlock& block)
{
    t=0.0;
    const RawBeamData* d = block[n];
    if (d && d->GetDataLength()) t=d->GetData()[0];
    if (t>50)
        cerr << n << " " << t << endl;
}

bool fill_bdtest(BDTestData& td, const RawBeamMonBlock& block,
                 const RawBeamMonHeaderBlock& head)
{
    td.Reset();

    const int trig_event = block.TclkTriggerEvent();
    if (trig_event != 0xa9) return false;

    // time stamp
    BDEarliest earliest;
    earliest.SetSpill(head,block);
    double dae=0, vme=0;
    earliest.GetTimestamps(dae,vme);
    if (vme) td.earliest = vme;
    else td.earliest = dae;

    // toroids
    set_toroid("E:TORTGT",td.tortgt,block);
    set_toroid("E:TRTGTD",td.trtgtd,block);
    set_toroid("E:TOR101",td.tor101,block);
    set_toroid("E:TR101D",td.tr101d,block);

    BDTarget targ;
    targ.SetSpill(head,block);
    vector<double> xp,yp,xi,yi;
    targ.BpmProjection(xp,yp,xi,yi);
    td.nbunch = xp.size();
    for (int ind=0; ind<td.nbunch && ind<8; ++ind) {
        td.hitgt[ind] = static_cast<float>(xi[ind]);
        td.vitgt[ind] = static_cast<float>(yi[ind]);
        td.xtgt[ind] =  static_cast<float>(xp[ind]);
        td.ytgt[ind] =  static_cast<float>(yp[ind]);
    }
    double xm=0, ym=0, xsig=0, ysig=0;

    targ.ProfileProjection(xm,ym,xsig,ysig);
    td.xsigma = xsig; 
    td.ysigma = ysig; 
    bool is_in=false;
    double z_loc = 0;
    BDTarget::BeamType beam_type = targ.TargetIn(is_in,z_loc);
    if (beam_type == BDTarget::kUnknown || !is_in) {
        cerr << "Unknown target location: is_in="<<is_in<<" at " << z_loc << endl;
    }
    if (!is_in)
        td.tgtdist = -1e6;
    else td.tgtdist = z_loc;

    // horn
    const char* dev[] = {"E:NSLINA","E:NSLINB","E:NSLINC","E:NSLIND",0};
    float horni = 0;
    for (int ind=0; dev[ind]; ++ind) {
        const RawBeamData* d = block[dev[ind]];
        if (d && d->GetDataLength()) horni += d->GetData()[0];
    }
    td.horni = horni;

    td.spilltimend =
      SpillTimeFinder::Instance().GetTimeOfNearestSpill(head.GetVldContext());

    return true;
}

void convert_raw_to_bdtest(const char* raw_file, const char* output_file)
{
    TFile infile(raw_file,"READ");
    TTree* intree=(TTree*)(infile.Get("BeamMon"));

    TFile outfile(output_file,"RECREATE");
    TTree* outtree = new TTree("bdtest", "Beam Monitoring Data");
    BDTestData* bdtd = new BDTestData;
    outtree->Branch("pot","BDTestData",&bdtd);

    RawRecord* record = 0;
    int nentries = intree->GetEntries();
    for ( Int_t ient = 0; ient < nentries; ient++ ) {
        intree -> SetBranchAddress("RawRecord",&record);
        intree->GetEntry(ient);

        const RawBeamMonBlock* rbmb = 0;
        const RawBeamMonHeaderBlock* rbmhb = 0;
        TObject *obj = 0;
        TIter itr = record->GetRawBlockIter();
        while ((obj = itr())) {
            const RawBeamMonBlock* b = dynamic_cast<const RawBeamMonBlock*>(obj);
            if (b) rbmb = b;
            const RawBeamMonHeaderBlock* h = dynamic_cast<const RawBeamMonHeaderBlock*>(obj);
            if (h) rbmhb = h;
        }

        if (!(rbmb && rbmhb)) continue;

        if (fill_bdtest(*bdtd,*rbmb,*rbmhb))
            outtree->Fill();
    }

    outfile.cd();
    outtree->Write();
    outfile.Close();

    infile.Close();
}

#include <TGraph.h>
#include <TCanvas.h>
#include <TH1F.h>
#include <TAxis.h>
#include <TGaxis.h>

class BDTDPlotter {
    TGraph inst_all,inst_le,inst_pme,inst_phe;
    TGraph accu_all,accu_le,accu_pme,accu_phe;
    double count_all, count_le, count_pme, count_phe;

public:
    BDTDPlotter();

    void fill(BDTestData& d);
    void write(const char* name);
};

BDTDPlotter::BDTDPlotter()
{
    count_all = count_le = count_pme = count_phe = 0.0;
}

void BDTDPlotter::fill(BDTestData& bdtd)
{
    double pot = 0.0;
    if (pot <= 0.0) pot = bdtd.trtgtd;
    if (pot <= 0.0) pot = bdtd.tortgt;
    if (pot <= 0.0) pot = bdtd.tr101d;
    if (pot <= 0.0) pot = bdtd.tor101;
    
    double ts = bdtd.earliest;

    if (ts == 0.0) {
        cerr << "Skipping zero timestamp - pot=" << pot << endl;
        return;
    }
    if (pot < 0.0) {
        if (pot > -0.1) return;
        TDatime dt((time_t)ts);
        cerr << "Skipping negative PoT = " << pot
             << " at " << Form("%.6f",ts)
             << " " << dt.AsSQLString()
             << endl;
        return;
    }

    inst_all.SetPoint(inst_all.GetN(),ts,pot);

    double accum_pot=0;
    if (!accu_all.GetN()) {
        accu_all.SetPoint(0,ts,pot);
        accum_pot=pot;
    }
    else {
        double accum_ts=0;
        accu_all.GetPoint(accu_all.GetN()-1,accum_ts,accum_pot);
        accum_pot += pot;
        accu_all.SetPoint(accu_all.GetN(),ts,accum_pot);
    }    
    count_all += pot;

    double tgtdist = bdtd.tgtdist;

    // Ripped from BDTarget::TargetIn

    const float dist_cut = 0.2*Munits::meter;
    const float le_pos   = 0.0*Munits::meter;
    const float pme_pos  = 1.0*Munits::meter;
    const float phe_pos  = 2.5*Munits::meter;
    if (fabs(tgtdist-le_pos) < dist_cut) { // kLE;
        inst_le.SetPoint(inst_le.GetN(),ts,pot);
        accu_le.SetPoint(accu_le.GetN(),ts,accum_pot);
        count_le += pot;
    }
    else if (fabs(tgtdist-pme_pos) < dist_cut) { // kPsME;
        inst_pme.SetPoint(inst_pme.GetN(),ts,pot);
        accu_pme.SetPoint(accu_pme.GetN(),ts,accum_pot);
        count_pme += pot;
    }
    else if (fabs(tgtdist-phe_pos) < dist_cut) { // kPsHE;
        inst_phe.SetPoint(inst_phe.GetN(),ts,pot);
        accu_phe.SetPoint(accu_phe.GetN(),ts,accum_pot);
        count_phe += pot;
    }

    if (pot > 50)
        cerr << accu_all.GetN()
             << ": pot=" << pot << " accum_pot=" << accum_pot
             << " @ " << Form("%.6f",ts)
             << endl;
}

#include <TDatime.h>
#include <TText.h>

void footer()
{
    TDatime t;
    double x=0.01, y=0.01;
    const char* str = Form("Brett Viren, %d/%02d/%02d",t.GetYear(),t.GetMonth(),t.GetDay());
    TText tt(x,y,str);
    float siz = tt.GetTextSize();
    tt.SetTextSize(0.5*siz);
    tt.DrawTextNDC(x,y,str);

}

#include <TVirtualPad.h>
#include <TMarker.h>
#include <TPaveText.h>
#include <TImageDump.h>

void legend(double X1, double Y1, double X2, double Y2,
            int nlines, int colors[], string keys[],
            char type = 'l', int padcolor=0);

void legend(double X1, double Y1, double X2, double Y2,
            int nlines, int colors[], string keys[],
            char type, int padcolor)
{
    TVirtualPad* oldpad = gPad;

    TPad* newpad = new TPad("newpad","Legend",X1,Y1,X2,Y2,0,1);
    newpad->SetBorderSize(1);
    newpad->SetBorderMode(0);
    newpad->SetFillColor(padcolor);
    newpad->SetFillStyle(4090);
    newpad->Draw();
    newpad->cd();

    double t = 1.0/nlines;

    int ind;
    for (ind=0; ind < nlines; ++ind) {
        double y2   = (1.0 + ind)*t;
        double ymid = (0.5 + ind)*t;
        double y1   = (0.0 + ind)*t;

        switch (type) {
        case 'm': {             // markers
            TMarker marker;
            marker.SetMarkerColor(colors[ind]);
            if (colors[ind] == 1)
                marker.SetMarkerStyle(8);
            else
                marker.SetMarkerStyle(colors[ind]);
            marker.DrawMarker(0.125,ymid);
            break;
        }
        case 'l': default: {    // lines
            TLine line;
            line.SetLineColor(colors[ind]);
            //line.SetLineStyle(colors[ind]);
            line.SetLineWidth(5);
            line.DrawLine(0.0,ymid,0.25,ymid);
            break;
        }
        }
        TPaveText* pt = new TPaveText(0.25,y1,1.0,y2,"NDC");
        pt->SetBorderSize(0);
        pt->SetFillColor(0);
        pt->SetFillStyle(0);
        pt->AddText(keys[ind].c_str());
        pt->Draw();

        cerr << ind << ": " << colors[ind] << ", " << keys[ind] << endl;
    }
    newpad->Modified();
    newpad->Update();

    oldpad->cd();
}

void BDTDPlotter::write(const char* /*filename*/)
{
    // mostly ripped from BeamData/ana/npot.C

    const float min = 0.0;      // *scale for instantaneous
    const float max = 40.0;
    TCanvas* c1 = new TCanvas("c1","npot",640,480);

    TPad* pad = new TPad("pad","",0,0,1,1);
    //pad->SetFillColor(42);
    pad->SetGrid();
    pad->Draw();
    pad->cd();

    TGraph* inst = &inst_all;
    TH1F* hinst = inst->GetHistogram();
    TAxis *ainst = hinst->GetXaxis();

    TH1F* hfinst = c1->DrawFrame(ainst->GetXmin(), min, ainst->GetXmax(), max);
    TAxis* fainst = hfinst->GetXaxis();
    fainst->SetTimeDisplay(1);
    fainst->SetTimeFormat("%m/%d");
    fainst->SetTitle("Date (2005)");
    hfinst->GetYaxis()->SetTitle("Protons per Spill (1E12)");
    int ndiv = hfinst->GetNdivisions("Y");
    hfinst->SetTitle("NuMI Protons");

    inst->Draw("B");
    c1->cd();
    TPad* overlay = new TPad("overlay","",0,0,1,1);
    overlay->SetFillStyle(4000);
    overlay->SetFillColor(0);
    overlay->SetFrameFillStyle(4000);
    overlay->Draw();
    overlay->cd();

    const float large_marker_size = 1;
    const float small_marker_size = 0.8;
    const int small_line_size = 3;
    const char* drawopt = "P";

    
    accu_all.SetLineColor(0);
    accu_all.SetLineWidth(5);
    accu_all.SetMarkerSize(large_marker_size);
    accu_all.SetMarkerStyle(20);

    int n_accumulated = accu_all.GetN();
    double maxaccum = accu_all.GetY()[n_accumulated-1];
    cerr << "Max = " << maxaccum  << endl;
    double maxscale=1.60;       // highest value for accum
    int expo=20;                // accum exponent 
    double max_pad_scale = maxscale*1e8;

    cerr << "using maxscale = " << maxscale << ", expo = "
         << expo << " max_pad_scale = " << max_pad_scale << endl;

    TH1F* hfaccum = overlay->DrawFrame(fainst->GetXmin(),0,
                                       fainst->GetXmax(),max_pad_scale);
    hfaccum->SetNdivisions(0,"X");
    hfaccum->SetNdivisions(0,"Y");
    TAxis* oainst = hfaccum->GetXaxis();
    oainst->SetTimeDisplay(1);
    oainst->SetTimeFormat("%m/%d");
    oainst->SetLabelOffset(99);

    hfaccum->GetYaxis()->SetLabelOffset(99);

    accu_all.Draw(drawopt);

    TGaxis* axis = new TGaxis(oainst->GetXmax(), 0,
                              oainst->GetXmax(), max_pad_scale,
                              0,maxscale, ndiv,"+L");

    axis->SetTitle(Form("Total Delivered Protons (1E%d)",expo));
    axis->SetTextColor(2);
    axis->SetLineColor(2);
    axis->SetLabelColor(2);
    axis->Draw();

    inst_le.SetMarkerColor(2);
    inst_le.SetLineWidth(1);
    accu_le.SetMarkerColor(2);
    accu_le.SetLineColor(2);
    accu_le.SetLineWidth(small_line_size);
    accu_le.SetMarkerSize(small_marker_size);
    accu_le.SetMarkerStyle(20);

    inst_pme.SetMarkerColor(4);
    inst_pme.SetLineWidth(1);
    accu_pme.SetMarkerColor(4);
    accu_pme.SetLineColor(4);
    accu_pme.SetLineWidth(small_line_size);
    accu_pme.SetMarkerSize(small_marker_size);
    accu_pme.SetMarkerStyle(20);

    inst_phe.SetMarkerColor(6);
    inst_phe.SetLineWidth(1);
    accu_phe.SetMarkerColor(6);
    accu_phe.SetLineColor(6);
    accu_phe.SetLineWidth(small_line_size);
    accu_phe.SetMarkerSize(small_marker_size);
    accu_phe.SetMarkerStyle(20);

    if (inst_le.GetN()) {
        //inst_le.Draw("B");
        accu_le.Draw(drawopt);
    }

    if (inst_pme.GetN()) {
        //inst_pme.Draw("B");
        accu_pme.Draw(drawopt);
    }

    if (inst_phe.GetN()) {
        //inst_phe.Draw("B");
        accu_phe.Draw(drawopt);
    }

    footer();

    int colors[] = { 1,2,4,6 };
    string keys[] = { 
        Form("All (%5.2E PoT)",count_all*1e12),
        Form("LE  (%5.2E PoT)",count_le*1e12),
        Form("pME (%5.2E PoT)",count_pme*1e12),
        Form("pHE (%5.2E PoT)",count_phe*1e12)
    };
    legend(.2,.8,.5,.6,4,colors,keys);

    c1->Modified();
    c1->Update();


#if 0
    cerr << "Printing eps\n";
    c1->Print("npot.eps");
    cerr << "Printing pdf\n";
    c1->Print("npot.pdf");
#endif

#if 0
    cerr << "Dumping npot.png" << endl;
    TImageDump* di = new TImageDump("npot.png");
    c1->Paint();
    di->Close();
    //delete di;
#else
    cerr << "Printing npot.png" << endl;
    c1->Print("npot.png");
#endif
#if 0
    cerr << "Dumping npot.gif" << endl;
    di = new TImageDump("npot.gif");
    c1->Paint();
    di->Close();
    //delete di;
#else
    cerr << "Printing npot.gif" << endl;
    c1->Print("npot.gif");
#endif


    cerr << "Writing root file\n";
    TFile file("npot.root","RECREATE");
    file.cd();
    inst->SetName("pot");
    inst->Write();
    file.Close();

}

void generate_bdtd_plots(const char* input_directory, const char* outfile)
{
    void* dir = gSystem->OpenDirectory(input_directory);
    if (!dir) {
        cerr << "Failed to open directory \"" << input_directory << "\"\n";
        return;
    }
    
    BDTDPlotter* plotter = new BDTDPlotter; // leak it!

    list<string> file_names;
    const char* cptr=0;
    while ( (cptr = gSystem->GetDirEntry(dir)) ) if (cptr) file_names.push_back(cptr);
    if (!file_names.size()) return;
    file_names.sort();
    list<string>::iterator it, done = file_names.end();
    for (it=file_names.begin(); it != done; ++it) {
        string file = *it;
        if (string::npos == file.find(".bdtd.root")) {
            cerr << "Skipping unmatched file \"" << file << "\"\n";
            continue;
        }

        string path = input_directory;
        path += "/" + file;
        cerr << "Processing " << path << endl;
        TFile tfile(path.c_str(),"READ");
        TTree* tree = (TTree*)(tfile.Get("bdtest"));
        BDTestData* bd_test_data = 0;
        int ninds = tree->GetEntries();
        for (int ind=0; ind<ninds; ++ind) {
            tree->SetBranchAddress("pot",&bd_test_data);
            tree->GetEntry(ind);

            plotter->fill(*bd_test_data);
        }
        tfile.Close();
    }

    plotter->write(outfile);
}

ClassImp(BDTestData)

---