SiStripTools/bin/StatisticsPlots.cc

0001 #include "StatisticsPlots.h"
0002 #include <vector>
0003 #include <algorithm>
0004 #include <string>
0005 #include <map>
0006 #include <iostream>
0007 #include "TPad.h"
0008 #include "TH1D.h"
0009 #include "TFile.h"
0010 #include "TH2F.h"
0011 #include "TH1F.h"
0012 #include "TText.h"
0013 #include "DPGAnalysis/SiStripTools/interface/CommonAnalyzer.h"
0014 #include "TCanvas.h"
0015 #include "TGraphAsymmErrors.h"
0016
0017 void DeadTimeAPVCycle(TH1F* hist, const std::vector<int>& bins) {
0018   float ntotodd = 0;
0019   float ntoteven = 0;
0020   float nspecialodd = 0;
0021   float nspecialeven = 0;
0022   unsigned int nbintotodd = 0;
0023   unsigned int nbintoteven = 0;
0024   unsigned int nbinspecialodd = 0;
0025   unsigned int nbinspecialeven = 0;
0026
0027   for (int i = 1; i < hist->GetNbinsX() + 1; ++i) {
0028     bool isSpecial = false;
0029     for (unsigned int special = 0; special < bins.size(); ++special) {
0030       if (i == bins[special]) {
0031         isSpecial = true;
0032         break;
0033       }
0034     }
0035
0036     if (i % 2 == 0) {
0037       if (isSpecial) {
0038         ++nbinspecialeven;
0039         nspecialeven += hist->GetBinContent(i);
0040       } else {
0041         ++nbintoteven;
0042         ntoteven += hist->GetBinContent(i);
0043       }
0044     } else {
0045       if (isSpecial) {
0046         ++nbinspecialodd;
0047         nspecialodd += hist->GetBinContent(i);
0048       } else {
0049         ++nbintotodd;
0050         ntotodd += hist->GetBinContent(i);
0051       }
0052     }
0053   }
0054   std::cout << "Summary" << std::endl;
0055   std::cout << "Odd events " << ntotodd << " special " << nspecialodd << std::endl;
0056   std::cout << "Odd bins " << nbintotodd << " special " << nbinspecialodd << std::endl;
0057   std::cout << "Odd bins populations" << float(ntotodd) / nbintotodd << " special "
0058             << float(nspecialodd) / nbinspecialodd << std::endl;
0059   std::cout << "Even events " << ntoteven << " special " << nspecialeven << std::endl;
0060   std::cout << "Even bins " << nbintoteven << " special " << nbinspecialeven << std::endl;
0061   std::cout << "Even bins populations" << float(ntoteven) / nbintoteven << " special "
0062             << float(nspecialeven) / nbinspecialeven << std::endl;
0063
0064   float oddloss = nspecialodd - nbinspecialodd * ntotodd / nbintotodd;
0065   float evenloss = nspecialeven - nbinspecialeven * ntoteven / nbintoteven;
0066
0067   float fracloss = (oddloss + evenloss) / (ntotodd + ntoteven + nspecialodd + nspecialeven);
0068
0069   std::cout << "Loss " << oddloss << " " << evenloss << " " << fracloss << std::endl;
0070 }
0071
0072 TH1F* CombinedHisto(TFile& ff, const char* module, const char* histname) {
0073   CommonAnalyzer castat(&ff, "", module);
0074
0075   TH1F* result = nullptr;
0076
0077   std::vector<unsigned int> runs = castat.getRunList();
0078   std::sort(runs.begin(), runs.end());
0079
0080   {
0081     for (unsigned int i = 0; i < runs.size(); ++i) {
0082       char runlabel[100];
0083       sprintf(runlabel, "%d", runs[i]);
0084       char runpath[100];
0085       sprintf(runpath, "run_%d", runs[i]);
0086       castat.setPath(runpath);
0087
0088       TH1F* hist = (TH1F*)castat.getObject(histname);
0089
0090       if (hist) {
0091         if (result == nullptr) {
0092           result = new TH1F(*hist);
0093           result->Reset();
0094         }
0095         result->Add(hist);
0096         std::cout << hist->GetTitle() << " added: " << hist->GetEntries() << " " << result->GetEntries() << std::endl;
0097       }
0098     }
0099   }
0100
0101   return result;
0102 }
0103
0104 TH1F* TimeRatio(TFile& ff, const char* modulen, const char* moduled, const int irun, const int rebin) {
0105   CommonAnalyzer castatn(&ff, "", modulen);
0106   CommonAnalyzer castatd(&ff, "", moduled);
0107
0108   char runlabel[100];
0109   sprintf(runlabel, "%d", irun);
0110   char runpath[100];
0111   sprintf(runpath, "run_%d", irun);
0112   castatn.setPath(runpath);
0113   castatd.setPath(runpath);
0114
0115   TH1F* ratio = nullptr;
0116
0117   TH1F* orbitn = nullptr;
0118   if (orbitn == nullptr)
0119     orbitn = (TH1F*)castatn.getObject("orbit");
0120   TH1F* orbitd = nullptr;
0121   if (orbitd == nullptr)
0122     orbitd = (TH1F*)castatd.getObject("orbit");
0123   if (orbitn != nullptr && orbitd != nullptr) {
0124     orbitn->Rebin(rebin);
0125     orbitd->Rebin(rebin);
0126     ratio = new TH1F(*orbitd);
0127     ratio->Reset();
0128     ratio->Divide(orbitn, orbitd);
0129   }
0130   return ratio;
0131 }
0132
0133 TH1D* SummaryHisto(TFile& ff, const char* module) {
0134   CommonAnalyzer castat(&ff, "", module);
0135
0136   TH1D* nevents = new TH1D("nevents", "Number of events vs run", 10, 0., 10.);
0137   nevents->SetCanExtend(TH1::kXaxis);
0138
0139   std::vector<unsigned int> runs = castat.getRunList();
0140   std::sort(runs.begin(), runs.end());
0141
0142   {
0143     for (unsigned int i = 0; i < runs.size(); ++i) {
0144       char runlabel[100];
0145       sprintf(runlabel, "%d", runs[i]);
0146       char runpath[100];
0147       sprintf(runpath, "run_%d", runs[i]);
0148       castat.setPath(runpath);
0149
0150       TH1F* orbit = nullptr;
0151       if (orbit == nullptr)
0152         orbit = (TH1F*)castat.getObject("orbit");
0153       if (orbit) {
0154         // fill the summary
0155         nevents->Fill(runlabel, orbit->GetEntries());
0156       }
0157     }
0158   }
0159   return nevents;
0160 }
0161
0162 TH1D* SummaryHistoRatio(TFile& f1, const char* mod1, TFile& f2, const char* mod2, const char* hname) {
0163   TH1D* denom = SummaryHisto(f1, mod1);
0164   TH1D* numer = SummaryHisto(f2, mod2);
0165
0166   TH1D* ratio = (TH1D*)denom->Clone(hname);
0167   ratio->SetTitle("Fraction of events vs run");
0168   ratio->Sumw2();
0169   ratio->Reset();
0170   ratio->SetDirectory(nullptr);
0171   ratio->Divide(numer, denom, 1, 1, "b");
0172
0173   return ratio;
0174 }
0175
0176 TGraphAsymmErrors* SummaryHistoRatioGraph(
0177     TFile& f1, const char* mod1, TFile& f2, const char* mod2, const char* /* hname */) {
0178   TH1D* denom = SummaryHisto(f1, mod1);
0179   TH1D* numer = SummaryHisto(f2, mod2);
0180
0181   TGraphAsymmErrors* ratio = new TGraphAsymmErrors;
0182   ;
0183
0184   ratio->BayesDivide(numer, denom);
0185
0186   return ratio;
0187 }
0188
0189 TH2F* Combined2DHisto(TFile& ff, const char* module, const char* histname) {
0190   CommonAnalyzer castat(&ff, "", module);
0191
0192   TH2F* result = nullptr;
0193
0194   std::vector<unsigned int> runs = castat.getRunList();
0195   std::sort(runs.begin(), runs.end());
0196
0197   {
0198     for (unsigned int i = 0; i < runs.size(); ++i) {
0199       char runlabel[100];
0200       sprintf(runlabel, "%d", runs[i]);
0201       char runpath[100];
0202       sprintf(runpath, "run_%d", runs[i]);
0203       castat.setPath(runpath);
0204
0205       TH2F* hist = (TH2F*)castat.getObject(histname);
0206
0207       if (hist) {
0208         if (result == nullptr) {
0209           result = new TH2F(*hist);
0210           result->Reset();
0211         }
0212         result->Add(hist);
0213         std::cout << hist->GetTitle() << " added: " << hist->GetEntries() << " " << result->GetEntries() << std::endl;
0214       }
0215     }
0216   }
0217
0218   return result;
0219 }
0220
0221 void StatisticsPlots(const char* fullname,
0222                      const char* module,
0223                      const char* label,
0224                      const char* postfix,
0225                      const char* shortname,
0226                      const char* outtrunk) {
0227   char modfull[300];
0228   sprintf(modfull, "%s%s", module, postfix);
0229   char labfull[300];
0230   sprintf(labfull, "%s%s", label, postfix);
0231
0232   char dirname[400];
0233   sprintf(dirname, "%s", shortname);
0234
0235   //  char fullname[300];
0236   //  if(strlen(family)==0) {  sprintf(fullname,"rootfiles/Tracking_PFG_%s.root",filename);}
0237   //  else {  sprintf(fullname,"rootfiles/%s.root",dirname); }
0238
0239   TFile ff(fullname);
0240
0241   // Colliding events
0242
0243   /*
0244   char modname[300];
0245   sprintf(modname,"trackcount%s",postfix);
0246
0247   CommonAnalyzer castatall(&ff,"",modname);
0248
0249   TH1F* ntrk = (TH1F*)castatall.getObject("ntrk");
0250   if (ntrk) {
0251     std::cout << " All runs " << ntrk->GetEntries() << std::endl;
0252     delete ntrk;
0253   }
0254   */
0255
0256   //  sprintf(modname,"eventtimedistribution%s",postfix);
0257
0258   CommonAnalyzer castat(&ff, "", modfull);
0259
0260   // Summary histograms
0261
0262   TH1D* nevents = new TH1D("nevents", "Number of events vs run", 10, 0., 10.);
0263   nevents->SetCanExtend(TH1::kXaxis);
0264
0265   std::vector<unsigned int> runs = castat.getRunList();
0266   std::sort(runs.begin(), runs.end());
0267
0268   {
0269     std::cout << "Collision Events" << std::endl;
0270
0271     for (unsigned int i = 0; i < runs.size(); ++i) {
0272       char runlabel[100];
0273       sprintf(runlabel, "%d", runs[i]);
0274       char runpath[100];
0275       sprintf(runpath, "run_%d", runs[i]);
0276       castat.setPath(runpath);
0277
0278       TH1* orbit = nullptr;
0279       TH2F* orbitvsbx = (TH2F*)castat.getObject("orbitvsbxincycle");
0280       if (orbitvsbx == nullptr)
0281         orbitvsbx = (TH2F*)castat.getObject("orbitvsbx");
0282       if (orbitvsbx) {
0283         std::cout << runpath << " " << orbitvsbx->GetEntries() << std::endl;
0284         // prepare plot
0285         new TCanvas;
0286         if (orbitvsbx->GetEntries() > 0) {
0287           orbit = orbitvsbx->ProjectionY();
0288         }
0289         delete orbitvsbx;
0290       }
0291       if (orbit == nullptr)
0292         orbit = (TH1F*)castat.getObject("orbit");
0293       if (orbit) {
0294         orbit->GetXaxis()->SetTitle("orbit");
0295         orbit->SetTitle(runpath);
0296         //normalize to get the rate
0297         orbit->Scale(11223, "width");
0298         orbit->GetYaxis()->SetTitle("Rate (Hz)");
0299         // fill the summary
0300         nevents->Fill(runlabel, orbit->GetEntries());
0301
0302         //
0303         orbit->Draw();
0304         std::string plotfilename;
0305         plotfilename += outtrunk;
0306         plotfilename += dirname;
0307         plotfilename += "/orbit_";
0308         plotfilename += labfull;
0309         plotfilename += "_";
0310         plotfilename += dirname;
0311         plotfilename += "_";
0312         plotfilename += runpath;
0313         plotfilename += ".gif";
0314         gPad->Print(plotfilename.c_str());
0315         delete orbit;
0316       }
0317
0318       TH1F* dbx = (TH1F*)castat.getObject("dbx");
0319       if (dbx) {
0320         // prepare plot
0321         if (dbx->GetEntries() > 0) {
0322           dbx->Draw();
0323           gPad->SetLogy(1);
0324           std::string plotfilename;
0325           plotfilename += outtrunk;
0326           plotfilename += dirname;
0327           plotfilename += "/dbx_";
0328           plotfilename += labfull;
0329           plotfilename += "_";
0330           plotfilename += dirname;
0331           plotfilename += "_";
0332           plotfilename += runpath;
0333           plotfilename += ".gif";
0334           gPad->Print(plotfilename.c_str());
0335           delete dbx;
0336         }
0337         gPad->SetLogy(0);
0338       }
0339       TH1F* bx = (TH1F*)castat.getObject("bx");
0340       if (bx) {
0341         // prepare plot
0342         if (bx->GetEntries() > 0) {
0343           bx->SetLineColor(kRed);
0344           bx->Draw();
0345           gPad->SetLogy(1);
0346           std::string plotfilename;
0347           plotfilename += outtrunk;
0348           plotfilename += dirname;
0349           plotfilename += "/bx_";
0350           plotfilename += labfull;
0351           plotfilename += "_";
0352           plotfilename += dirname;
0353           plotfilename += "_";
0354           plotfilename += runpath;
0355           plotfilename += ".gif";
0356           gPad->Print(plotfilename.c_str());
0357           delete bx;
0358         }
0359         gPad->SetLogy(0);
0360       }
0361     }
0362   }
0363
0364   if (!runs.empty()) {
0365     std::string plotfilename;
0366     plotfilename = outtrunk;
0367     plotfilename += dirname;
0368     plotfilename += "/nevents_";
0369     plotfilename += labfull;
0370     plotfilename += "_";
0371     plotfilename += dirname;
0372     plotfilename += ".gif";
0373
0374     TCanvas* cwide = new TCanvas(plotfilename.c_str(), plotfilename.c_str(), 1500, 500);
0375     nevents->Draw();
0376
0377     char nentries[100];
0378     sprintf(nentries, "%.0f", nevents->GetSumOfWeights());
0379     TText ttlabel;
0380     ttlabel.DrawTextNDC(.8, .7, nentries);
0381
0382     std::cout << nentries << std::endl;
0383
0384     gPad->Print(plotfilename.c_str());
0385     delete cwide;
0386   }
0387   delete nevents;
0388
0389   ff.Close();
0390 }