StatisticsPlots.cc

Go to the documentation of this file.

#include "StatisticsPlots.h"
#include <vector>
#include <algorithm>
#include <string>
#include <map>
#include <iostream>
#include "TPad.h"
#include "TH1D.h"
#include "TFile.h"
#include "TH2F.h"
#include "TH1F.h"
#include "TText.h"
#include "DPGAnalysis/SiStripTools/interface/CommonAnalyzer.h"
#include "TCanvas.h"
#include "TGraphAsymmErrors.h"

void DeadTimeAPVCycle(TH1F* hist, const std::vector<int>& bins) {

  float ntotodd = 0;
  float ntoteven = 0;
  float nspecialodd = 0;
  float nspecialeven = 0;
  unsigned int nbintotodd = 0;
  unsigned int nbintoteven = 0;
  unsigned int nbinspecialodd = 0;
  unsigned int nbinspecialeven = 0;

  for(int i = 1; i < hist->GetNbinsX()+1; ++i) {

    bool isSpecial = false;
    for(unsigned int special = 0; special < bins.size() ; ++special) {
      if(i==bins[special]) {
    isSpecial=true;
    break;
      }
    }

    if(i%2==0) {
      if(isSpecial) {
    ++nbinspecialeven;
    nspecialeven += hist->GetBinContent(i);
      }
      else {
    ++nbintoteven;
    ntoteven += hist->GetBinContent(i);
      }
    }
    else {
      if(isSpecial) {
    ++nbinspecialodd;
    nspecialodd += hist->GetBinContent(i);
      }
      else {
    ++nbintotodd;
    ntotodd += hist->GetBinContent(i);
      }
    }
  }
  std::cout << "Summary" << std::endl;
  std::cout << "Odd events " << ntotodd << " special " << nspecialodd << std::endl;
  std::cout << "Odd bins " << nbintotodd << " special " << nbinspecialodd << std::endl;
  std::cout << "Odd bins populations" << float(ntotodd)/nbintotodd << " special " << float(nspecialodd)/nbinspecialodd << std::endl;
  std::cout << "Even events " << ntoteven << " special " << nspecialeven << std::endl;
  std::cout << "Even bins " << nbintoteven << " special " << nbinspecialeven << std::endl;
  std::cout << "Even bins populations" << float(ntoteven)/nbintoteven << " special " << float(nspecialeven)/nbinspecialeven << std::endl;

  float oddloss = nspecialodd -nbinspecialodd*ntotodd/nbintotodd;
  float evenloss = nspecialeven -nbinspecialeven*ntoteven/nbintoteven;

  float fracloss = (oddloss+evenloss)/(ntotodd + ntoteven + nspecialodd + nspecialeven);

  std::cout << "Loss " << oddloss << " " << evenloss << " " << fracloss << std::endl;

}

TH1F* CombinedHisto(TFile& ff, const char* module, const char* histname) {
  
  CommonAnalyzer castat(&ff,"",module);
  
  TH1F* result = nullptr;
  
  std::vector<unsigned int> runs = castat.getRunList();
  std::sort(runs.begin(),runs.end());
  
  {
    for(unsigned int i=0;i<runs.size();++i) {
      
      char runlabel[100];
      sprintf(runlabel,"%d",runs[i]);
      char runpath[100];
      sprintf(runpath,"run_%d",runs[i]);
      castat.setPath(runpath);
      
      TH1F* hist = (TH1F*)castat.getObject(histname);
      
      if(hist) {
    if(result==nullptr) {
      result = new TH1F(*hist);
      result->Reset();
    }
    result->Add(hist);
    std::cout << hist->GetTitle() << " added: " << hist->GetEntries() << " " << result->GetEntries() << std::endl;
      }
      
    }
  }
  
  return result;
  
}

TH1F* TimeRatio(TFile& ff, const char* modulen, const char* moduled, const int irun, const int rebin) {
  
  CommonAnalyzer castatn(&ff,"",modulen);
  CommonAnalyzer castatd(&ff,"",moduled);
  
  
  char runlabel[100];
  sprintf(runlabel,"%d",irun);
  char runpath[100];
  sprintf(runpath,"run_%d",irun);
  castatn.setPath(runpath);
  castatd.setPath(runpath);
  
  TH1F* ratio = nullptr;
  
  TH1F* orbitn=nullptr;
  if(orbitn==nullptr) orbitn = (TH1F*)castatn.getObject("orbit");
  TH1F* orbitd=nullptr;
  if(orbitd==nullptr) orbitd = (TH1F*)castatd.getObject("orbit");
  if(orbitn != nullptr && orbitd != nullptr) {
    orbitn->Rebin(rebin);
    orbitd->Rebin(rebin);
    ratio = new TH1F(*orbitd);
    ratio->Reset();
    ratio->Divide(orbitn,orbitd);
  }
  return ratio;
  
}

