Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:56:58 

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