Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-15 23:06:23 

0001 #include "TFile.h"
0002 #include "TMath.h"
0003 #include "TF1.h"
0004 #include "TH1F.h"
0005 #include "TH2F.h"
0006 #include "TKey.h"
0007 #include "TDirectory.h"
0008 #include "TCanvas.h"
0009 #include "TObject.h"
0010 #include "TStyle.h"
0011 #include "TSystem.h"
0012 #include "TClass.h"
0013 #include "TLegend.h"
0014 #include "TObjString.h"
0015 #include <string>
0016 #include <iomanip>
0017 #include "TPaveText.h"
0018 #include <fstream>  // std::ofstream
0019 #include <iostream>
0020 #include <algorithm>
0021 #include <vector>
0022 #include <map>
0023 #include <regex>
0024 #include <unordered_map>
0025 #define PLOTTING_MACRO  // to remove message logger
0026 #include "Alignment/OfflineValidation/interface/PVValidationHelpers.h"
0027 #include "Alignment/OfflineValidation/macros/CMS_lumi.h"
0028 
0029 /* 
0030    Store here the globals binning settings
0031 */
0032 
0033 float SUMPTMIN = 1.;
0034 float SUMPTMAX = 1e3;
0035 int TRACKBINS = 120;
0036 int VTXBINS = 60;
0037 
0038 /* 
0039    This is an auxilliary class to store the list of files
0040    to be used to plot
0041 */
0042 
0043 bool debugMode = false;
0044 
0045 class PVResolutionVariables {
0046 public:
0047   PVResolutionVariables(TString fileName, TString baseDir, TString legName = "", int color = 1, int style = 1);
0048   int getLineColor() { return lineColor; }
0049   int getLineStyle() { return lineStyle; }
0050   TString getName() { return legendName; }
0051   TFile* getFile() { return file; }
0052   TString getFileName() { return fname; }
0053 
0054 private:
0055   TFile* file;
0056   int lineColor;
0057   int lineStyle;
0058   TString legendName;
0059   TString fname;
0060 };
0061 
0062 PVResolutionVariables::PVResolutionVariables(
0063     TString fileName, TString baseDir, TString legName, int lColor, int lStyle) {
0064   fname = fileName;
0065   lineColor = lColor;
0066   lineStyle = lStyle % 100;
0067   if (legName == "") {
0068     std::string s_fileName = fileName.Data();
0069     int start = 0;
0070     if (s_fileName.find('/'))
0071       start = s_fileName.find_last_of('/') + 1;
0072     int stop = s_fileName.find_last_of('.');
0073     legendName = s_fileName.substr(start, stop - start);
0074   } else {
0075     legendName = legName;
0076   }
0077 
0078   // check if the base dir exists
0079   file = TFile::Open(fileName.Data(), "READ");
0080 
0081   if (!file) {
0082     std::cout << "ERROR! file " << fileName.Data() << " does not exist!" << std::endl;
0083     assert(false);
0084   }
0085 
0086   if (file->Get(baseDir.Data())) {
0087     std::cout << fileName.Data() << ", found base directory: " << baseDir.Data() << std::endl;
0088   } else {
0089     std::cout << fileName.Data() << ", no directory named: " << baseDir.Data() << std::endl;
0090     assert(false);
0091   }
0092 }
0093 
0094 namespace PVResolution {
0095   std::vector<PVResolutionVariables*> sourceList;
0096 
0097   // fill the list of files
0098   //*************************************************************
0099   void loadFileList(const char* inputFile, TString baseDir, TString legendName, int lineColor, int lineStyle)
0100   //*************************************************************
0101   {
0102     gErrorIgnoreLevel = kFatal;
0103     sourceList.push_back(new PVResolutionVariables(inputFile, baseDir, legendName, lineColor, lineStyle));
0104   }
0105 
0106   //*************************************************************
0107   void clearFileList()
0108   //*************************************************************
0109   {
0110     sourceList.clear();
0111   }
0112 }  // namespace PVResolution
0113 
0114 namespace statmode {
0115   using fitParams = std::pair<std::pair<double, double>, std::pair<double, double> >;
0116 }
0117 
0118 Int_t def_markers[9] = {kFullSquare,
0119                         kFullCircle,
0120                         kFullTriangleDown,
0121                         kOpenSquare,
0122                         kDot,
0123                         kOpenCircle,
0124                         kFullTriangleDown,
0125                         kFullTriangleUp,
0126                         kOpenTriangleDown};
0127 Int_t def_colors[9] = {kBlack, kRed, kBlue, kMagenta, kGreen, kCyan, kViolet, kOrange, kGreen + 2};
0128 
0129 // auxilliary functions
0130 void setPVResolStyle();
0131 void fillTrendPlotByIndex(TH1F* trendPlot,
0132                           std::unordered_map<std::string, TH1F*>& h,
0133                           std::regex toMatch,
0134                           PVValHelper::estimator fitPar_);
0135 statmode::fitParams fitResolutions(TH1* hist, bool singleTime = false);
0136 void makeNiceTrendPlotStyle(TH1* hist, Int_t color, Int_t style);
0137 void adjustMaximum(TH1F* histos[], int size);
0138 
0139 // inline function
0140 namespace pvresol {
0141   int check(const double a[], int n) {
0142     //while (--n > 0 && a[n] == a[0])    // exact match
0143     while (--n > 0 && (a[n] - a[0]) < 0.01)  // merged input files, protection agains numerical precision
0144       ;
0145     return n != 0;
0146   }
0147 }  // namespace pvresol
0148 
0149 // MAIN
0150 //*************************************************************
0151 void FitPVResolution(TString namesandlabels, TString theDate = "", bool isStrict = false) {
0152   //*************************************************************
0153 
0154   bool fromLoader = false;
0155   setPVResolStyle();
0156 
0157   // check if the loader is empty
0158   if (!PVResolution::sourceList.empty()) {
0159     fromLoader = true;
0160   }
0161 
0162   // if enters here, whatever is passed from command line is neglected
0163   if (fromLoader) {
0164     std::cout << "FitPVResiduals::FitPVResiduals(): file list specified from loader" << std::endl;
0165     std::cout << "======================================================" << std::endl;
0166     std::cout << "!!    arguments passed from CLI will be neglected   !!" << std::endl;
0167     std::cout << "======================================================" << std::endl;
0168     for (std::vector<PVResolutionVariables*>::iterator it = PVResolution::sourceList.begin();
0169          it != PVResolution::sourceList.end();
0170          ++it) {
0171       std::cout << "name:  " << std::setw(20) << (*it)->getName() << " |file:  " << std::setw(15) << (*it)->getFile()
0172                 << " |color: " << std::setw(5) << (*it)->getLineColor() << " |style: " << std::setw(5)
0173                 << (*it)->getLineStyle() << std::endl;
0174     }
0175     std::cout << "======================================================" << std::endl;
0176   }
0177 
0178   Int_t theFileCount = 0;
0179   TList* FileList = new TList();
0180   TList* LabelList = new TList();
0181 
0182   if (!fromLoader) {
0183     namesandlabels.Remove(TString::kTrailing, ',');
0184     TObjArray* nameandlabelpairs = namesandlabels.Tokenize(",");
0185     for (Int_t i = 0; i < nameandlabelpairs->GetEntries(); ++i) {
0186       TObjArray* aFileLegPair = TString(nameandlabelpairs->At(i)->GetName()).Tokenize("=");
0187 
0188       if (aFileLegPair->GetEntries() == 2) {
0189         FileList->Add(TFile::Open(aFileLegPair->At(0)->GetName(), "READ"));  // 2
0190         LabelList->Add(aFileLegPair->At(1));
0191       } else {
0192         std::cout << "Please give file name and legend entry in the following form:\n"
0193                   << " filename1=legendentry1,filename2=legendentry2\n";
0194         exit(EXIT_FAILURE);
0195       }
0196     }
0197 
0198     theFileCount = FileList->GetSize();
0199   } else {
0200     for (std::vector<PVResolutionVariables*>::iterator it = PVResolution::sourceList.begin();
0201          it != PVResolution::sourceList.end();
0202          ++it) {
0203       //FileList->Add((*it)->getFile()); // was extremely slow
0204       FileList->Add(TFile::Open((*it)->getFileName(), "READ"));
0205     }
0206     theFileCount = PVResolution::sourceList.size();
0207   }
0208 
0209   const Int_t nFiles_ = theFileCount;
0210   TString LegLabels[10];
0211   TFile* fins[nFiles_];
0212   Int_t markers[9];
0213   Int_t colors[9];
0214 
0215   for (Int_t j = 0; j < nFiles_; j++) {
0216     // Retrieve files
0217     fins[j] = (TFile*)FileList->At(j);
0218 
0219     if (!fromLoader) {
0220       TObjString* legend = (TObjString*)LabelList->At(j);
0221       LegLabels[j] = legend->String();
0222       markers[j] = def_markers[j];
0223       colors[j] = def_colors[j];
0224     } else {
0225       LegLabels[j] = PVResolution::sourceList[j]->getName();
0226       markers[j] = PVResolution::sourceList[j]->getLineStyle();
0227       colors[j] = PVResolution::sourceList[j]->getLineColor();
0228     }
0229 
0230     LegLabels[j].ReplaceAll("_", " ");
0231     std::cout << "FitPVResolution::FitPVResolution(): label[" << j << "] " << LegLabels[j] << std::endl;
0232   }
0233 
0234   // get the binnings
0235   TH1F* theBinHistos[nFiles_];
0236   double theSumPtMin_[nFiles_];
0237   double theSumPtMax_[nFiles_];
0238   double theTrackBINS_[nFiles_];
0239   double theVtxBINS_[nFiles_];
0240 
0241   for (Int_t i = 0; i < nFiles_; i++) {
0242     fins[i]->cd("PrimaryVertexResolution/BinningFeatures/");
0243     if (gDirectory->GetListOfKeys()->Contains("h_profileBinnings")) {
0244       gDirectory->GetObject("h_profileBinnings", theBinHistos[i]);
0245       theSumPtMin_[i] =
0246           theBinHistos[i]->GetBinContent(1) / (theBinHistos[i]->GetEntries() / theBinHistos[i]->GetNbinsX());
0247       std::cout << "File n. " << i << " has theSumPtMin[" << i << "] = " << theSumPtMin_[i] << std::endl;
0248       theSumPtMax_[i] =
0249           theBinHistos[i]->GetBinContent(2) / (theBinHistos[i]->GetEntries() / theBinHistos[i]->GetNbinsX());
0250       std::cout << "File n. " << i << " has theSumPtMax[" << i << "] = " << theSumPtMax_[i] << std::endl;
0251       theTrackBINS_[i] =
0252           theBinHistos[i]->GetBinContent(3) / (theBinHistos[i]->GetEntries() / theBinHistos[i]->GetNbinsX());
0253       std::cout << "File n. " << i << " has theTrackBINS[" << i << "] = " << theTrackBINS_[i] << std::endl;
0254       theVtxBINS_[i] =
0255           theBinHistos[i]->GetBinContent(4) / (theBinHistos[i]->GetEntries() / theBinHistos[i]->GetNbinsX());
0256       std::cout << "File n. " << i << " has theVtxBINS[" << i << "] = " << theVtxBINS_[i] << std::endl;
0257     } else {
0258       theSumPtMin_[i] = 1.;
0259       std::cout << "File n. " << i << " getting the default minimum sum pT range: " << theSumPtMin_[i] << std::endl;
0260       theSumPtMax_[i] = 1e3;
0261       std::cout << "File n. " << i << " getting the default maxmum sum pT range: " << theSumPtMax_[i] << std::endl;
0262       theTrackBINS_[i] = 120.;
0263       std::cout << "File n. " << i << " getting the default number of tracks bins: " << theTrackBINS_[i] << std::endl;
0264       theTrackBINS_[i] = 60.;
0265       std::cout << "File n. " << i << " getting the default number of vertices bins: " << theVtxBINS_[i] << std::endl;
0266     }
0267   }
0268 
0269   // checks if all minimum sum pT ranges coincide
0270   // if not, exits
0271   if (pvresol::check(theSumPtMin_, nFiles_)) {
0272     std::cout << "======================================================" << std::endl;
0273     std::cout << "FitPVResolution::FitPVResolution(): the minimum sum pT is different" << std::endl;
0274     std::cout << "exiting..." << std::endl;
0275     exit(EXIT_FAILURE);
0276   } else {
0277     SUMPTMIN = theSumPtMin_[0];
0278     std::cout << "======================================================" << std::endl;
0279     std::cout << "FitPVResolution::FitPVResolution(): the minimum sum pT is: " << SUMPTMIN << std::endl;
0280     std::cout << "======================================================" << std::endl;
0281   }
0282 
0283   // checks if all maximum sum pT ranges coincide
0284   // if not, exits
0285   if (pvresol::check(theSumPtMax_, nFiles_)) {
0286     std::cout << "======================================================" << std::endl;
0287     std::cout << "FitPVResolution::FitPVResolution(): the maximum sum pT is different" << std::endl;
0288     std::cout << "exiting..." << std::endl;
0289     exit(EXIT_FAILURE);
0290   } else {
0291     SUMPTMAX = theSumPtMax_[0];
0292     std::cout << "======================================================" << std::endl;
0293     std::cout << "FitPVResolution::FitPVResolution(): the maximum sum pT is: " << SUMPTMAX << std::endl;
0294     std::cout << "======================================================" << std::endl;
0295   }
0296 
0297   // checks if all number of tracks bins coincide
0298   // if not, exits
0299   if (pvresol::check(theTrackBINS_, nFiles_)) {
0300     std::cout << "======================================================" << std::endl;
0301     std::cout << "FitPVResolution::FitPVResolution(): the number of track bins is different" << std::endl;
0302     if (isStrict) {
0303       std::cout << "exiting..." << std::endl;
0304       std::cout << "======================================================" << std::endl;
0305       exit(EXIT_FAILURE);  //this is stricter
0306     } else {
0307       TRACKBINS = *std::max_element(theTrackBINS_, theTrackBINS_ + nFiles_);
0308       std::cout << "chosen the maximum: " << TRACKBINS << std::endl;
0309       std::cout << "======================================================" << std::endl;
0310     }
0311   } else {
0312     TRACKBINS = int(theTrackBINS_[0]);
0313     std::cout << "======================================================" << std::endl;
0314     std::cout << "FitPVResolution::FitPVResolution(): the number of track bins is: " << TRACKBINS << std::endl;
0315     std::cout << "======================================================" << std::endl;
0316   }
0317 
0318   // checks if all number of vertices bins coincide
0319   // if not, exits
0320   if (pvresol::check(theVtxBINS_, nFiles_)) {
0321     std::cout << "======================================================" << std::endl;
0322     std::cout << "FitPVResolution::FitPVResolution(): the number of vertices bins is different" << std::endl;
0323     if (isStrict) {
0324       std::cout << "exiting..." << std::endl;
0325       std::cout << "======================================================" << std::endl;
0326       exit(EXIT_FAILURE);  //this is stricter
0327     } else {
0328       VTXBINS = *std::max_element(theVtxBINS_, theVtxBINS_ + nFiles_);
0329       std::cout << "chosen the maximum: " << VTXBINS << std::endl;
0330       std::cout << "======================================================" << std::endl;
0331     }
0332   } else {
0333     VTXBINS = int(theVtxBINS_[0]);
0334     std::cout << "======================================================" << std::endl;
0335     std::cout << "FitPVResolution::FitPVResolution(): the number of vertices bins is: " << VTXBINS << std::endl;
0336     std::cout << "======================================================" << std::endl;
0337   }
0338 
0339   // max vertices
0340   const int max_n_vertices = std::min(60, VTXBINS);  // take the minimum to avoid overflow
0341   std::vector<float> myNVtx_bins_;
0342   for (float i = 0; i <= max_n_vertices; i++) {
0343     myNVtx_bins_.push_back(i - 0.5f);
0344   }
0345 
0346   // max track
0347   const int max_n_tracks = std::min(120, TRACKBINS);  // take the minimum to avoid overflow
0348   std::vector<float> myNTrack_bins_;
0349   for (float i = 0; i <= max_n_tracks; i++) {
0350     myNTrack_bins_.push_back(i - 0.5f);
0351   }
0352 
0353   // max sumPt
0354   const int max_sum_pt = 30;
0355   std::array<float, max_sum_pt + 1> mypT_bins_ = PVValHelper::makeLogBins<float, max_sum_pt>(SUMPTMIN, SUMPTMAX);
0356 
0357   // define the maps
0358   std::unordered_map<std::string, TH1F*> hpulls_;
0359   std::unordered_map<std::string, TH1F*> hdiffs_;
0360 
0361   // summary plots
0362 
0363   TH1F* p_resolX_vsSumPt_[nFiles_];
0364   TH1F* p_resolY_vsSumPt_[nFiles_];
0365   TH1F* p_resolZ_vsSumPt_[nFiles_];
0366 
0367   TH1F* p_resolX_vsNtracks_[nFiles_];
0368   TH1F* p_resolY_vsNtracks_[nFiles_];
0369   TH1F* p_resolZ_vsNtracks_[nFiles_];
0370 
0371   TH1F* p_resolX_vsNVtx_[nFiles_];
0372   TH1F* p_resolY_vsNVtx_[nFiles_];
0373   TH1F* p_resolZ_vsNVtx_[nFiles_];
0374 
0375   TH1F* p_pullX_vsSumPt_[nFiles_];
0376   TH1F* p_pullY_vsSumPt_[nFiles_];
0377   TH1F* p_pullZ_vsSumPt_[nFiles_];
0378 
0379   TH1F* p_pullX_vsNtracks_[nFiles_];
0380   TH1F* p_pullY_vsNtracks_[nFiles_];
0381   TH1F* p_pullZ_vsNtracks_[nFiles_];
0382 
0383   TH1F* p_pullX_vsNVtx_[nFiles_];
0384   TH1F* p_pullY_vsNVtx_[nFiles_];
0385   TH1F* p_pullZ_vsNVtx_[nFiles_];
0386 
0387   // load in the map all the relevant histograms
0388 
0389   for (Int_t j = 0; j < nFiles_; j++) {
0390     // vs n. tracks
0391 
0392     p_pullX_vsNtracks_[j] = new TH1F(Form("p_pullX_vsNtracks_%i", j),
0393                                      "x-pull vs n_{tracks};n_{tracks} in vertex; x vertex pull",
0394                                      myNTrack_bins_.size() - 1,
0395                                      myNTrack_bins_.data());
0396     p_pullY_vsNtracks_[j] = new TH1F(Form("p_pullY_vsNtracks_%i", j),
0397                                      "y-pull vs n_{tracks};n_{tracks} in vertex; y vertex pull",
0398                                      myNTrack_bins_.size() - 1,
0399                                      myNTrack_bins_.data());
0400     p_pullZ_vsNtracks_[j] = new TH1F(Form("p_pullZ_vsNtracks_%i", j),
0401                                      "z-pull vs n_{tracks};n_{tracks} in vertex; z vertex pull",
0402                                      myNTrack_bins_.size() - 1,
0403                                      myNTrack_bins_.data());
0404 
0405     p_resolX_vsNtracks_[j] = new TH1F(Form("p_resolX_vsNtracks_%i", j),
0406                                       "x-resolution vs n_{tracks};n_{tracks} in vertex; x vertex resolution [#mum]",
0407                                       myNTrack_bins_.size() - 1,
0408                                       myNTrack_bins_.data());
0409     p_resolY_vsNtracks_[j] = new TH1F(Form("p_resolY_vsNtracks_%i", j),
0410                                       "y-resolution vs n_{tracks};n_{tracks} in vertex; y vertex resolution [#mum]",
0411                                       myNTrack_bins_.size() - 1,
0412                                       myNTrack_bins_.data());
0413     p_resolZ_vsNtracks_[j] = new TH1F(Form("p_resolZ_vsNtracks_%i", j),
0414                                       "z-resolution vs n_{tracks};n_{tracks} in vertex; z vertex resolution [#mum]",
0415                                       myNTrack_bins_.size() - 1,
0416                                       myNTrack_bins_.data());
0417 
0418     for (int i = 0; i < max_n_tracks; i++) {
0419       hpulls_[Form("pullX_%dTrks_file_%i", i, j)] =
0420           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/xPullNtracks/histo_pullX_Ntracks_plot%i", i));
0421       hpulls_[Form("pullY_%dTrks_file_%i", i, j)] =
0422           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/yPullNtracks/histo_pullY_Ntracks_plot%i", i));
0423       hpulls_[Form("pullZ_%dTrks_file_%i", i, j)] =
0424           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/zPullNtracks/histo_pullZ_Ntracks_plot%i", i));
0425       hdiffs_[Form("diffX_%dTrks_file_%i", i, j)] =
0426           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/xResolNtracks/histo_resolX_Ntracks_plot%i", i));
0427       hdiffs_[Form("diffY_%dTrks_file_%i", i, j)] =
0428           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/yResolNtracks/histo_resolY_Ntracks_plot%i", i));
0429       hdiffs_[Form("diffZ_%dTrks_file_%i", i, j)] =
0430           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/zResolNtracks/histo_resolZ_Ntracks_plot%i", i));
0431     }
0432 
0433     // vs sumPt
0434 
0435     p_pullX_vsSumPt_[j] = new TH1F(Form("p_pullX_vsSumPt_%i", j),
0436                                    "x-pull vs #Sigma p_{T};#sum p_{T} [GeV]; x vertex pull",
0437                                    mypT_bins_.size() - 1,
0438                                    mypT_bins_.data());
0439     p_pullY_vsSumPt_[j] = new TH1F(Form("p_pullY_vsSumPt_%i", j),
0440                                    "y-pull vs #Sigma p_{T};#sum p_{T} [GeV]; y vertex pull",
0441                                    mypT_bins_.size() - 1,
0442                                    mypT_bins_.data());
0443     p_pullZ_vsSumPt_[j] = new TH1F(Form("p_pullZ_vsSumPt_%i", j),
0444                                    "z-pull vs #Sigma p_{T};#sum p_{T} [GeV]; z vertex pull",
0445                                    mypT_bins_.size() - 1,
0446                                    mypT_bins_.data());
0447 
0448     p_resolX_vsSumPt_[j] = new TH1F(Form("p_resolX_vsSumPt_%i", j),
0449                                     "x-resolution vs #Sigma p_{T};#sum p_{T} [GeV]; x vertex resolution [#mum]",
0450                                     mypT_bins_.size() - 1,
0451                                     mypT_bins_.data());
0452     p_resolY_vsSumPt_[j] = new TH1F(Form("p_resolY_vsSumPt_%i", j),
0453                                     "y-resolution vs #Sigma p_{T};#sum p_{T} [GeV]; y vertex resolution [#mum]",
0454                                     mypT_bins_.size() - 1,
0455                                     mypT_bins_.data());
0456     p_resolZ_vsSumPt_[j] = new TH1F(Form("p_resolZ_vsSumPt_%i", j),
0457                                     "z-resolution vs #Sigma p_{T};#sum p_{T} [GeV]; z vertex resolution [#mum]",
0458                                     mypT_bins_.size() - 1,
0459                                     mypT_bins_.data());
0460 
0461     for (int i = 0; i < max_sum_pt; i++) {
0462       hpulls_[Form("pullX_%dsumPt_file_%i", i, j)] =
0463           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/xPullSumPt/histo_pullX_sumPt_plot%i", i));
0464       hpulls_[Form("pullY_%dsumPt_file_%i", i, j)] =
0465           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/yPullSumPt/histo_pullY_sumPt_plot%i", i));
0466       hpulls_[Form("pullZ_%dsumPt_file_%i", i, j)] =
0467           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/zPullSumPt/histo_pullZ_sumPt_plot%i", i));
0468       hdiffs_[Form("diffX_%dsumPt_file_%i", i, j)] =
0469           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/xResolSumPt/histo_resolX_sumPt_plot%i", i));
0470       hdiffs_[Form("diffY_%dsumPt_file_%i", i, j)] =
0471           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/yResolSumPt/histo_resolY_sumPt_plot%i", i));
0472       hdiffs_[Form("diffZ_%dsumPt_file_%i", i, j)] =
0473           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/zResolSumPt/histo_resolZ_sumPt_plot%i", i));
0474     }
0475 
0476     // vs n. vertices
0477 
0478     p_pullX_vsNVtx_[j] = new TH1F(Form("p_pullX_vsNVtx_%i", j),
0479                                   "x-pull vs n_{vertices};n_{vertices} in event; x vertex pull",
0480                                   myNVtx_bins_.size() - 1,
0481                                   myNVtx_bins_.data());
0482     p_pullY_vsNVtx_[j] = new TH1F(Form("p_pullY_vsNVtx_%i", j),
0483                                   "y-pull vs n_{vertices};n_{vertices} in event; y vertex pull",
0484                                   myNVtx_bins_.size() - 1,
0485                                   myNVtx_bins_.data());
0486     p_pullZ_vsNVtx_[j] = new TH1F(Form("p_pullZ_vsNVtx_%i", j),
0487                                   "z-pull vs n_{vertices};n_{vertices} in event; z vertex pull",
0488                                   myNVtx_bins_.size() - 1,
0489                                   myNVtx_bins_.data());
0490 
0491     p_resolX_vsNVtx_[j] = new TH1F(Form("p_resolX_vsNVtx_%i", j),
0492                                    "x-resolution vs n_{vertices};n_{vertices} in event; x vertex resolution [#mum]",
0493                                    myNVtx_bins_.size() - 1,
0494                                    myNVtx_bins_.data());
0495     p_resolY_vsNVtx_[j] = new TH1F(Form("p_resolY_vsNVtx_%i", j),
0496                                    "y-resolution vs n_{vertices};n_{vertices} in event; y vertex resolution [#mum]",
0497                                    myNVtx_bins_.size() - 1,
0498                                    myNVtx_bins_.data());
0499     p_resolZ_vsNVtx_[j] = new TH1F(Form("p_resolZ_vsNVtx_%i", j),
0500                                    "z-resolution vs n_{vertices};n_{vertices} in event; z vertex resolution [#mum]",
0501                                    myNVtx_bins_.size() - 1,
0502                                    myNVtx_bins_.data());
0503 
0504     for (int i = 0; i < max_n_vertices; i++) {
0505       hpulls_[Form("pullX_%dNVtx_file_%i", i, j)] =
0506           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/xPullNvtx/histo_pullX_Nvtx_plot%i", i));
0507       hpulls_[Form("pullY_%dNVtx_file_%i", i, j)] =
0508           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/yPullNvtx/histo_pullY_Nvtx_plot%i", i));
0509       hpulls_[Form("pullZ_%dNVtx_file_%i", i, j)] =
0510           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/zPullNvtx/histo_pullZ_Nvtx_plot%i", i));
0511       hdiffs_[Form("diffX_%dNVtx_file_%i", i, j)] =
0512           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/xResolNvtx/histo_resolX_Nvtx_plot%i", i));
0513       hdiffs_[Form("diffY_%dNVtx_file_%i", i, j)] =
0514           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/yResolNvtx/histo_resolY_Nvtx_plot%i", i));
0515       hdiffs_[Form("diffZ_%dNVtx_file_%i", i, j)] =
0516           (TH1F*)fins[j]->Get(Form("PrimaryVertexResolution/zResolNvtx/histo_resolZ_Nvtx_plot%i", i));
0517     }
0518   }
0519 
0520   // dump the list of keys and check all needed histograms are available
0521   for (const auto& key : hpulls_) {
0522     if (key.second == nullptr) {
0523       std::cout << "!!!WARNING!!! FitPVResolution::FitPVResolution(): missing histograms for key " << key.first
0524                 << ". I am going to exit. Good-bye!" << std::endl;
0525       exit(EXIT_FAILURE);
0526     }
0527   }
0528 
0529   // compute and store the trend plots
0530   for (Int_t j = 0; j < nFiles_; j++) {
0531     // resolutions
0532     fillTrendPlotByIndex(p_resolX_vsSumPt_[j],
0533                          hdiffs_,
0534                          std::regex(("diffX_(.*)sumPt_file_" + std::to_string(j)).c_str()),
0535                          PVValHelper::WIDTH);
0536     fillTrendPlotByIndex(p_resolY_vsSumPt_[j],
0537                          hdiffs_,
0538                          std::regex(("diffY_(.*)sumPt_file_" + std::to_string(j)).c_str()),
0539                          PVValHelper::WIDTH);
0540     fillTrendPlotByIndex(p_resolZ_vsSumPt_[j],
0541                          hdiffs_,
0542                          std::regex(("diffZ_(.*)sumPt_file_" + std::to_string(j)).c_str()),
0543                          PVValHelper::WIDTH);
0544 
0545     fillTrendPlotByIndex(p_resolX_vsNtracks_[j],
0546                          hdiffs_,
0547                          std::regex(("diffX_(.*)Trks_file_" + std::to_string(j)).c_str()),
0548                          PVValHelper::WIDTH);
0549     fillTrendPlotByIndex(p_resolY_vsNtracks_[j],
0550                          hdiffs_,
0551                          std::regex(("diffY_(.*)Trks_file_" + std::to_string(j)).c_str()),
0552                          PVValHelper::WIDTH);
0553     fillTrendPlotByIndex(p_resolZ_vsNtracks_[j],
0554                          hdiffs_,
0555                          std::regex(("diffZ_(.*)Trks_file_" + std::to_string(j)).c_str()),
0556                          PVValHelper::WIDTH);
0557 
0558     fillTrendPlotByIndex(p_resolX_vsNVtx_[j],
0559                          hdiffs_,
0560                          std::regex(("diffX_(.*)NVtx_file_" + std::to_string(j)).c_str()),
0561                          PVValHelper::WIDTH);
0562     fillTrendPlotByIndex(p_resolY_vsNVtx_[j],
0563                          hdiffs_,
0564                          std::regex(("diffY_(.*)NVtx_file_" + std::to_string(j)).c_str()),
0565                          PVValHelper::WIDTH);
0566     fillTrendPlotByIndex(p_resolZ_vsNVtx_[j],
0567                          hdiffs_,
0568                          std::regex(("diffZ_(.*)NVtx_file_" + std::to_string(j)).c_str()),
0569                          PVValHelper::WIDTH);
0570 
0571     // pulls
0572 
0573     fillTrendPlotByIndex(p_pullX_vsSumPt_[j],
0574                          hpulls_,
0575                          std::regex(("pullX_(.*)sumPt_file_" + std::to_string(j)).c_str()),
0576                          PVValHelper::WIDTH);
0577     fillTrendPlotByIndex(p_pullY_vsSumPt_[j],
0578                          hpulls_,
0579                          std::regex(("pullY_(.*)sumPt_file_" + std::to_string(j)).c_str()),
0580                          PVValHelper::WIDTH);
0581     fillTrendPlotByIndex(p_pullZ_vsSumPt_[j],
0582                          hpulls_,
0583                          std::regex(("pullZ_(.*)sumPt_file_" + std::to_string(j)).c_str()),
0584                          PVValHelper::WIDTH);
0585 
0586     fillTrendPlotByIndex(p_pullX_vsNtracks_[j],
0587                          hpulls_,
0588                          std::regex(("pullX_(.*)Trks_file_" + std::to_string(j)).c_str()),
0589                          PVValHelper::WIDTH);
0590     fillTrendPlotByIndex(p_pullY_vsNtracks_[j],
0591                          hpulls_,
0592                          std::regex(("pullY_(.*)Trks_file_" + std::to_string(j)).c_str()),
0593                          PVValHelper::WIDTH);
0594     fillTrendPlotByIndex(p_pullZ_vsNtracks_[j],
0595                          hpulls_,
0596                          std::regex(("pullZ_(.*)Trks_file_" + std::to_string(j)).c_str()),
0597                          PVValHelper::WIDTH);
0598 
0599     fillTrendPlotByIndex(p_pullX_vsNVtx_[j],
0600                          hpulls_,
0601                          std::regex(("pullX_(.*)NVtx_file_" + std::to_string(j)).c_str()),
0602                          PVValHelper::WIDTH);
0603     fillTrendPlotByIndex(p_pullY_vsNVtx_[j],
0604                          hpulls_,
0605                          std::regex(("pullY_(.*)NVtx_file_" + std::to_string(j)).c_str()),
0606                          PVValHelper::WIDTH);
0607     fillTrendPlotByIndex(p_pullZ_vsNVtx_[j],
0608                          hpulls_,
0609                          std::regex(("pullZ_(.*)NVtx_file_" + std::to_string(j)).c_str()),
0610                          PVValHelper::WIDTH);
0611 
0612     // beautify
0613 
0614     makeNiceTrendPlotStyle(p_resolX_vsSumPt_[j], colors[j], markers[j]);
0615     makeNiceTrendPlotStyle(p_resolY_vsSumPt_[j], colors[j], markers[j]);
0616     makeNiceTrendPlotStyle(p_resolZ_vsSumPt_[j], colors[j], markers[j]);
0617 
0618     makeNiceTrendPlotStyle(p_resolX_vsNtracks_[j], colors[j], markers[j]);
0619     makeNiceTrendPlotStyle(p_resolY_vsNtracks_[j], colors[j], markers[j]);
0620     makeNiceTrendPlotStyle(p_resolZ_vsNtracks_[j], colors[j], markers[j]);
0621 
0622     makeNiceTrendPlotStyle(p_resolX_vsNVtx_[j], colors[j], markers[j]);
0623     makeNiceTrendPlotStyle(p_resolY_vsNVtx_[j], colors[j], markers[j]);
0624     makeNiceTrendPlotStyle(p_resolZ_vsNVtx_[j], colors[j], markers[j]);
0625 
0626     // pulls
0627 
0628     makeNiceTrendPlotStyle(p_pullX_vsSumPt_[j], colors[j], markers[j]);
0629     makeNiceTrendPlotStyle(p_pullY_vsSumPt_[j], colors[j], markers[j]);
0630     makeNiceTrendPlotStyle(p_pullZ_vsSumPt_[j], colors[j], markers[j]);
0631 
0632     makeNiceTrendPlotStyle(p_pullX_vsNtracks_[j], colors[j], markers[j]);
0633     makeNiceTrendPlotStyle(p_pullY_vsNtracks_[j], colors[j], markers[j]);
0634     makeNiceTrendPlotStyle(p_pullZ_vsNtracks_[j], colors[j], markers[j]);
0635 
0636     makeNiceTrendPlotStyle(p_pullX_vsNVtx_[j], colors[j], markers[j]);
0637     makeNiceTrendPlotStyle(p_pullY_vsNVtx_[j], colors[j], markers[j]);
0638     makeNiceTrendPlotStyle(p_pullZ_vsNVtx_[j], colors[j], markers[j]);
0639   }
0640 
0641   // adjust the maxima
0642 
0643   adjustMaximum(p_resolX_vsSumPt_, nFiles_);
0644   adjustMaximum(p_resolY_vsSumPt_, nFiles_);
0645   adjustMaximum(p_resolZ_vsSumPt_, nFiles_);
0646 
0647   adjustMaximum(p_resolX_vsNtracks_, nFiles_);
0648   adjustMaximum(p_resolY_vsNtracks_, nFiles_);
0649   adjustMaximum(p_resolZ_vsNtracks_, nFiles_);
0650 
0651   adjustMaximum(p_resolX_vsNVtx_, nFiles_);
0652   adjustMaximum(p_resolY_vsNVtx_, nFiles_);
0653   adjustMaximum(p_resolZ_vsNVtx_, nFiles_);
0654 
0655   adjustMaximum(p_pullX_vsSumPt_, nFiles_);
0656   adjustMaximum(p_pullY_vsSumPt_, nFiles_);
0657   adjustMaximum(p_pullZ_vsSumPt_, nFiles_);
0658 
0659   adjustMaximum(p_pullX_vsNtracks_, nFiles_);
0660   adjustMaximum(p_pullY_vsNtracks_, nFiles_);
0661   adjustMaximum(p_pullZ_vsNtracks_, nFiles_);
0662 
0663   adjustMaximum(p_pullX_vsNVtx_, nFiles_);
0664   adjustMaximum(p_pullY_vsNVtx_, nFiles_);
0665   adjustMaximum(p_pullZ_vsNVtx_, nFiles_);
0666 
0667   TCanvas* c1 = new TCanvas("VertexResolVsSumPt", "Vertex Resolution vs #sum p_{T} [GeV]", 1500, 700);
0668   c1->Divide(3, 1);
0669   TCanvas* c2 = new TCanvas("VertexPullVsSumPt", "Vertex Resolution vs #sum p_{T} [GeV]", 1500, 700);
0670   c2->Divide(3, 1);
0671   TCanvas* c3 = new TCanvas("VertexResolVsNTracks", "Vertex Resolution vs n. tracks", 1500, 700);
0672   c3->Divide(3, 1);
0673   TCanvas* c4 = new TCanvas("VertexPullVsNTracks", "Vertex Resolution vs n. tracks", 1500, 700);
0674   c4->Divide(3, 1);
0675   TCanvas* c5 = new TCanvas("VertexResolVsNVtx", "Vertex Resolution vs n. vertices", 1500, 700);
0676   c5->Divide(3, 1);
0677   TCanvas* c6 = new TCanvas("VertexPullVsNVtx", "Vertex Resolution vs n. vertices", 1500, 700);
0678   c6->Divide(3, 1);
0679 
0680   for (Int_t c = 1; c <= 3; c++) {
0681     c1->cd(c)->SetBottomMargin(0.14);
0682     c1->cd(c)->SetLeftMargin(0.15);
0683     c1->cd(c)->SetRightMargin(0.05);
0684     c1->cd(c)->SetTopMargin(0.05);
0685 
0686     c2->cd(c)->SetBottomMargin(0.14);
0687     c2->cd(c)->SetLeftMargin(0.15);
0688     c2->cd(c)->SetRightMargin(0.05);
0689     c2->cd(c)->SetTopMargin(0.05);
0690 
0691     c3->cd(c)->SetBottomMargin(0.14);
0692     c3->cd(c)->SetLeftMargin(0.15);
0693     c3->cd(c)->SetRightMargin(0.05);
0694     c3->cd(c)->SetTopMargin(0.05);
0695 
0696     c4->cd(c)->SetBottomMargin(0.14);
0697     c4->cd(c)->SetLeftMargin(0.15);
0698     c4->cd(c)->SetRightMargin(0.05);
0699     c4->cd(c)->SetTopMargin(0.05);
0700 
0701     c5->cd(c)->SetBottomMargin(0.14);
0702     c5->cd(c)->SetLeftMargin(0.15);
0703     c5->cd(c)->SetRightMargin(0.05);
0704     c5->cd(c)->SetTopMargin(0.05);
0705 
0706     c6->cd(c)->SetBottomMargin(0.14);
0707     c6->cd(c)->SetLeftMargin(0.15);
0708     c6->cd(c)->SetRightMargin(0.05);
0709     c6->cd(c)->SetTopMargin(0.05);
0710   }
0711 
0712   TLegend* lego = new TLegend(0.18, 0.80, 0.79, 0.93);
0713   // might be useful if many objects are compared
0714   if (nFiles_ > 4) {
0715     lego->SetNColumns(2);
0716   }
0717 
0718   lego->SetFillColor(10);
0719   if (nFiles_ > 3) {
0720     lego->SetTextSize(0.032);
0721   } else {
0722     lego->SetTextSize(0.042);
0723   }
0724   lego->SetTextFont(42);
0725   lego->SetFillColor(10);
0726   lego->SetLineColor(10);
0727   lego->SetShadowColor(10);
0728 
0729   TPaveText* ptDate = new TPaveText(0.17, 0.96, 0.50, 0.99, "blNDC");
0730   ptDate->SetFillColor(kYellow);
0731   //ptDate->SetFillColor(10);
0732   ptDate->SetBorderSize(1);
0733   ptDate->SetLineColor(kBlue);
0734   ptDate->SetLineWidth(1);
0735   ptDate->SetTextFont(32);
0736   TText* textDate = ptDate->AddText(theDate);
0737   textDate->SetTextSize(0.04);
0738   textDate->SetTextColor(kBlue);
0739 
0740   for (Int_t j = 0; j < nFiles_; j++) {
0741     // first canvas
0742 
0743     //p_resolX_vsSumPt_[j]->GetXaxis()->SetRangeUser(10., 400.);
0744     //p_resolY_vsSumPt_[j]->GetXaxis()->SetRangeUser(10., 400.);
0745     //p_resolZ_vsSumPt_[j]->GetXaxis()->SetRangeUser(10., 400.);
0746 
0747     c1->cd(1);
0748     j == 0 ? p_resolX_vsSumPt_[j]->Draw("E1") : p_resolX_vsSumPt_[j]->Draw("E1same");
0749     lego->AddEntry(p_resolX_vsSumPt_[j], LegLabels[j]);
0750 
0751     if (j == nFiles_ - 1)
0752       lego->Draw("same");
0753     if (theDate.Length() != 0)
0754       ptDate->Draw("same");
0755     TPad* current_pad = static_cast<TPad*>(c1->GetPad(1));
0756     CMS_lumi(current_pad, 6, 33);
0757 
0758     c1->cd(2);
0759     j == 0 ? p_resolY_vsSumPt_[j]->Draw("E1") : p_resolY_vsSumPt_[j]->Draw("E1same");
0760 
0761     if (j == nFiles_ - 1)
0762       lego->Draw("same");
0763     if (theDate.Length() != 0)
0764       ptDate->Draw("same");
0765     current_pad = static_cast<TPad*>(c1->GetPad(2));
0766     CMS_lumi(current_pad, 6, 33);
0767 
0768     c1->cd(3);
0769     j == 0 ? p_resolZ_vsSumPt_[j]->Draw("E1") : p_resolZ_vsSumPt_[j]->Draw("E1same");
0770 
0771     if (j == nFiles_ - 1)
0772       lego->Draw("same");
0773     if (theDate.Length() != 0)
0774       ptDate->Draw("same");
0775     current_pad = static_cast<TPad*>(c1->GetPad(3));
0776     CMS_lumi(current_pad, 6, 33);
0777 
0778     // second canvas
0779 
0780     c2->cd(1);
0781     j == 0 ? p_pullX_vsSumPt_[j]->Draw("E1") : p_pullX_vsSumPt_[j]->Draw("E1same");
0782 
0783     if (j == nFiles_ - 1)
0784       lego->Draw("same");
0785     if (theDate.Length() != 0)
0786       ptDate->Draw("same");
0787     current_pad = static_cast<TPad*>(c2->GetPad(1));
0788     CMS_lumi(current_pad, 6, 33);
0789 
0790     c2->cd(2);
0791     j == 0 ? p_pullY_vsSumPt_[j]->Draw("E1") : p_pullY_vsSumPt_[j]->Draw("E1same");
0792 
0793     if (j == nFiles_ - 1)
0794       lego->Draw("same");
0795     if (theDate.Length() != 0)
0796       ptDate->Draw("same");
0797     current_pad = static_cast<TPad*>(c2->GetPad(2));
0798     CMS_lumi(current_pad, 6, 33);
0799 
0800     c2->cd(3);
0801     j == 0 ? p_pullZ_vsSumPt_[j]->Draw("E1") : p_pullZ_vsSumPt_[j]->Draw("E1same");
0802 
0803     if (j == nFiles_ - 1)
0804       lego->Draw("same");
0805     if (theDate.Length() != 0)
0806       ptDate->Draw("same");
0807     current_pad = static_cast<TPad*>(c2->GetPad(3));
0808     CMS_lumi(current_pad, 6, 33);
0809 
0810     // third canvas
0811 
0812     c3->cd(1);
0813     j == 0 ? p_resolX_vsNtracks_[j]->Draw("E1") : p_resolX_vsNtracks_[j]->Draw("E1same");
0814 
0815     if (j == nFiles_ - 1)
0816       lego->Draw("same");
0817     if (theDate.Length() != 0)
0818       ptDate->Draw("same");
0819     current_pad = static_cast<TPad*>(c3->GetPad(1));
0820     CMS_lumi(current_pad, 6, 33);
0821 
0822     c3->cd(2);
0823     j == 0 ? p_resolY_vsNtracks_[j]->Draw("E1") : p_resolY_vsNtracks_[j]->Draw("E1same");
0824 
0825     if (j == nFiles_ - 1)
0826       lego->Draw("same");
0827     if (theDate.Length() != 0)
0828       ptDate->Draw("same");
0829     current_pad = static_cast<TPad*>(c3->GetPad(2));
0830     CMS_lumi(current_pad, 6, 33);
0831 
0832     c3->cd(3);
0833     j == 0 ? p_resolZ_vsNtracks_[j]->Draw("E1") : p_resolZ_vsNtracks_[j]->Draw("E1same");
0834 
0835     if (j == nFiles_ - 1)
0836       lego->Draw("same");
0837     if (theDate.Length() != 0)
0838       ptDate->Draw("same");
0839     current_pad = static_cast<TPad*>(c3->GetPad(3));
0840     CMS_lumi(current_pad, 6, 33);
0841 
0842     // fourth canvas
0843 
0844     c4->cd(1);
0845     j == 0 ? p_pullX_vsNtracks_[j]->Draw("E1") : p_pullX_vsNtracks_[j]->Draw("E1same");
0846 
0847     if (j == nFiles_ - 1)
0848       lego->Draw("same");
0849     if (theDate.Length() != 0)
0850       ptDate->Draw("same");
0851     current_pad = static_cast<TPad*>(c4->GetPad(1));
0852     CMS_lumi(current_pad, 6, 33);
0853 
0854     c4->cd(2);
0855     j == 0 ? p_pullY_vsNtracks_[j]->Draw("E1") : p_pullY_vsNtracks_[j]->Draw("E1same");
0856 
0857     if (j == nFiles_ - 1)
0858       lego->Draw("same");
0859     if (theDate.Length() != 0)
0860       ptDate->Draw("same");
0861     current_pad = static_cast<TPad*>(c4->GetPad(2));
0862     CMS_lumi(current_pad, 6, 33);
0863 
0864     c4->cd(3);
0865     j == 0 ? p_pullZ_vsNtracks_[j]->Draw("E1") : p_pullZ_vsNtracks_[j]->Draw("E1same");
0866 
0867     if (j == nFiles_ - 1)
0868       lego->Draw("same");
0869     if (theDate.Length() != 0)
0870       ptDate->Draw("same");
0871     current_pad = static_cast<TPad*>(c4->GetPad(3));
0872     CMS_lumi(current_pad, 6, 33);
0873 
0874     // fifth canvas
0875 
0876     c5->cd(1);
0877     j == 0 ? p_resolX_vsNVtx_[j]->Draw("E1") : p_resolX_vsNVtx_[j]->Draw("E1same");
0878 
0879     if (j == nFiles_ - 1)
0880       lego->Draw("same");
0881     if (theDate.Length() != 0)
0882       ptDate->Draw("same");
0883     current_pad = static_cast<TPad*>(c5->GetPad(1));
0884     CMS_lumi(current_pad, 6, 33);
0885 
0886     c5->cd(2);
0887     j == 0 ? p_resolY_vsNVtx_[j]->Draw("E1") : p_resolY_vsNVtx_[j]->Draw("E1same");
0888 
0889     if (j == nFiles_ - 1)
0890       lego->Draw("same");
0891     if (theDate.Length() != 0)
0892       ptDate->Draw("same");
0893     current_pad = static_cast<TPad*>(c5->GetPad(2));
0894     CMS_lumi(current_pad, 6, 33);
0895 
0896     c5->cd(3);
0897     j == 0 ? p_resolZ_vsNVtx_[j]->Draw("E1") : p_resolZ_vsNVtx_[j]->Draw("E1same");
0898 
0899     if (j == nFiles_ - 1)
0900       lego->Draw("same");
0901     if (theDate.Length() != 0)
0902       ptDate->Draw("same");
0903     current_pad = static_cast<TPad*>(c5->GetPad(3));
0904     CMS_lumi(current_pad, 6, 33);
0905 
0906     // sixth canvas
0907 
0908     c6->cd(1);
0909     j == 0 ? p_pullX_vsNVtx_[j]->Draw("E1") : p_pullX_vsNVtx_[j]->Draw("E1same");
0910 
0911     if (j == nFiles_ - 1)
0912       lego->Draw("same");
0913     if (theDate.Length() != 0)
0914       ptDate->Draw("same");
0915     current_pad = static_cast<TPad*>(c6->GetPad(1));
0916     CMS_lumi(current_pad, 6, 33);
0917 
0918     c6->cd(2);
0919     j == 0 ? p_pullY_vsNVtx_[j]->Draw("E1") : p_pullY_vsNVtx_[j]->Draw("E1same");
0920 
0921     if (j == nFiles_ - 1)
0922       lego->Draw("same");
0923     if (theDate.Length() != 0)
0924       ptDate->Draw("same");
0925     current_pad = static_cast<TPad*>(c6->GetPad(2));
0926     CMS_lumi(current_pad, 6, 33);
0927 
0928     c6->cd(3);
0929     j == 0 ? p_pullZ_vsNVtx_[j]->Draw("E1") : p_pullZ_vsNVtx_[j]->Draw("E1same");
0930 
0931     if (j == nFiles_ - 1)
0932       lego->Draw("same");
0933     if (theDate.Length() != 0)
0934       ptDate->Draw("same");
0935     current_pad = static_cast<TPad*>(c6->GetPad(3));
0936     CMS_lumi(current_pad, 6, 33);
0937   }
0938 
0939   if (theDate.Length() != 0)
0940     theDate.Prepend("_");
0941 
0942   TString theStrAlignment = LegLabels[0];
0943   for (Int_t j = 1; j < nFiles_; j++) {
0944     theStrAlignment += ("_vs_" + LegLabels[j]);
0945   }
0946   theStrAlignment.ReplaceAll(" ", "_");
0947 
0948   c1->SaveAs("VertexResolVsSumPt_" + theStrAlignment + theDate + ".png");
0949   c2->SaveAs("VertexPullVsSumPt_" + theStrAlignment + theDate + ".png");
0950   c3->SaveAs("VertexResolVsNTracks_" + theStrAlignment + theDate + ".png");
0951   c4->SaveAs("VertexPullVsNTracks_" + theStrAlignment + theDate + ".png");
0952   c5->SaveAs("VertexResolVsNVertices_" + theStrAlignment + theDate + ".png");
0953   c6->SaveAs("VertexPullVsNVertices_" + theStrAlignment + theDate + ".png");
0954 
0955   c1->SaveAs("VertexResolVsSumPt_" + theStrAlignment + theDate + ".pdf");
0956   c2->SaveAs("VertexPullVsSumPt_" + theStrAlignment + theDate + ".pdf");
0957   c3->SaveAs("VertexResolVsNTracks_" + theStrAlignment + theDate + ".pdf");
0958   c4->SaveAs("VertexPullVsNTracks_" + theStrAlignment + theDate + ".pdf");
0959   c5->SaveAs("VertexResolVsNVertices_" + theStrAlignment + theDate + ".pdf");
0960   c6->SaveAs("VertexPullVsNVertices_" + theStrAlignment + theDate + ".pdf");
0961 
0962   delete c1;
0963   delete c2;
0964   delete c3;
0965   delete c4;
0966   delete c5;
0967   delete c6;
0968 
0969   // delete everything in the source list
0970   for (std::vector<PVResolutionVariables*>::iterator it = PVResolution::sourceList.begin();
0971        it != PVResolution::sourceList.end();
0972        ++it) {
0973     delete (*it);
0974   }
0975 }
0976 
0977 /*--------------------------------------------------------------------*/
0978 void fillTrendPlotByIndex(TH1F* trendPlot,
0979                           std::unordered_map<std::string, TH1F*>& h,
0980                           std::regex toMatch,
0981                           PVValHelper::estimator fitPar_)
0982 /*--------------------------------------------------------------------*/
0983 {
0984   for (const auto& iterator : h) {
0985     statmode::fitParams myFit = fitResolutions(iterator.second, false);
0986 
0987     int bin = -1;
0988     std::string result;
0989     try {
0990       std::smatch match;
0991       if (std::regex_search(iterator.first, match, toMatch) && match.size() > 1) {
0992         result = match.str(1);
0993         bin = std::stoi(result) + 1;
0994       } else {
0995         result = std::string("");
0996         continue;
0997       }
0998     } catch (std::regex_error& e) {
0999       // Syntax error in the regular expression
1000     }
1001 
1002     switch (fitPar_) {
1003       case PVValHelper::MEAN: {
1004         float mean_ = myFit.first.first;
1005         float meanErr_ = myFit.first.second;
1006         trendPlot->SetBinContent(bin, mean_);
1007         trendPlot->SetBinError(bin, meanErr_);
1008         break;
1009       }
1010       case PVValHelper::WIDTH: {
1011         float width_ = myFit.second.first;
1012         float widthErr_ = myFit.second.second;
1013         trendPlot->SetBinContent(bin, width_);
1014         trendPlot->SetBinError(bin, widthErr_);
1015         break;
1016       }
1017       case PVValHelper::MEDIAN: {
1018         float median_ = PVValHelper::getMedian(iterator.second).value();
1019         float medianErr_ = PVValHelper::getMedian(iterator.second).error();
1020         trendPlot->SetBinContent(bin, median_);
1021         trendPlot->SetBinError(bin, medianErr_);
1022         break;
1023       }
1024       case PVValHelper::MAD: {
1025         float mad_ = PVValHelper::getMAD(iterator.second).value();
1026         float madErr_ = PVValHelper::getMAD(iterator.second).error();
1027         trendPlot->SetBinContent(bin, mad_);
1028         trendPlot->SetBinError(bin, madErr_);
1029         break;
1030       }
1031       default:
1032         std::cout << "fillTrendPlotByIndex() " << fitPar_ << " unknown estimator!" << std::endl;
1033         break;
1034     }
1035   }
1036 }
1037 
1038 /*--------------------------------------------------------------------*/
1039 statmode::fitParams fitResolutions(TH1* hist, bool singleTime)
1040 /*--------------------------------------------------------------------*/
1041 {
1042   if (hist->GetEntries() < 10) {
1043     // std::cout<<"hist name: "<<hist->GetName() << std::endl;
1044     return std::make_pair(std::make_pair(0., 0.), std::make_pair(0., 0.));
1045   }
1046 
1047   float maxHist = hist->GetXaxis()->GetXmax();
1048   float minHist = hist->GetXaxis()->GetXmin();
1049   float mean = hist->GetMean();
1050   float sigma = hist->GetRMS();
1051 
1052   if (TMath::IsNaN(mean) || TMath::IsNaN(sigma)) {
1053     mean = 0;
1054     //sigma= - hist->GetXaxis()->GetBinLowEdge(1) + hist->GetXaxis()->GetBinLowEdge(hist->GetNbinsX()+1);
1055     sigma = -minHist + maxHist;
1056     std::cout << "FitPVResolution::fitResolutions(): histogram" << hist->GetName() << " mean or sigma are NaN!!"
1057               << std::endl;
1058   }
1059 
1060   TF1 func("tmp", "gaus", mean - 2. * sigma, mean + 2. * sigma);
1061   if (0 == hist->Fit(&func, "QNR")) {  // N: do not blow up file by storing fit!
1062     mean = func.GetParameter(1);
1063     sigma = func.GetParameter(2);
1064 
1065     if (!singleTime) {
1066       // Check if histogram is weighted
1067       double sumWeights = hist->GetSumOfWeights();
1068       double effectiveEntries = hist->GetEffectiveEntries();
1069       bool isWeighted = !(sumWeights == effectiveEntries);
1070 
1071       if (isWeighted && debugMode) {
1072         std::cout << "A weighted input histogram has been provided, will use least squares fit instead of likelihood!"
1073                   << " Sum of weights: " << sumWeights << " effective entries: " << hist->GetEffectiveEntries()
1074                   << std::endl;
1075       }
1076       // If histogram is weighted, exclude the "L" option (Likelihood fit)
1077       std::string fitOptions = isWeighted ? "Q0R" : "Q0LR";
1078 
1079       // second fit: three sigma of first fit around mean of first fit
1080       func.SetRange(std::max(mean - 3 * sigma, minHist), std::min(mean + 3 * sigma, maxHist));
1081 
1082       // Perform fit with the appropriate options
1083       // I: integral gives more correct results if binning is too wide
1084       // L: Likelihood can treat empty bins correctly (if hist not weighted...)
1085       if (0 == hist->Fit(&func, fitOptions.c_str())) {
1086         if (hist->GetFunction(func.GetName())) {  // Take care that it is later on drawn:
1087           hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
1088         }
1089       }
1090     }
1091   }
1092 
1093   return std::make_pair(std::make_pair(func.GetParameter(1), func.GetParError(1)),
1094                         std::make_pair(func.GetParameter(2), func.GetParError(2)));
1095 }
1096 
1097 /*--------------------------------------------------------------------*/
1098 void makeNiceTrendPlotStyle(TH1* hist, Int_t color, Int_t style)
1099 /*--------------------------------------------------------------------*/
1100 {
1101   hist->SetStats(kFALSE);
1102   hist->SetLineWidth(2);
1103   hist->GetXaxis()->CenterTitle(true);
1104   hist->GetYaxis()->CenterTitle(true);
1105   hist->GetXaxis()->SetTitleFont(42);
1106   hist->GetYaxis()->SetTitleFont(42);
1107   hist->GetXaxis()->SetTitleSize(0.065);
1108   hist->GetYaxis()->SetTitleSize(0.065);
1109   hist->GetXaxis()->SetTitleOffset(1.0);
1110   hist->GetYaxis()->SetTitleOffset(1.2);
1111   hist->GetXaxis()->SetLabelFont(42);
1112   hist->GetYaxis()->SetLabelFont(42);
1113   hist->GetYaxis()->SetLabelSize(.055);
1114   hist->GetXaxis()->SetLabelSize(.055);
1115   hist->GetXaxis()->SetNdivisions(505);
1116   if (color != 8) {
1117     hist->SetMarkerSize(1.2);
1118   } else {
1119     hist->SetLineWidth(3);
1120     hist->SetMarkerSize(0.0);
1121   }
1122   hist->SetMarkerStyle(style);
1123   hist->SetLineColor(color);
1124   hist->SetMarkerColor(color);
1125 
1126   if (TString(hist->GetName()).Contains("pull"))
1127     hist->GetYaxis()->SetRangeUser(0., 2.);
1128 }
1129 
1130 /*--------------------------------------------------------------------*/
1131 void adjustMaximum(TH1F* histos[], int size)
1132 /*--------------------------------------------------------------------*/
1133 {
1134   float theMax(-999.);
1135   for (int i = 0; i < size; i++) {
1136     if (histos[i]->GetMaximum() > theMax)
1137       theMax = histos[i]->GetMaximum();
1138   }
1139 
1140   for (int i = 0; i < size; i++) {
1141     histos[i]->SetMaximum(theMax * 1.25);
1142   }
1143 }
1144 
1145 /*--------------------------------------------------------------------*/
1146 void setPVResolStyle() {
1147   /*--------------------------------------------------------------------*/
1148 
1149   writeExtraText = true;  // if extra text
1150   lumi_13p6TeV = "pp collisions";
1151   lumi_13TeV = "pp collisions";
1152   lumi_0p9TeV = "pp collisions";
1153   extraText = "Internal";
1154 
1155   TH1::StatOverflows(kTRUE);
1156   gStyle->SetOptTitle(0);
1157   gStyle->SetOptStat("e");
1158   //gStyle->SetPadTopMargin(0.05);
1159   //gStyle->SetPadBottomMargin(0.15);
1160   //gStyle->SetPadLeftMargin(0.17);
1161   //gStyle->SetPadRightMargin(0.02);
1162   gStyle->SetPadBorderMode(0);
1163   gStyle->SetTitleFillColor(10);
1164   gStyle->SetTitleFont(42);
1165   gStyle->SetTitleColor(1);
1166   gStyle->SetTitleTextColor(1);
1167   gStyle->SetTitleFontSize(0.06);
1168   gStyle->SetTitleBorderSize(0);
1169   gStyle->SetStatColor(kWhite);
1170   gStyle->SetStatFont(42);
1171   gStyle->SetStatFontSize(0.05);  ///---> gStyle->SetStatFontSize(0.025);
1172   gStyle->SetStatTextColor(1);
1173   gStyle->SetStatFormat("6.4g");
1174   gStyle->SetStatBorderSize(1);
1175   gStyle->SetPadTickX(1);  // To get tick marks on the opposite side of the frame
1176   gStyle->SetPadTickY(1);
1177   gStyle->SetPadBorderMode(0);
1178   gStyle->SetOptFit(1);
1179   gStyle->SetNdivisions(510);
1180 
1181   // this is the standard palette
1182   const Int_t NRGBs = 5;
1183   const Int_t NCont = 255;
1184 
1185   Double_t stops[NRGBs] = {0.00, 0.34, 0.61, 0.84, 1.00};
1186   Double_t red[NRGBs] = {0.00, 0.00, 0.87, 1.00, 0.51};
1187   Double_t green[NRGBs] = {0.00, 0.81, 1.00, 0.20, 0.00};
1188   Double_t blue[NRGBs] = {0.51, 1.00, 0.12, 0.00, 0.00};
1189   TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
1190   gStyle->SetNumberContours(NCont);
1191 
1192   /*
1193   // try an alternative palette
1194   const Int_t NRGBs = 6;
1195   const Int_t NCont = 999;
1196 
1197   Double_t stops[NRGBs] = { 0.00, 0.1, 0.34, 0.61, 0.84, 1.00 };
1198   Double_t red[NRGBs]   = { 0.99, 0.0, 0.00, 0.87, 1.00, 0.51 };
1199   Double_t green[NRGBs] = { 0.00, 0.0, 0.81, 1.00, 0.20, 0.00 };
1200   Double_t blue[NRGBs]  = { 0.99, 0.0, 1.00, 0.12, 0.00, 0.00 };
1201 
1202   TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
1203   gStyle->SetNumberContours(NCont);
1204 
1205   */
1206 
1207   /*
1208   const Int_t NRGBs = 9;
1209   const Int_t NCont = 255;
1210  
1211   Double_t stops[NRGBs] = { 0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000};
1212 
1213   // dark body radiator
1214   // Double_t red[NRGBs]   = { 0./255., 45./255., 99./255., 156./255., 212./255., 230./255., 237./255., 234./255., 242./255.};
1215   // Double_t green[NRGBs] = { 0./255.,  0./255.,  0./255.,  45./255., 101./255., 168./255., 238./255., 238./255., 243./255.};
1216   // Double_t blue[NRGBs]  = { 0./255.,  1./255.,  1./255.,   3./255.,   9./255.,   8./255.,  11./255.,  95./255., 230./255.};
1217   
1218   // printable on grey
1219   //Double_t red[9]   = { 0./255.,   0./255.,   0./255.,  70./255., 148./255., 231./255., 235./255., 237./255., 244./255.};
1220   //Double_t green[9] = { 0./255.,   0./255.,   0./255.,   0./255.,   0./255.,  69./255.,  67./255., 216./255., 244./255.};
1221   //Double_t blue[9]  = { 0./255., 102./255., 228./255., 231./255., 177./255., 124./255., 137./255.,  20./255., 244./255.};
1222 
1223   // thermometer
1224   //Double_t red[9]   = {  34./255.,  70./255., 129./255., 187./255., 225./255., 226./255., 216./255., 193./255., 179./255.};
1225   //Double_t green[9] = {  48./255.,  91./255., 147./255., 194./255., 226./255., 229./255., 196./255., 110./255.,  12./255.};
1226   //Double_t blue[9]  = { 234./255., 212./255., 216./255., 224./255., 206./255., 110./255.,  53./255.,  40./255.,  29./255.};
1227 
1228   // visible spectrum
1229   Double_t red[9]   = { 18./255.,  72./255.,   5./255.,  23./255.,  29./255., 201./255., 200./255., 98./255., 29./255.};
1230   Double_t green[9] = {  0./255.,   0./255.,  43./255., 167./255., 211./255., 117./255.,   0./255.,  0./255.,  0./255.};
1231   Double_t blue[9]  = { 51./255., 203./255., 177./255.,  26./255.,  10./255.,   9./255.,   8./255.,  3./255.,  0./255.};
1232 
1233   TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
1234   gStyle->SetNumberContours(NCont);
1235   */
1236 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!