Draw_L1_histo.C

Go to the documentation of this file.

#include <unistd.h>   // for dir navigation
 
const int textFont = 22;   // TNewRoman
const bool divide = 0;     // = 0 - each histo in separate file\screen. = 1 - all in one file\screen
 
void FitHistoGaus(TH1* hist, float* sigma = 0)
{
    if (hist && (hist->GetEntries() != 0)) {
        TF1* fit = new TF1("fit", "gaus");
        fit->SetLineColor(2);
        fit->SetLineWidth(3);
        hist->Fit("fit", "", "", hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax());
        if (sigma)
            *sigma = fit->GetParameter(2);
    } else {
        std::cout << " E: Read hists error! " << std::endl;
    }
}
 
void FitHistoLine(TH1* hist)
{
    if (hist && (hist->GetEntries() != 0)) {
        TF1* fit = new TF1("fit", "[0]");
        fit->SetLineColor(2);
        fit->SetLineWidth(3);
        hist->Fit("fit", "", "", hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax());
    } else {
        std::cout << " E: Read hists error! " << std::endl;
    }
}
 
void makeUpHisto(TH1* hist, TString title)
{
    if (hist && (hist->GetEntries() != 0)) {
 
        hist->GetXaxis()->SetLabelFont(textFont);
        hist->GetXaxis()->SetTitleFont(textFont);
        hist->GetYaxis()->SetLabelFont(textFont);
        hist->GetYaxis()->SetTitleFont(textFont);
 
        hist->GetXaxis()->SetTitle(title);
        hist->GetXaxis()->SetTitleOffset(1);
        hist->GetYaxis()->SetTitle("Entries");
        hist->GetYaxis()->SetTitleOffset(1.05);   // good then entries per bit <= 9999
    } else {
        std::cout << " E: Read hists error! " << std::endl;
    }
}
 
void makeUpHisto(TH2* hist, TString title)
{
    if (hist && (hist->GetEntries() != 0)) {
        // hist->GetXaxis()->SetLabelFont(textFont);
        // hist->GetXaxis()->SetTitleFont(textFont);
        // hist->GetYaxis()->SetLabelFont(textFont);
        // hist->GetYaxis()->SetTitleFont(textFont);
 
        // hist->GetXaxis()->SetTitle(title);
        // hist->GetXaxis()->SetTitleOffset(1);
        // hist->GetYaxis()->SetTitle("Entries");
        // hist->GetYaxis()->SetTitleOffset(1.05); // good then entries per bit <= 9999
 
        hist->SetOption("CONTZ");
    } else {
        std::cout << " E: Read hists error! " << std::endl;
    }
}
 
void makeUpProfile(TProfile* prof, TString titleX, TString titleY)
{
    if (prof && (prof->GetEntries() != 0)) {
        prof->SetMarkerColor(kRed);
        prof->SetLineColor(kRed);
 
        prof->GetXaxis()->SetTitle(titleX);
        prof->GetXaxis()->SetTitleOffset(1);
        prof->GetYaxis()->SetTitle(titleY);
        prof->GetYaxis()->SetTitleOffset(0.85);
 
        prof->GetXaxis()->SetLabelFont(textFont);
        prof->GetXaxis()->SetTitleFont(textFont);
        prof->GetYaxis()->SetLabelFont(textFont);
        prof->GetYaxis()->SetTitleFont(textFont);
    } else {
        std::cout << " E: Read profile error! " << std::endl;
    }
}
 
