Skip to content
Snippets Groups Projects
Draw_L1_histo.C 14.43 KiB
/* Copyright (C) 2010-2013 Frankfurt Institute for Advanced Studies, Goethe-Universität Frankfurt, Frankfurt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Igor Kulakov [committer] */

#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("..");
}