TH1D* SummaryHisto(TFile& ff, const char* module) {
  
  CommonAnalyzer castat(&ff,"",module);
  
  TH1D* nevents = new TH1D("nevents","Number of events vs run",10,0.,10.);
  nevents->SetCanExtend(TH1::kXaxis);
  
  std::vector<unsigned int> runs = castat.getRunList();
  std::sort(runs.begin(),runs.end());
  
  {
    for(unsigned int i=0;i<runs.size();++i) {
      
      char runlabel[100];
      sprintf(runlabel,"%d",runs[i]);
      char runpath[100];
      sprintf(runpath,"run_%d",runs[i]);
      castat.setPath(runpath);
      
      
      TH1F* orbit=nullptr;
      if(orbit==nullptr) orbit = (TH1F*)castat.getObject("orbit");
      if(orbit) {
    // fill the summary
    nevents->Fill(runlabel,orbit->GetEntries());
      }
    }
  } 
  return nevents;
  
}

TH1D* SummaryHistoRatio(TFile& f1, const char* mod1, TFile& f2, const char* mod2, const char* hname) {

  TH1D* denom = SummaryHisto(f1,mod1);
  TH1D* numer = SummaryHisto(f2,mod2);

  TH1D* ratio = (TH1D*)denom->Clone(hname);
  ratio->SetTitle("Fraction of events vs run");
  ratio->Sumw2();
  ratio->Reset();
  ratio->SetDirectory(nullptr);
  ratio->Divide(numer,denom,1,1,"b");

  return ratio;

}

TGraphAsymmErrors* SummaryHistoRatioGraph(TFile& f1, const char* mod1, TFile& f2, const char* mod2, const char* /* hname */) {

  TH1D* denom = SummaryHisto(f1,mod1);
  TH1D* numer = SummaryHisto(f2,mod2);

  TGraphAsymmErrors* ratio = new TGraphAsymmErrors;;

  ratio->BayesDivide(numer,denom);

  return ratio;

}

TH2F* Combined2DHisto(TFile& ff, const char* module, const char* histname) {
  
  CommonAnalyzer castat(&ff,"",module);
  
  TH2F* result = nullptr;
  
  std::vector<unsigned int> runs = castat.getRunList();
  std::sort(runs.begin(),runs.end());
  
  {
    for(unsigned int i=0;i<runs.size();++i) {
      
      char runlabel[100];
      sprintf(runlabel,"%d",runs[i]);
      char runpath[100];
      sprintf(runpath,"run_%d",runs[i]);
      castat.setPath(runpath);
      
      TH2F* hist = (TH2F*)castat.getObject(histname);
      
      if(hist) {
    if(result==nullptr) {
      result = new TH2F(*hist);
      result->Reset();
    }
    result->Add(hist);
    std::cout << hist->GetTitle() << " added: " << hist->GetEntries() << " " << result->GetEntries() << std::endl;
      }
      
    }
  }
  
  return result;
  
}