void Draw_L1_histo()
{
 
    // ============================ Set Styles ============================
    TStyle* histoStyle = new TStyle("histoStyle", "Plain Style(no colors/fill areas)");
    //   TStyle *histoStyle = gStyle;
 
    histoStyle->SetTextFont(textFont);
    histoStyle->SetPadColor(0);
    histoStyle->SetCanvasColor(0);
    histoStyle->SetTitleColor(0);
    histoStyle->SetStatColor(0);
 
    histoStyle->SetOptTitle(0);   // without main up title
    histoStyle->SetOptStat(1000111010);
    //   The parameter mode can be = ksiourmen  (default = 000001111)
    //   k = 1;  kurtosis printed
    //   k = 2;  kurtosis and kurtosis error printed
    //   s = 1;  skewness printed
    //   s = 2;  skewness and skewness error printed
    //   i = 1;  integral of bins printed
    //   o = 1;  number of overflows printed
    //   u = 1;  number of underflows printed
    //   r = 1;  rms printed
    //   r = 2;  rms and rms error printed
    //   m = 1;  mean value printed
    //   m = 2;  mean and mean error values printed
    //   e = 1;  number of entries printed
    //   n = 1;  name of histogram is printed
 
    histoStyle->SetOptFit(10001);
    //   The parameter mode can be = pcev  (default = 0111)
    //   p = 1;  print Probability
    //   c = 1;  print Chisquare/Number of degress of freedom
    //   e = 1;  print errors (if e=1, v must be 1)
    //   v = 1;  print name/values of parameters
 
    histoStyle->SetStatW(0.175);
    histoStyle->SetStatH(0.02);
    histoStyle->SetStatX(0.95);
    histoStyle->SetStatY(0.97);
    histoStyle->SetStatFontSize(0.05);
 
    histoStyle->SetStatFont(textFont);
    histoStyle->SetPadBorderMode(0);
    histoStyle->SetFrameBorderMode(0);
 
    histoStyle->SetPalette(1);   // for 2D
 
    TStyle* profStyle = new TStyle("profStyle", "Plain Style(no colors/fill areas)");
    //   TStyle *profStyle = gStyle;
 
    profStyle->SetTextFont(textFont);
    profStyle->SetPadColor(0);
    profStyle->SetCanvasColor(0);
    profStyle->SetTitleColor(0);
    profStyle->SetStatColor(0);
 
    profStyle->SetOptTitle(0);   // without main up title
    profStyle->SetOptStat(0);
    profStyle->SetOptFit(0);
    profStyle->SetGridStyle(2);
 
    profStyle->SetPadBorderMode(0);
    profStyle->SetFrameBorderMode(0);
 
    // ============================ Open file ============================
 
    const TString fileName = "L1_histo.root";
 
    TFile* fileIn = new TFile(fileName.Data(), "read");
 
    FILE *ress, *pulls;
    ress = fopen("resolution.txt", "w");
    pulls = fopen("pull.txt", "w");
 
    // ------------ make histos -----------
 
    // ------------ tracks -----------
 
    TString dirFitName = "/L1/Fit";
    dirFitName = fileName + ":" + dirFitName;
    TDirectory* dirFit = fileIn->Get(dirFitName);
 
    const int nHisto = 40;
 
    struct THistoData
    {
        char* name;    // one, which used in input file
        char* title;   // for X-axis
    };
    const THistoData histoData[nHisto] = {
        {"fst_x", "Residual (x^{reco}-x^{mc}) [#mum]"} {"fst_y", "Residual (y^{reco}-y^{mc}) [#mum]"} {
            "fst_tx", "Residual (t_{x}^{reco}-t_{x}^{mc}) [10^{-3}]"} {
            "fst_ty", "Residual (t_{y}^{reco}-t_{y}^{mc}) [10^{-3}]"} {"fst_P", "Resolution (p^{reco}-p^{mc})/p^{mc}"} {
            "fst_px", "Pull x"} {"fst_py", "Pull y"} {"fst_ptx", "Pull t_{x}"} {"fst_pty", "Pull t_{y}"} {"fst_pQP",
                                                                                                          "Pull q/p"}
 
        {"lst_pQP", "Pull q/p"} {"lst_ptx", "Pull t_{x}"} {"lst_pty", "Pull t_{y}"} {"lst_px", "Pull x"} {"lst_py",
                                                                                                          "Pull y"} {
            "lst_P", "Resolution (p^{reco}-p^{mc})/p^{mc}"} {"lst_tx", "Residual (t_{x}^{reco}-t_{x}^{mc}) [mrad]"} {
            "lst_ty", "Residual (t_{y}^{reco}-t_{y}^{mc}) [mrad]"} {"lst_x", "Residual (x^{reco}-x^{mc}) [#mum]"} {
            "lst_y", "Residual (y^{reco}-y^{mc}) [#mum]"} {"svrt_pQP", "Pull q/p"} {"svrt_ptx", "Pull t_{x}"} {
            "svrt_pty", "Pull t_{y}"} {"svrt_px", "Pull x"} {"svrt_py", "Pull y"} {
            "svrt_P", "Resolution (p^{reco}-p^{mc})/p^{mc}"} {"svrt_tx", "Residual (t_{x}^{reco}-t_{x}^{mc}) [mrad]"} {
            "svrt_ty", "Residual (t_{y}^{reco}-t_{y}^{mc}) [mrad]"} {"svrt_x", "Residual (x^{reco}-x^{mc}) [#mum]"} {
            "svrt_y", "Residual (y^{reco}-y^{mc}) [#mum]"} {"pvrt_pQP", "Pull q/p"} {"pvrt_ptx", "Pull t_{x}"} {
            "pvrt_pty", "Pull t_{y}"} {"pvrt_px", "Pull x"} {"pvrt_py", "Pull y"} {
            "pvrt_P", "Resolution (p^{reco}-p^{mc})/p^{mc}"} {"pvrt_tx", "Residual (t_{x}^{reco}-t_{x}^{mc}) [mrad]"} {
            "pvrt_ty", "Residual (t_{y}^{reco}-t_{y}^{mc}) [mrad]"} {"pvrt_x", "Residual (x^{reco}-x^{mc}) [#mum]"} {
            "pvrt_y", "Residual (y^{reco}-y^{mc}) [#mum]"}};
 
    TH1F* histo[nHisto];
 
    float charat[nHisto][2];   // 0 - sigma, 1 - RMS
    for (int i = 0; i < nHisto; i++) {
        histo[i] = (TH1F*)dirFit->Get(histoData[i].name);
        makeUpHisto(histo[i], histoData[i].title);
        FitHistoGaus(histo[i], &(charat[i][0]));
        charat[i][1] = histo[i]->GetRMS();
    }
 
    cout.setf(std::ios::scientific, std::ios::floatfield);
    cout.precision(2);
    cout << endl << "Name   " << "   \t" << " Sigma   " << "\t" << " RMS " << endl;
    for (int i = 0; i < nHisto; i++) {
        cout << histoData[i].name << "   \t" << charat[i][0] << "\t" << charat[i][1] << endl;
    }
    // ------------ hits, strips -----------
 
    TString dirInputName = "/L1/Input";
    dirInputName = fileName + ":" + dirInputName;
    TDirectory* dirInput = fileIn->Get(dirInputName);
 
    const int nHistoInput = 4;
 
    const THistoData histoDataInput[nHisto] = {{"Px", "Pull x"} {"Py", "Pull y"} {
        "x", "Residual (x^{hit}-x^{mc}) [#mum]"} {"y", "Residual (y^{hit}-y^{mc}) [#mum]"}};
 
    TH1F* histoInput[nHistoInput];
 
    for (int i = 0; i < nHistoInput; i++) {
        histoInput[i] = (TH1F*)dirInput->Get(histoDataInput[i].name);
        makeUpHisto(histoInput[i], histoDataInput[i].title);
        FitHistoGaus(histoInput[i]);
    }
 
    // ------------ other histo -----------
 
    TString dirOtherName = "/L1";
    dirOtherName = fileName + ":" + dirOtherName;
    TDirectory* dirOther = fileIn->Get(dirOtherName);
 
    const int nHistoOther = 15;
 
    const THistoData histoDataOther[nHisto] = {
        {"h_acc_MCmom", "Momentum of tracks"} {"h_reco_MCmom", "Momentum of tracks"} {
            "h_unreco_MCmom", "Momentum of tracks"} {"h_ghost_Rmom", "Momentum of tracks"} {
            "h_acc_prim_MCmom", "Momentum of tracks"} {"h_reco_prim_MCmom", "Momentum of tracks"} {
            "h_unreco_prim_MCmom", "Momentum of tracks"} {"h_ghost_Rmom", "Momentum of tracks"} {
            "h_acc_sec_MCmom", "Momentum of tracks"} {"h_reco_sec_MCmom", "Momentum of tracks"} {
            "h_unreco_sec_MCmom", "Momentum of tracks"} {"h_ghost_Rmom", "Momentum of tracks"}
        //  reg,acc,reco,ghost distribution always have to be first
        {"h_reco_prob", "Prob of reco track"} {"h_ghost_prob", "Prob of ghost track"} {"h_rest_prob",
                                                                                       "Prob of reco rest track"}};
 
    TH1F* histoOther[nHistoOther];
 
    for (int i = 0; i < 12; i++) {
        if (i % 4 != 2)
            histoOther[i] = (TH1F*)dirOther->Get(histoDataOther[i].name);
        else {
            histoOther[i] = new TH1F(histoDataOther[i].name, "Not reconstructed tracks", 100, 0.0, 5.0);
            for (j = 0; j < histoOther[i]->GetNbinsX(); j++) {
                histoOther[i]->SetBinContent(j,
                                             histoOther[i - 2]->GetBinContent(j) - histoOther[i - 1]->GetBinContent(j));
            }
        }
        makeUpHisto(histoOther[i], histoDataOther[i].title);
    }
 
    for (int i = 12; i < nHistoOther; i++) {
        histoOther[i] = (TH1F*)dirOther->Get(histoDataOther[i].name);
        makeUpHisto(histoOther[i], histoDataOther[i].title);
        FitHistoLine(histoOther[i]);
    }
 
    // ------------ other 2D histo -----------
 
    TString dirOther2DName = "/L1";
    dirOther2DName = fileName + ":" + dirOther2DName;
    TDirectory* dirOther2D = fileIn->Get(dirOther2DName);
 
    const int nHistoOther2D = 1;
 
    const THistoData histoDataOther2D[nHisto] = {
        {"h2_reco_nhits_vs_mom_prim", "NHits and momentum for primary"}   // {"h2_reco_prob_ndf",  "Prob and NDF" }
                                                                          // {"h2_reco_prob_mom",  "Prob and momentum" }
                                                                          // {"h2_reco_prob_teta",  "Prob and teta" }
                                                                          // {"h2_teta_prob",  "Prob" }
                                                                          // {"h2_teta_probI",  "ProbI" }
                                                                          // {"h2_teta_probO",  "ProbO" }
    };
 
    TH2F* histoOther2D[nHistoOther2D];
 
    for (int i = 0; i < nHistoOther2D; i++) {
        histoOther2D[i] = (TH2F*)dirOther2D->Get(histoDataOther2D[i].name);
        makeUpHisto(histoOther2D[i], histoDataOther2D[i].title);
    }
 
    // ------------ make profiles -----------
 
    TString dirName = "/L1";
    dirName = fileName + ":" + dirName;
    TDirectory* dir = fileIn->Get(dirName);
 
    const int nProf = 4;
    struct TProfData
    {
        char* name;     // one, which used in input file
        char* titleX;   // for X-axis
        char* titleY;   // for Y-axis
    };
    const TProfData profData[nProf] = {
        {"p_eff_all_vs_mom", "Momentum [GeV/c]", "Efficiency [%]"} {"p_eff_all_vs_nhits", "NMCHits", "Efficiency [%]"} {
            "p_eff_prim_vs_nhits", "NMCHits", "Efficiency [%]"} {"p_eff_sec_vs_nhits", "NMCHits", "Efficiency [%]"}};
 
    TProfile* profile[nProf];
 
    for (int i = 0; i < nProf; i++) {
        profile[i] = (TProfile*)dir->Get(profData[i].name);
        makeUpProfile(profile[i], profData[i].titleX, profData[i].titleY);
    }
 
    // ============================ Draw  ============================
    //   TCanvas *c1;
    //   c1 = new TCanvas("c1","c1",0,0,1200,800);
    //   c1 -> Divide(nHisto);
    //
    //   c1->UseCurrentStyle();
    //
    //   for (int i = 0; i < nHisto; i++){
    //     c1->cd(i+1);
    //     histo[i]->Draw();
    //   }
 
    //   c1->SaveAs("histo.png");
 
    system("mkdir L1_histo -p");
    chdir("L1_histo");
 
    histoStyle->cd();
    TCanvas *c2, *c3;
    c2 = new TCanvas("c2", "c2", 0, 0, 900, 600);
    c3 = new TCanvas("c3", "c3", 0, 0, 900, 600);
    c3->Divide(5, 2);
 
    int iPart = 1;
    for (int i = 0; i < 10; i++) {
        c3->cd(iPart++);
        histo[i]->Draw();
    }
    c3->SaveAs("fst.pdf");
 
    c2->cd();
    for (int i = 0; i < nHisto; i++) {
        histo[i]->Draw();
        TString name = TString(histoData[i].name) + ".pdf";
        c2->SaveAs(name);
    }
 
    for (int i = 0; i < nHistoInput; i++) {
        histoInput[i]->Draw();
        TString name = TString(histoDataInput[i].name) + ".pdf";
        c2->SaveAs(name);
    }
 
    // draw reg,acc,reco,ghost distribution
    histoStyle->SetOptStat(000000000);
    for (int j = 0, iH = 0; j < 3; j++) {
        TLegend* leg1 = new TLegend(0.2 + 0.3, .62, 0.55 + 0.3, .75);
        leg1->SetFillColor(0);
        leg1->SetLineWidth(0);
        leg1->SetBorderSize(0);
        leg1->SetTextFont(22);
        for (int i = 0; i < 4; i++, iH++) {
            switch (i) {
                case 0:
                    histoOther[iH]->SetLineColor(kBlue);
                    break;
                case 1:
                    histoOther[iH]->SetLineColor(kRed);
                    break;
                case 2:
                    histoOther[iH]->SetLineColor(kBlue);
                    histoOther[iH]->SetLineStyle(2);
                    break;
                case 3:
                    histoOther[iH]->SetLineColor(kBlack);
                    break;
            }
            if (i == 0)
                histoOther[iH]->Draw();
            else
                histoOther[iH]->Draw("same");
            leg1->AddEntry(histoOther[iH], histoOther[iH]->GetTitle());
        }
        leg1->Draw();
        switch (j) {
            case 0:
                c2->SaveAs("MultiplicityMom.pdf");
                break;
            case 1:
                c2->SaveAs("MultiplicityMomPrim.pdf");
                break;
            case 2:
                c2->SaveAs("MultiplicityMomSec.pdf");
                break;
        }
    }
    histoStyle->SetOptStat(1000111010);
 
    for (int i = 12; i < nHistoOther; i++) {
        histoOther[i]->Draw();
        TString name = TString(histoDataOther[i].name) + ".pdf";
        c2->SaveAs(name);
    }
 
    histoStyle->SetOptStat(1000000010);
    for (int i = 0; i < nHistoOther2D; i++) {
        histoOther2D[i]->Draw();
        TString name = TString(histoDataOther2D[i].name) + ".pdf";
        c2->SaveAs(name);
    }
    histoStyle->SetOptStat(1000111010);
 
    profStyle->cd();
    c2->cd();
    for (int i = 0; i < nProf; i++) {
        c2->SetGridx();
        c2->SetGridy();
        profile[i]->Draw();
        TString name = TString(profData[i].name) + ".pdf";
        c2->SaveAs(name);
    }
 
    chdir("..");
}