SiStripTools/bin/OOTMultiplicityPlotMacros.cc

0001 #include "OOTMultiplicityPlotMacros.h"
0002 #include "DPGAnalysis/SiStripTools/interface/CommonAnalyzer.h"
0003 #include "TFile.h"
0004 #include "TCanvas.h"
0005 #include "TH1F.h"
0006 #include "TF1.h"
0007 #include "TProfile.h"
0008 #include <iostream>
0009 #include <algorithm>
0010 #include <cmath>
0011
0012 OOTSummary* ComputeOOTFractionvsFill(TFile* ff,
0013                                      const char* itmodule,
0014                                      const char* ootmodule,
0015                                      const char* etmodule,
0016                                      const char* hname,
0017                                      OOTSummary* ootsumm) {
0018   if (ootsumm == nullptr) {
0019     ootsumm = new OOTSummary;
0020   }
0021
0022   if (ff) {
0023     CommonAnalyzer ca(ff, "", itmodule);
0024     std::vector<unsigned int> runs = ca.getFillList();
0025     std::sort(runs.begin(), runs.end());
0026     for (unsigned int i = 0; i < runs.size(); ++i) {
0027       char runlabel[100];
0028       sprintf(runlabel, "%d", runs[i]);
0029
0030       OOTResult* res = ComputeOOTFraction(ff, itmodule, ootmodule, etmodule, runs[i], hname, true);
0031
0032       if (res->ngoodbx != res->hratio->GetEntries())
0033         std::cout << "Inconsistency in number of good bx" << std::endl;
0034       ootsumm->hngoodbx->Fill(runlabel, res->ngoodbx);
0035       int ibin = ootsumm->hootfracsum->Fill(runlabel, res->ootfracsum);
0036       ootsumm->hootfracsum->SetBinError(ibin, res->ootfracsumerr);
0037       int bin = ootsumm->hootfrac->Fill(runlabel, res->ootfrac);
0038       ootsumm->hootfrac->SetBinError(bin, res->ootfracerr);
0039       delete res;
0040     }
0041   } else {
0042     std::cout << "File is not ok" << std::endl;
0043   }
0044
0045   return ootsumm;
0046 }
0047 OOTSummary* ComputeOOTFractionvsRun(TFile* ff,
0048                                     const char* itmodule,
0049                                     const char* ootmodule,
0050                                     const char* etmodule,
0051                                     const char* hname,
0052                                     OOTSummary* ootsumm) {
0053   if (ootsumm == nullptr) {
0054     ootsumm = new OOTSummary;
0055   }
0056
0057   if (ff) {
0058     CommonAnalyzer ca(ff, "", itmodule);
0059     std::vector<unsigned int> runs = ca.getRunList();
0060     std::sort(runs.begin(), runs.end());
0061     for (unsigned int i = 0; i < runs.size(); ++i) {
0062       char runlabel[100];
0063       sprintf(runlabel, "%d", runs[i]);
0064
0065       OOTResult* res = ComputeOOTFraction(ff, itmodule, ootmodule, etmodule, runs[i], hname);
0066
0067       if (res->ngoodbx != res->hratio->GetEntries())
0068         std::cout << "Inconsistency in number of good bx" << std::endl;
0069       ootsumm->hngoodbx->Fill(runlabel, res->ngoodbx);
0070       int ibin = ootsumm->hootfracsum->Fill(runlabel, res->ootfracsum);
0071       ootsumm->hootfracsum->SetBinError(ibin, res->ootfracsumerr);
0072       int bin = ootsumm->hootfrac->Fill(runlabel, res->ootfrac);
0073       ootsumm->hootfrac->SetBinError(bin, res->ootfracerr);
0074       delete res;
0075     }
0076   } else {
0077     std::cout << "File is not ok" << std::endl;
0078   }
0079
0080   return ootsumm;
0081 }
0082
0083 OOTResult* ComputeOOTFraction(TFile* ff,
0084                               const char* itmodule,
0085                               const char* ootmodule,
0086                               const char* etmodule,
0087                               const int run,
0088                               const char* hname,
0089                               const bool& perFill) {
0090   if (perFill) {
0091     std::cout << "Processing fill " << run << std::endl;
0092   } else {
0093     std::cout << "Processing run " << run << std::endl;
0094   }
0095
0096   char itpath[100];
0097   char ootpath[100];
0098   char etpath[100];
0099   if (perFill) {
0100     sprintf(itpath, "%s/fill_%d", itmodule, run);
0101     sprintf(ootpath, "%s/fill_%d", ootmodule, run);
0102     sprintf(etpath, "%s/fill_%d", etmodule, run);
0103   } else {
0104     sprintf(itpath, "%s/run_%d", itmodule, run);
0105     sprintf(ootpath, "%s/run_%d", ootmodule, run);
0106     sprintf(etpath, "%s/run_%d", etmodule, run);
0107   }
0108
0109   OOTResult* res = new OOTResult;
0110
0111   std::vector<int> filledbx = FillingSchemeFromProfile(ff, itpath, hname);
0112   res->nfilledbx = filledbx.size();
0113
0114   if (!perFill) {
0115     std::vector<int> filledbxtest = FillingScheme(ff, etpath);
0116     if (filledbx.size() != filledbxtest.size())
0117       std::cout << "Inconsistency in number of filled BX " << run << " " << filledbx.size() << " "
0118                 << filledbxtest.size() << std::endl;
0119   }
0120
0121   TProfile* itmult = nullptr;
0122   TProfile* ootmult = nullptr;
0123   res->hratio = new TH1F("ratio", "ratio", 200, 0., 2.);
0124
0125   float rclzb = 0;
0126   float rclrandom = 0;
0127   float errclzb = 0;
0128   float errclrandom = 0;
0129   float nzb = 0;
0130   float nrandom = 0;
0131
0132   if (ff) {
0133     if (ff->cd(itpath)) {
0134       itmult = (TProfile*)gDirectory->Get(hname);
0135     } else {
0136       std::cout << "In time path is not ok" << std::endl;
0137     }
0138     if (ff->cd(ootpath)) {
0139       ootmult = (TProfile*)gDirectory->Get(hname);
0140     } else {
0141       std::cout << "out of time path is not ok" << std::endl;
0142     }
0143     if (itmult && ootmult) {
0144       ootmult->SetLineColor(kRed);
0145       ootmult->SetMarkerColor(kRed);
0146       //      ootmult->Draw();
0147       //      itmult->Draw("same");
0148       for (std::vector<int>::const_iterator fbx = filledbx.begin(); fbx != filledbx.end(); ++fbx) {
0149         nzb += itmult->GetBinEntries(*fbx);
0150         nrandom += ootmult->GetBinEntries(*fbx + 1);
0151       }
0152       for (std::vector<int>::const_iterator fbx = filledbx.begin(); fbx != filledbx.end(); ++fbx) {
0153         if (nzb > 0 && nrandom > 0) {
0154           rclzb += (itmult->GetBinContent(*fbx) * itmult->GetBinEntries(*fbx)) / nzb;
0155           errclzb += (itmult->GetBinError(*fbx) * itmult->GetBinEntries(*fbx)) *
0156                      (itmult->GetBinError(*fbx) * itmult->GetBinEntries(*fbx)) / (nzb * nzb);
0157           rclrandom += (ootmult->GetBinContent(*fbx + 1) * ootmult->GetBinEntries(*fbx + 1)) / nrandom;
0158           errclrandom += (ootmult->GetBinError(*fbx + 1) * ootmult->GetBinEntries(*fbx + 1)) *
0159                          (ootmult->GetBinError(*fbx + 1) * ootmult->GetBinEntries(*fbx + 1)) / (nrandom * nrandom);
0160         }
0161         if (itmult->GetBinContent(*fbx) == 0) {
0162           std::cout << "No cluster in filled BX! " << *fbx << std::endl;
0163         } else if (ootmult->GetBinEntries(*fbx + 1) ==
0164                    0) { /* std::cout << "No entry in OOT BX " << *fbx+1 << std::endl; */
0165         } else {
0166           float rat = ootmult->GetBinContent(*fbx + 1) / itmult->GetBinContent(*fbx);
0167           res->hratio->Fill(rat);
0168         }
0169       }
0170     } else {
0171       std::cout << "histograms not found" << std::endl;
0172     }
0173   } else {
0174     std::cout << "Input file pointer is not ok" << std::endl;
0175   }
0176
0177   res->ngoodbx = res->hratio->GetEntries();
0178
0179   if (nzb > 0 && nrandom > 0 && rclzb > 0 && rclrandom > 0) {
0180     res->ootfracsum = rclrandom / rclzb;
0181     res->ootfracsumerr = rclrandom / rclzb * sqrt(errclzb / (rclzb * rclzb) + errclrandom / (rclrandom * rclrandom));
0182   }
0183   if (res->ngoodbx) {
0184     res->hratio->Fit("gaus", "Q0", "", .01, 1.99);
0185     if (res->hratio->GetFunction("gaus")) {
0186       res->ootfrac = res->hratio->GetFunction("gaus")->GetParameter(1);
0187       res->ootfracerr = res->hratio->GetFunction("gaus")->GetParError(1);
0188     } else {
0189       std::cout << "Missing fitting function" << std::endl;
0190     }
0191   } else {
0192     std::cout << "No filled BX or strange filling scheme" << std::endl;
0193   }
0194
0195   return res;
0196 }
0197
0198 std::vector<int> FillingScheme(TFile* ff, const char* path, const float thr) {
0199   TH1F* bx = nullptr;
0200   std::vector<int> filledbx;
0201   if (ff) {
0202     if (ff->cd(path)) {
0203       bx = (TH1F*)gDirectory->Get("bx");
0204       if (bx) {
0205         //  bx->Draw();
0206         std::cout << "Number of entries " << bx->GetEntries() << " threshold " << thr / 3564. * bx->GetEntries()
0207                   << std::endl;
0208         for (int i = 1; i < bx->GetNbinsX() + 1; ++i) {
0209           if (bx->GetBinContent(i) > thr / 3564. * bx->GetEntries()) {
0210             if (!filledbx.empty() && i == filledbx[filledbx.size() - 1] + 1) {
0211               std::cout << "This is not a 50ns run ! " << std::endl;
0212               filledbx.clear();
0213               return filledbx;
0214             }
0215             filledbx.push_back(i);
0216           }
0217         }
0218       } else {
0219         std::cout << "Histogram not found" << std::endl;
0220       }
0221     } else {
0222       std::cout << "module path is not ok" << std::endl;
0223     }
0224   } else {
0225     std::cout << "Input file pointer is not ok" << std::endl;
0226   }
0227
0228   //  std::cout << filledbx.size() << " filled bunch crossings" << std::endl;
0229   //  for(std::vector<int>::const_iterator fbx=filledbx.begin();fbx!=filledbx.end();++fbx) { std::cout << *fbx << std::endl;}
0230   return filledbx;
0231 }
0232 std::vector<int> FillingSchemeFromProfile(TFile* ff, const char* path, const char* hname, const float thr) {
0233   TProfile* bx = nullptr;
0234   std::vector<int> filledbx;
0235   if (ff) {
0236     if (ff->cd(path)) {
0237       bx = (TProfile*)gDirectory->Get(hname);
0238       if (bx) {
0239         //  bx->Draw();
0240         std::cout << "Number of entries " << bx->GetEntries() << " threshold " << thr / 3564. * bx->GetEntries()
0241                   << std::endl;
0242         for (int i = 1; i < bx->GetNbinsX() + 1; ++i) {
0243           if (bx->GetBinEntries(i) > thr / 3564. * bx->GetEntries()) {
0244             if (!filledbx.empty() && i == filledbx[filledbx.size() - 1] + 1) {
0245               std::cout << "This is not a 50ns run ! " << std::endl;
0246               filledbx.clear();
0247               return filledbx;
0248             }
0249             filledbx.push_back(i);
0250           }
0251         }
0252       } else {
0253         std::cout << "Histogram not found" << std::endl;
0254       }
0255     } else {
0256       std::cout << "module path is not ok" << std::endl;
0257     }
0258   } else {
0259     std::cout << "Input file pointer is not ok" << std::endl;
0260   }
0261
0262   //  std::cout << filledbx.size() << " filled bunch crossings" << std::endl;
0263   //  for(std::vector<int>::const_iterator fbx=filledbx.begin();fbx!=filledbx.end();++fbx) { std::cout << *fbx << std::endl;}
0264   return filledbx;
0265 }