void StatisticsPlots(const char* fullname, const char* module, const char* label, const char* postfix, const char* shortname,
             const char* outtrunk) {
  
  char modfull[300];
  sprintf(modfull,"%s%s",module,postfix);
  char labfull[300];
  sprintf(labfull,"%s%s",label,postfix);

  char dirname[400];
  sprintf(dirname,"%s",shortname);

  //  char fullname[300];
  //  if(strlen(family)==0) {  sprintf(fullname,"rootfiles/Tracking_PFG_%s.root",filename);}
  //  else {  sprintf(fullname,"rootfiles/%s.root",dirname); }

  TFile ff(fullname);

  // Colliding events

  /*
  char modname[300];
  sprintf(modname,"trackcount%s",postfix);

  CommonAnalyzer castatall(&ff,"",modname);

  TH1F* ntrk = (TH1F*)castatall.getObject("ntrk");
  if (ntrk) {
    std::cout << " All runs " << ntrk->GetEntries() << std::endl;
    delete ntrk;
  }
  */
  
  //  sprintf(modname,"eventtimedistribution%s",postfix);

  CommonAnalyzer castat(&ff,"",modfull);

  // Summary histograms

  TH1D* nevents = new TH1D("nevents","Number of events vs run",10,0.,10.);
  nevents->SetCanExtend(TH1::kXaxis);


  std::vector<unsigned int> runs = castat.getRunList();
  std::sort(runs.begin(),runs.end());

  {
    
    std::cout << "Collision Events" << std::endl;
    
    for(unsigned int i=0;i<runs.size();++i) {
      
      char runlabel[100];
      sprintf(runlabel,"%d",runs[i]);
      char runpath[100];
      sprintf(runpath,"run_%d",runs[i]);
      castat.setPath(runpath);

      TH1* orbit=nullptr;
      TH2F* orbitvsbx = (TH2F*)castat.getObject("orbitvsbxincycle");
      if(orbitvsbx==nullptr) orbitvsbx = (TH2F*)castat.getObject("orbitvsbx");
      if (orbitvsbx) {
    std::cout << runpath << " " << orbitvsbx->GetEntries() << std::endl;
    // prepare plot
    new TCanvas;
    if(orbitvsbx->GetEntries()>0) {
      orbit = orbitvsbx->ProjectionY();
    }
    delete orbitvsbx;
      }
      if(orbit==nullptr) orbit = (TH1F*)castat.getObject("orbit");
      if(orbit) {
    orbit->GetXaxis()->SetTitle("orbit");
    orbit->SetTitle(runpath);
    //normalize to get the rate
    orbit->Scale(11223,"width");
    orbit->GetYaxis()->SetTitle("Rate (Hz)");
    // fill the summary
    nevents->Fill(runlabel,orbit->GetEntries());

    //
    orbit->Draw();
    std::string plotfilename;
    plotfilename += outtrunk;
    plotfilename += dirname;
    plotfilename += "/orbit_";
    plotfilename += labfull;
    plotfilename += "_";
    plotfilename += dirname;
    plotfilename += "_";
    plotfilename += runpath;
    plotfilename += ".gif";
    gPad->Print(plotfilename.c_str());
    delete orbit;
      }

      TH1F* dbx = (TH1F*)castat.getObject("dbx");
      if (dbx) {
    // prepare plot
    if(dbx->GetEntries()>0) {
      dbx->Draw();
      gPad->SetLogy(1);
      std::string plotfilename;
      plotfilename += outtrunk;
      plotfilename += dirname;
      plotfilename += "/dbx_";
      plotfilename += labfull;
      plotfilename += "_";
      plotfilename += dirname;
      plotfilename += "_";
      plotfilename += runpath;
      plotfilename += ".gif";
      gPad->Print(plotfilename.c_str());
      delete dbx;
    }
    gPad->SetLogy(0);
      }
      TH1F* bx = (TH1F*)castat.getObject("bx");
      if (bx) {
    // prepare plot
    if(bx->GetEntries()>0) {
      bx->SetLineColor(kRed);
      bx->Draw();
      gPad->SetLogy(1);
      std::string plotfilename;
      plotfilename += outtrunk;
      plotfilename += dirname;
      plotfilename += "/bx_";
      plotfilename += labfull;
      plotfilename += "_";
      plotfilename += dirname;
      plotfilename += "_";
      plotfilename += runpath;
      plotfilename += ".gif";
      gPad->Print(plotfilename.c_str());
      delete bx;
    }
    gPad->SetLogy(0);
      }
    }
  }

  if(!runs.empty()) {
    std::string plotfilename;
    plotfilename = outtrunk;
    plotfilename += dirname;
    plotfilename += "/nevents_";
    plotfilename += labfull;
    plotfilename += "_";
    plotfilename += dirname;
    plotfilename += ".gif";

    TCanvas * cwide = new TCanvas(plotfilename.c_str(),plotfilename.c_str(),1500,500);
    nevents->Draw();
    
    char nentries[100];
    sprintf(nentries,"%.0f",nevents->GetSumOfWeights());
    TText ttlabel;
    ttlabel.DrawTextNDC(.8,.7,nentries);

    std::cout << nentries << std::endl;

    gPad->Print(plotfilename.c_str());
    delete cwide;
  }
  delete nevents;

  ff.Close();

}