plots.cxx

Go to the documentation of this file.

// Make some plots with the first beam data
// usage: root plots.C
// must run from the directory that contains
// - first-data.tree.root
// - {ped,swic}Constants.txt
// ie, BeamData/ana

// Convert data files to flatter TTrees.

#include <HistMan/HistMan.h>
#include <BeamData/ana/bmnt/ProfMon.h>
#include <BeamData/ana/bmnt/Knot.h>
#include <BeamData/ana/bmnt/Pedestals.h>
#include <BeamData/ana/bmnt/AcnetDevice.h>

#include <TCanvas.h>
#include <TSystem.h>
#include <TTree.h>
#include <TFile.h>
#include <TH1D.h>
#include <TH2D.h>

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

const double MIN_ETOR101 = 0.15;



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

const char* profile_monitor[] = {
    "E:M101DS",
    "E:M105DS",
    "E:M107DS",
    "E:M108DS",
    "E:M112DS",
    "E:M114DS",
    "E:M115DS",
    "E:M117DS",
    "E:M121DS",
    "E:MTGTDS",
    0
};

void book_once()
{
    HistMan hm("bv");
    // ...

    hm.Book<TH1D>("max_dae_dt","Maximum difference between 2 DAEs in data block",1200,-60,60);

}
void fill(Knot& k, ProfMon& , Pedestals& , int )
{
    double max_dt = 0;
    
    for (int ind1=0; profile_monitor[ind1]; ++ind1) {
        const AcnetDevice* ad1 = k.GetDevice(profile_monitor[ind1]);
        for (int ind2=ind1; profile_monitor[ind2]; ++ind2) {
            const AcnetDevice* ad2 = k.GetDevice(profile_monitor[ind2]);
            double diff = fabs(ad1->timestamp - ad2->timestamp);
            if (diff > max_dt) max_dt = diff;
        }
    }
    HistMan hm("bv");
    hm.Fill1d("max_dae_dt",max_dt);
    
}

void fill_profile(ProfMon& pm, int count,
                  TH1D& xh, TH1D& yh, TH2D& xa, TH2D& ya)
{
    for (int ch=pm.MIN_CHANNEL; ch<=pm.MAX_CHANNEL; ++ch) {
        xh.Fill(pm.WirePosition(ch),pm.GetVoltage(pm.Xindex(ch)));
        xa.Fill(count,pm.WirePosition(ch),pm.GetVoltage(pm.Xindex(ch)));
        yh.Fill(pm.WirePosition(ch),pm.GetVoltage(pm.Yindex(ch)));
        ya.Fill(count,pm.WirePosition(ch),pm.GetVoltage(pm.Yindex(ch)));

    }
}


void book_and_fill(Knot& k, ProfMon& pm, Pedestals& peds, int count)
{
    HistMan hm(Form("bv/spill%03d",count));

    double min = ProfMon::WirePosition(ProfMon::MIN_CHANNEL);
    double max = ProfMon::WirePosition(ProfMon::MAX_CHANNEL);
    min -= ProfMon::WIRE_SPACING;
    max -= ProfMon::WIRE_SPACING;
    int nbins = ProfMon::MAX_CHANNEL - ProfMon::MIN_CHANNEL + 1;

    int npms = 0;
    for (; profile_monitor[npms]; ++npms);
    TH2D* xall =
        hm.Book<TH2D>(Form("xallpm_%03d",count),
                      Form("All X Profiles for spill %d",count),
                      npms,0,npms,nbins,min,max);
    TH2D* yall =
        hm.Book<TH2D>(Form("yallpm_%03d",count),
                      Form("All Y Profiles for spill %d",count),
                      npms,0,npms,nbins,min,max);

    for (int ind=0; profile_monitor[ind]; ++ind) {
        const AcnetDevice* ad = k.GetDevice(profile_monitor[ind]);
        if (!pm.SetPeds(peds.GetPeds(profile_monitor[ind]))) {
            cerr << "Failed to set peds for " << profile_monitor[ind]
                 << " skipping\n";
            continue;
        }
        int ret = pm.SetAcnetDevice(*ad);
        if (! ret) {
            cerr << "Failed to set acnet device for spill " << count << endl;
            continue;
        }

        string pms = profile_monitor[ind];
        pms[1]='_';
        TH1D* xhist =
            hm.Book<TH1D>(Form("%sX_%03d",pms.c_str(),count),
                          Form("ProfMon %s X view spill %d",pms.c_str(),count),
                          nbins,min,max);
        TH1D* yhist =
            hm.Book<TH1D>(Form("%sY_%03d",pms.c_str(),count),
                          Form("ProfMon %s Y view spill %d",pms.c_str(),count),
                          nbins,min,max);

        fill_profile(pm,ind,*xhist,*yhist,*xall,*yall);
        ++npms;
    }

}

void write_plots()
{
    HistMan hm("bv");
    TFile file("plots.root","RECREATE");
    hm.WriteOut(file);
    file.Close();
}

static int nentries=0;

void make_plots()
{
    //TFile f("first-data.tree.root");
    TFile f("second.bmnt.root");
    TTree* bd = dynamic_cast<TTree*>(f.FindObjectAny("bd"));
    assert(bd);
    Knot k(*bd);

    book_once();

    Pedestals peds;
    peds.LoadPedFile("peds.txt");

    nentries = k.GetSize();
    int count = 0;
    for (int entry=0; entry<nentries; ++entry) {
        cerr << entry << " ";
        k.GetEntry(entry);

        const AcnetDevice* ad = k.GetDevice("E:TOR101");
        if (ad->ndata == 0) {
            cerr << "No E:TOR101!\n";
            continue;
        }
        double cut = fabs(ad->data[0]);
        if (cut < MIN_ETOR101) continue;
        ++count;

        ProfMon pm;
        pm.SetScale(cut/MIN_ETOR101);
        book_and_fill(k,pm,peds,count);
        fill(k,pm,peds,count);
    }
    //write();
}

void make_prints() {
    TCanvas c("c","BeamData Plots",850,1100);

    c.Print("bdplots.ps[");

    HistMan hm("bv");
    hm.Get<TH1D>("max_dae_dt")->Draw();
    c.Print("bdplots.ps");
    
#if 0
    for (int ind=1; ind <= nentries; ++ind) {
        HistMan hm(Form("bv/spill%03d",ind));
        TH2D* hx = hm.Get<TH2D>(Form("xallpm_%03d",ind));
        TH2D* hy = hm.Get<TH2D>(Form("yallpm_%03d",ind));

        if (!hx || !hy) {
            cerr << "Failed to get hists for spill " << ind << endl;
            continue;
        }

        c.Clear();
        c.Divide(1,2);
        c.cd(1);
        hx->Draw("colz");
        c.cd(2);
        hy->Draw("colz");

        c.Print("bdplots.ps");
    }
#endif

    c.Print("bdplots.ps]");
    
}

void plots()
{
    make_plots();
    make_prints();
    write_plots();
}

---