File indexing completed on 2024-04-06 12:06:33
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
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* ) {
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
0236
0237
0238
0239 TFile ff(fullname);
0240
0241
0242
0243
0244
0245
0246
0247
0248
0249
0250
0251
0252
0253
0254
0255
0256
0257
0258 CommonAnalyzer castat(&ff, "", modfull);
0259
0260
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
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
0297 orbit->Scale(11223, "width");
0298 orbit->GetYaxis()->SetTitle("Rate (Hz)");
0299
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
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
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 }