OfflineValidation/macros/MultiRunAndPlotPVValidation.C

0001 #include "TArrow.h"
0002 #include "TAxis.h"
0003 #include "TCanvas.h"
0004 #include "TF1.h"
0005 #include "TFile.h"
0006 #include "TGaxis.h"
0007 #include "TGraph.h"
0008 #include "TGraphAsymmErrors.h"
0009 #include "TGraphErrors.h"
0010 #include "TH1F.h"
0011 #include "TH2F.h"
0012 #include "TLatex.h"
0013 #include "TLegend.h"
0014 #include "TList.h"
0015 #include "TMath.h"
0016 #include "TObjArray.h"
0017 #include "TObjString.h"
0018 #include "TPad.h"
0019 #include "TPaveText.h"
0020 #include "TProcPool.h"
0021 #include "TProfile.h"
0022 #include "TROOT.h"
0023 #include "TStyle.h"
0024 #include "TSystem.h"
0025 #include "TSystemDirectory.h"
0026 #include "TSystemFile.h"
0027 #include <TStopwatch.h>
0028 #include <algorithm>
0029 #include <bitset>
0030 #include <fstream>
0031 #include <functional>
0032 #include <iostream>
0033 #include <iterator>
0034 #include <map>
0035 #include <sstream>
0036 #include <vector>
0037
0038 /*!
0039  * \def some convenient I/O
0040  */
0041 #define logInfo std::cout << "INFO: "
0042 #define logWarning std::cout << "WARNING: "
0043 #define logError std::cout << "ERROR!!! "
0044
0045 /*!
0046  * \def boolean to decide if it is in debug mode
0047  */
0048 #define VERBOSE false
0049
0050 /*!
0051  * \def number of workers
0052  */
0053 const size_t nWorkers = 20;
0054
0055 /*!
0056  * \def basically the y-values of a TGraph
0057  */
0058 typedef std::map<TString, std::vector<double> > alignmentTrend;
0059
0060 namespace pv {
0061   enum view { dxyphi, dzphi, dxyeta, dzeta, pT, generic };
0062
0063   /*! \fn closest
0064    *  \brief method to find first value that doesn not compare left
0065    */
0066
0067   int closest(std::vector<int> const &vec, int value) {
0068     auto const it = std::lower_bound(vec.begin(), vec.end(), value);
0069     if (it == vec.end()) {
0070       return -1;
0071     }
0072     return *it;
0073   }
0074
0075   const Int_t markers[8] = {kFullSquare,
0076                             kFullCircle,
0077                             kFullTriangleDown,
0078                             kOpenSquare,
0079                             kOpenCircle,
0080                             kFullTriangleUp,
0081                             kOpenTriangleDown,
0082                             kOpenTriangleUp};
0083   const Int_t colors[8] = {kBlack, kRed, kBlue, kGreen + 2, kMagenta, kViolet, kCyan, kYellow};
0084
0085   /*! \struct biases
0086    *  \brief Structure biases
0087    *         Contains characterization of a single run PV bias plot
0088    *
0089    * @param m_mean:             mean value of the profile points
0090    * @param m_rms:              RMS value of the profle points
0091    * @param m_w_mean:           mean weighted on the errors
0092    * @param m_w_rms:            RMS weighted on the errors
0093    * @param m_min:              minimum of the profile
0094    * @param m_max:              maximum of the profile
0095    * @param m_chi2:             chi2 of a liner fit
0096    * @param m_ndf:              number of the dof of the linear fit
0097    * @param m_ks:               Kolmogorov-Smirnov score of comparison with flat line
0098    */
0099
0100   struct biases {
0101     // contructor
0102     biases(double mean, double rms, double wmean, double wrms, double min, double max, double chi2, int ndf, double ks) {
0103       m_mean = mean;
0104       m_rms = rms;
0105       m_w_mean = wmean;
0106       m_w_rms = wrms;
0107       m_min = min;
0108       m_max = max;
0109       m_chi2 = chi2;
0110       m_ndf = ndf;
0111       m_ks = ks;
0112     }
0113
0114     // empty constructor
0115     biases() { init(); }
0116
0117     /*! \fn init
0118      *  \brief initialising all members one by one
0119      */
0120
0121     void init() {
0122       m_mean = 0;
0123       m_rms = 0.;
0124       m_min = +999.;
0125       m_max = -999.;
0126       m_w_mean = 0.;
0127       m_w_rms = 0.;
0128       m_chi2 = -1.;
0129       m_ndf = 0.;
0130       m_ks = 9999.;
0131     }
0132
0133     double getMean() { return m_mean; }
0134     double getWeightedMean() { return m_w_mean; }
0135     double getRMS() { return m_rms; }
0136     double getWeightedRMS() { return m_w_rms; }
0137     double getMin() { return m_min; }
0138     double getMax() { return m_max; }
0139     double getChi2() { return m_chi2; }
0140     double getNDF() { return m_ndf; }
0141     double getNormChi2() { return double(m_chi2) / double(m_ndf); }
0142     double getChi2Prob() { return TMath::Prob(m_chi2, m_ndf); }
0143     double getKSScore() { return m_ks; }
0144
0145   private:
0146     double m_mean;
0147     double m_min;
0148     double m_max;
0149     double m_rms;
0150     double m_w_mean;
0151     double m_w_rms;
0152     double m_chi2;
0153     int m_ndf;
0154     double m_ks;
0155   };
0156
0157   /*! \struct wrappedTrends
0158    *  \brief Structure wrappedTrends
0159    *         Contains the ensemble vs run number of the alignmentTrend characterization
0160    *
0161    * @param mean:             alignmentTrend of the mean value of the profile points
0162    * @param low:              alignmentTrend of the lowest value of the profle points
0163    * @param high:             alignmentTrend of the highest value of the profile points
0164    * @param lowerr:           alignmentTrend of the difference between the lowest value and the mean of the profile points
0165    * @param higherr:          alignmentTrend of the difference between the highest value and the mean of the profile points
0166    * @param chi2:             alignmentTrend of the chi2 value of a linear fit to the profile points
0167    * @param KS:               alignmentTrend of the Kolmogorow-Smirnov score of the comarison of the profile points to a flat line
0168    */
0169
0170   struct wrappedTrends {
0171     /*! \fn wrappedTrends
0172      *  \brief Constructor of structure wrappedTrends, initialising all members from DMRs directly (with split)
0173      */
0174
0175     wrappedTrends(alignmentTrend mean,
0176                   alignmentTrend low,
0177                   alignmentTrend high,
0178                   alignmentTrend lowerr,
0179                   alignmentTrend higherr,
0180                   alignmentTrend chi2,
0181                   alignmentTrend KS) {
0182       logInfo << "pv::wrappedTrends c'tor" << std::endl;
0183
0184       m_mean = mean;
0185       m_low = low;
0186       m_high = high;
0187       m_lowerr = lowerr;
0188       m_higherr = higherr;
0189       m_chi2 = chi2;
0190       m_KS = KS;
0191     }
0192
0193     alignmentTrend getMean() const { return m_mean; }
0194     alignmentTrend getLow() const { return m_low; }
0195     alignmentTrend getHigh() const { return m_high; }
0196     alignmentTrend getLowErr() const { return m_lowerr; }
0197     alignmentTrend getHighErr() const { return m_higherr; }
0198     alignmentTrend getChi2() const { return m_chi2; }
0199     alignmentTrend getKS() const { return m_KS; }
0200
0201   private:
0202     alignmentTrend m_mean;
0203     alignmentTrend m_low;
0204     alignmentTrend m_high;
0205     alignmentTrend m_lowerr;
0206     alignmentTrend m_higherr;
0207     alignmentTrend m_chi2;
0208     alignmentTrend m_KS;
0209   };
0210
0211   /*! \struct bundle
0212    *  \brief Structure bundle
0213    *         Contains the ensemble of all the information to build the graphs alignmentTrends
0214    *
0215    * @param nObjects                     int, number of alignments to be considered
0216    * @param dataType                     TString, type of the data to be displayed (time, lumi)
0217    * @param dataTypeLabel                TString, x-axis label
0218    * @param lumiMapByrun                 std::map of the luminoisty by run number
0219    * @param times                        std::map of the date (UTC) by run number
0220    * @param lumiaxisformat               boolean, is the x-axis of the type lumi
0221    * @param timeaxisformat               boolean, is the x-axis of the type time
0222    */
0223
0224   struct bundle {
0225     bundle(int nObjects,
0226            const TString &dataType,
0227            const TString &dataTypeLabel,
0228            const std::map<int, double> &lumiMapByRun,
0229            const std::map<int, TDatime> &times,
0230            const bool &lumiaxisformat,
0231            const bool &timeaxisformat,
0232            const bool &useRMS) {
0233       m_nObjects = nObjects;
0234       m_datatype = dataType.Data();
0235       m_datatypelabel = dataTypeLabel.Data();
0236       m_lumiMapByRun = lumiMapByRun;
0237       m_times = times;
0238       m_useRMS = useRMS;
0239
0240       logInfo << "pv::bundle c'tor: " << dataTypeLabel << " member: " << m_datatypelabel << std::endl;
0241
0242       // make sure you don't use them at the same time
0243       if (lumiaxisformat || timeaxisformat) {
0244         assert(lumiaxisformat != timeaxisformat);
0245       }
0246
0247       if (lumiaxisformat) {
0248         logInfo << "is lumiaxis format" << std::endl;
0249         m_axis_types.set(0);
0250       } else if (timeaxisformat) {
0251         logInfo << "is timeaxis format" << std::endl;
0252         m_axis_types.set(1);
0253       } else {
0254         logInfo << "is runaxis format" << std::endl;
0255       }
0256
0257       logInfo << m_axis_types << std::endl;
0258
0259       m_totalLumi = lumiMapByRun.rbegin()->second;
0260     }
0261
0262     int getNObjects() const { return m_nObjects; }
0263     const char *getDataType() const { return m_datatype; }
0264     const char *getDataTypeLabel() const { return m_datatypelabel; }
0265     const std::map<int, double> getLumiMapByRun() const { return m_lumiMapByRun; }
0266     double getTotalLumi() const { return m_totalLumi; }
0267     const std::map<int, TDatime> getTimes() const { return m_times; }
0268     bool isLumiAxis() const { return m_axis_types.test(0); }
0269     bool isTimeAxis() const { return m_axis_types.test(1); }
0270     bool isUsingRMS() const { return m_useRMS; }
0271     void printAll() {
0272       logInfo << "dataType      " << m_datatype << std::endl;
0273       logInfo << "dataTypeLabel " << m_datatypelabel << std::endl;
0274       if (this->isLumiAxis())
0275         logInfo << "is lumi axis" << std::endl;
0276       if (this->isTimeAxis())
0277         logInfo << "is time axis" << std::endl;
0278     }
0279
0280   private:
0281     int m_nObjects;
0282     const char *m_datatype;
0283     const char *m_datatypelabel;
0284     float m_totalLumi;
0285     std::map<int, double> m_lumiMapByRun;
0286     std::map<int, TDatime> m_times;
0287     std::bitset<2> m_axis_types;
0288     bool m_useRMS;
0289   };
0290
0291 }  // namespace pv
0292
0293 // auxilliary struct to store
0294 // histogram features
0295 struct unrolledHisto {
0296   double m_y_min;
0297   double m_y_max;
0298   unsigned int m_n_bins;
0299   std::vector<double> m_bincontents;
0300
0301   unrolledHisto() {
0302     m_y_min = 0.;
0303     m_y_max = 0.;
0304     m_n_bins = 0;
0305     m_bincontents.clear();
0306   }  // first constructor empty
0307
0308   unrolledHisto(const double &y_min,
0309                 const double &y_max,
0310                 const unsigned int &n_bins,
0311                 const std::vector<double> &bincontents) {
0312     m_y_min = y_min;
0313     m_y_max = y_max;
0314     m_n_bins = n_bins;
0315     m_bincontents = bincontents;
0316   }  //look, a constructor
0317
0318   double get_y_min() { return m_y_min; }
0319
0320   double get_y_max() { return m_y_max; }
0321
0322   unsigned int get_n_bins() { return m_n_bins; }
0323
0324   std::vector<double> get_bin_contents() { return m_bincontents; }
0325
0326   double get_integral() {
0327     double ret(0.);
0328     for (const auto &binc : m_bincontents) {
0329       ret += binc;
0330     }
0331     return ret;
0332   }
0333 };
0334
0335 /*! \struct outTrends
0336    *  \brief Structure outTrends
0337    *         Contains the ensemble of all the alignmentTrends built by the functor
0338    *
0339    * @param m_index                     int, to keep track of which chunk of data has been processed
0340    * @param m_lumiSoFar                 double, luminosity in this section of the data
0341    * @param m_runs                      std::vector, list of the run processed in this section
0342    * @param m_lumiByRun                 std::vector, list of the luminisoties per run, indexed
0343    * @param m_lumiMarpByRun             std::map, map of the luminosities per run
0344    * @param m_dxyPhiMeans               alignmentTrend of the mean values of the profile dxy vs phi
0345    * @param m_dxyPhiChi2                alignmentTrend of chi2 of the linear fit per profile dxy vs phi
0346    * @param m_dxyPhiKS                  alignmentTrend of Kolmogorow-Smirnov score of comparison of dxy vs phi profile with flat line
0347    * @param m_dxyPhiHi                  alignmentTrend of the highest value of the profile dxy vs phi
0348    * @param m_dxyPhiLo                  alignmentTrend of the lowest value of the profile dxy vs phi
0349    * @param m_dxyEtaMeans               alignmentTrend of the mean values of the profile dxy vs eta
0350    * @param m_dxyEtaChi2                alignmentTrend of chi2 of the linear fit per profile dxy vs eta
0351    * @param m_dxyEtaKS                  alignmentTrend of Kolmogorow-Smirnov score of comparison of dxy vs eta profile with flat line
0352    * @param m_dxyEtaHi                  alignmentTrend of the highest value of the profile dxy vs eta
0353    * @param m_dxyEtaLo                  alignmentTrend of the lowest value of the profile dxy vs eta
0354    * @param m_dzPhiMeans                alignmentTrend of the mean values of the profile dz vs phi
0355    * @param m_dzPhiChi2                 alignmentTrend of chi2 of the linear fit per profile dz vs phi
0356    * @param m_dzPhiKS                   alignmentTrend of Kolmogorow-Smirnov score of comparison of dz vs phi profile with flat line
0357    * @param m_dzPhiHi                   alignmentTrend of the highest value of the profile dz vs phi
0358    * @param m_dzPhiLo                   alignmentTrend of the lowest value of the profile dz vs phi
0359    * @param m_dzEtaMeans                alignmentTrend of the mean values of the profile dz vs eta
0360    * @param m_dzEtaChi2                 alignmentTrend of chi2 of the linear fit per profile dz vs eta
0361    * @param m_dzEtaKS                   alignmentTrend of Kolmogorow-Smirnov score of comparison of dz vs eta profile with flat line
0362    * @param m_dzEtaHi                   alignmentTrend of the highest value of the profile dz vs eta
0363    * @param m_dzEtaLo                   alignmentTrend of the lowest value of the profile dz vs eta
0364    * @param m_dxyVect                   map of the unrolled histograms for dxy residuals
0365    * @param m_dzVect                    map of the unrolled histograms for dz residulas
0366    */
0367
0368 struct outTrends {
0369   int m_index;
0370   double m_lumiSoFar;
0371   std::vector<double> m_runs;
0372   std::vector<double> m_lumiByRun;
0373   std::map<int, double> m_lumiMapByRun;
0374   alignmentTrend m_dxyPhiMeans;
0375   alignmentTrend m_dxyPhiChi2;
0376   alignmentTrend m_dxyPhiKS;
0377   alignmentTrend m_dxyPhiHi;
0378   alignmentTrend m_dxyPhiLo;
0379   alignmentTrend m_dxyEtaMeans;
0380   alignmentTrend m_dxyEtaChi2;
0381   alignmentTrend m_dxyEtaKS;
0382   alignmentTrend m_dxyEtaHi;
0383   alignmentTrend m_dxyEtaLo;
0384   alignmentTrend m_dzPhiMeans;
0385   alignmentTrend m_dzPhiChi2;
0386   alignmentTrend m_dzPhiKS;
0387   alignmentTrend m_dzPhiHi;
0388   alignmentTrend m_dzPhiLo;
0389   alignmentTrend m_dzEtaMeans;
0390   alignmentTrend m_dzEtaChi2;
0391   alignmentTrend m_dzEtaKS;
0392   alignmentTrend m_dzEtaHi;
0393   alignmentTrend m_dzEtaLo;
0394   std::map<TString, std::vector<unrolledHisto> > m_dxyVect;
0395   std::map<TString, std::vector<unrolledHisto> > m_dzVect;
0396
0397   void init() {
0398     m_index = -1;
0399     m_lumiSoFar = 0.;
0400     m_runs.clear();
0401     m_lumiByRun.clear();
0402     m_lumiMapByRun.clear();
0403
0404     m_dxyPhiMeans.clear();
0405     m_dxyPhiChi2.clear();
0406     m_dxyPhiKS.clear();
0407     m_dxyPhiHi.clear();
0408     m_dxyPhiLo.clear();
0409
0410     m_dxyEtaMeans.clear();
0411     m_dxyEtaChi2.clear();
0412     m_dxyEtaKS.clear();
0413     m_dxyEtaHi.clear();
0414     m_dxyEtaLo.clear();
0415
0416     m_dzPhiMeans.clear();
0417     m_dzPhiChi2.clear();
0418     m_dzPhiKS.clear();
0419     m_dzPhiHi.clear();
0420     m_dzPhiLo.clear();
0421
0422     m_dzEtaMeans.clear();
0423     m_dzEtaChi2.clear();
0424     m_dzEtaKS.clear();
0425     m_dzEtaHi.clear();
0426     m_dzEtaLo.clear();
0427
0428     m_dxyVect.clear();
0429     m_dzVect.clear();
0430   }
0431 };
0432
0433 // forward declarations
0434 void MultiRunPVValidation(TString namesandlabels = "",
0435                           bool lumi_axis_format = false,
0436                           bool time_axis_format = false,
0437                           bool useRMS = true,
0438                           TString lumiInputFile = "");
0439 outTrends processData(size_t iter,
0440                       std::vector<int> intersection,
0441                       const Int_t nDirs_,
0442                       const char *dirs[10],
0443                       TString LegLabels[10],
0444                       bool useRMS);
0445
0446 void arrangeOutCanvas(TCanvas *canv,
0447                       TH1F *m_11Trend[100],
0448                       TH1F *m_12Trend[100],
0449                       TH1F *m_21Trend[100],
0450                       TH1F *m_22Trend[100],
0451                       Int_t nFiles,
0452                       TString LegLabels[10],
0453                       unsigned int theRun);
0454
0455 pv::biases getBiases(TH1F *hist);
0456 unrolledHisto getUnrolledHisto(TH1F *hist);
0457
0458 TH1F *DrawConstant(TH1F *hist, Int_t iter, Double_t theConst);
0459 TH1F *DrawConstantWithErr(TH1F *hist, Int_t iter, Double_t theConst);
0460 TH1F *DrawConstantGraph(TGraph *graph, Int_t iter, Double_t theConst);
0461 std::vector<int> list_files(const char *dirname = ".", const char *ext = ".root");
0462 TH1F *checkTH1AndReturn(TFile *f, TString address);
0463 void MakeNiceTrendPlotStyle(TH1 *hist, Int_t color);
0464 void cmsPrel(TPad *pad, size_t ipads = 1);
0465 void makeNewXAxis(TH1 *h);
0466 void beautify(TGraph *g);
0467 void beautify(TH1 *h);
0468 void adjustmargins(TCanvas *canv);
0469 void adjustmargins(TVirtualPad *canv);
0470 void setStyle();
0471 pv::view checkTheView(const TString &toCheck);
0472 template <typename T>
0473 void timify(T *mgr);
0474 Double_t getMaximumFromArray(TObjArray *array);
0475 void superImposeIOVBoundaries(TCanvas *c,
0476                               bool lumi_axis_format,
0477                               bool time_axis_format,
0478                               const std::map<int, double> &lumiMapByRun,
0479                               const std::map<int, TDatime> &timeMap,
0480                               bool drawText = true);
0481 void outputGraphs(const pv::wrappedTrends &allInputs,
0482                   const std::vector<double> &ticks,
0483                   const std::vector<double> &ex_ticks,
0484                   TCanvas *&canv,
0485                   TCanvas *&mean_canv,
0486                   TCanvas *&rms_canv,
0487                   TGraph *&g_mean,
0488                   TGraph *&g_chi2,
0489                   TGraph *&g_KS,
0490                   TGraph *&g_low,
0491                   TGraph *&g_high,
0492                   TGraphAsymmErrors *&g_asym,
0493                   TH1F *h_RMS[],
0494                   const pv::bundle &mybundle,
0495                   const pv::view &theView,
0496                   const int index,
0497                   TObjArray *&array,
0498                   const TString &label,
0499                   TLegend *&legend);
0500
0501 /*! \fn split
0502  *  \brief utility function to split strings
0503  */
0504
0505 /*--------------------------------------------------------------------*/
0506 std::vector<std::string> split(const std::string &s, char delimiter)
0507 /*--------------------------------------------------------------------*/
0508 {
0509   std::vector<std::string> tokens;
0510   std::string token;
0511   std::istringstream tokenStream(s);
0512   while (std::getline(tokenStream, token, delimiter)) {
0513     tokens.push_back(token);
0514   }
0515   return tokens;
0516 }
0517
0518 ///////////////////////////////////
0519 //
0520 //  Main function
0521 //
0522 ///////////////////////////////////
0523
0524 void MultiRunPVValidation(
0525     TString namesandlabels, bool lumi_axis_format, bool time_axis_format, bool useRMS, TString lumiInputFile) {
0526   TStopwatch timer;
0527   timer.Start();
0528
0529   using namespace std::placeholders;  // for _1, _2, _3...
0530   gROOT->ProcessLine("gErrorIgnoreLevel = kError;");
0531
0532   ROOT::EnableThreadSafety();
0533   TH1::AddDirectory(kFALSE);
0534
0535   // consistency check, we cannot do plot vs lumi if time_axis
0536   if (lumi_axis_format && time_axis_format) {
0537     logError << "##########################################################################################"
0538              << std::endl;
0539     logError << "msg-i: MultiRunPVValidation(): you're requesting both summary vs lumi and vs time, " << std::endl;
0540     logError << "       this combination is inconsistent --> exiting!" << std::endl;
0541     exit(EXIT_FAILURE);
0542   }
0543
0544   // preload the dates from file
0545   std::map<int, TDatime> times;
0546   if (time_axis_format) {
0547     std::ifstream infile("times.txt");
0548
0549     if (!infile) {
0550       logError << "Required time axis options, but missing input times file :(" << std::endl;
0551       logError << " -- exiting" << std::endl;
0552       exit(EXIT_FAILURE);
0553     }
0554
0555     std::string line;
0556     while (std::getline(infile, line)) {
0557       std::istringstream iss(line);
0558       std::string a, b, c;
0559       if (!(iss >> a >> b >> c)) {
0560         break;
0561       }  // error
0562
0563       //logInfo<<a<<"  "<<b<<"   "<<c<<"   "<<std::endl;
0564
0565       int run = std::stoi(a);
0566       auto tokens_b = split(b, '-');
0567       int year = std::stoi(tokens_b[0]);
0568       int month = std::stoi(tokens_b[1]);
0569       int day = std::stoi(tokens_b[2]);
0570
0571       auto tokens_c = split(c, '.');
0572       auto tokens_c1 = split(tokens_c[0], ':');
0573
0574       int hour = std::stoi(tokens_c1[0]);
0575       int minute = std::stoi(tokens_c1[2]);
0576       int second = std::stoi(tokens_c1[2]);
0577
0578       //logInfo<<run<<" "<<year<<" "<<month<<" "<<day<<" "<<hour<<" "<<minute<<" "<<second<<" "<<std::endl;
0579
0580       TDatime da(year, month, day, hour, minute, second);
0581       times[run] = da;
0582     }
0583   }  // if time axis in the plots
0584
0585   //std::ofstream outfile ("lumiByRun.txt");
0586   std::ofstream outfile("log.txt");
0587   setStyle();
0588
0589   TList *DirList = new TList();
0590   TList *LabelList = new TList();
0591
0592   TObjArray *nameandlabelpairs = namesandlabels.Tokenize(",");
0593   for (Int_t i = 0; i < nameandlabelpairs->GetEntries(); ++i) {
0594     TObjArray *aFileLegPair = TString(nameandlabelpairs->At(i)->GetName()).Tokenize("=");
0595
0596     if (aFileLegPair->GetEntries() == 2) {
0597       DirList->Add(aFileLegPair->At(0));
0598       LabelList->Add(aFileLegPair->At(1));
0599     } else {
0600       logError << "Please give file name and legend entry in the following form:\n"
0601                << " filename1=legendentry1,filename2=legendentry2\n";
0602     }
0603   }
0604
0605   const Int_t nDirs_ = DirList->GetSize();
0606   TString LegLabels[10];
0607   const char *dirs[10];
0608
0609   std::vector<int> intersection;
0610   std::vector<double> runs;
0611   std::vector<double> x_ticks;
0612   std::vector<double> ex_ticks = {0.};
0613
0614   std::vector<double> runtimes;
0615   if (time_axis_format) {
0616     for (const auto &element : times) {
0617       runtimes.push_back((element.second).Convert());
0618     }
0619   }
0620
0621   for (Int_t j = 0; j < nDirs_; j++) {
0622     // Retrieve labels
0623     TObjString *legend = (TObjString *)LabelList->At(j);
0624     TObjString *dir = (TObjString *)DirList->At(j);
0625     LegLabels[j] = legend->String();
0626     dirs[j] = (dir->String()).Data();
0627     logInfo << "MultiRunPVValidation(): label[" << j << "]" << LegLabels[j] << endl;
0628
0629     std::vector<int> currentList = list_files(dirs[j]);
0630     std::vector<int> tempSwap;
0631
0632     std::sort(currentList.begin(), currentList.end());
0633
0634     if (j == 0) {
0635       intersection = currentList;
0636     }
0637
0638     std::sort(intersection.begin(), intersection.end());
0639
0640     std::set_intersection(
0641         currentList.begin(), currentList.end(), intersection.begin(), intersection.end(), std::back_inserter(tempSwap));
0642
0643     intersection.clear();
0644     intersection = tempSwap;
0645     tempSwap.clear();
0646   }
0647
0648   // prelod the lumi from file
0649   std::vector<double> lumiByRun;
0650   std::map<int, double> lumiMapByRun;
0651
0652   bool useLumiByFile = (lumiInputFile.Length() > 0);
0653
0654   if (lumi_axis_format && useLumiByFile) {
0655     std::ifstream lumifile(lumiInputFile.Data());
0656     if (!lumifile) {
0657       logError << "Required luminosity from file, but missing input file :(" << std::endl;
0658       logError << " -- exiting" << std::endl;
0659       exit(EXIT_FAILURE);
0660     }
0661
0662     std::string line;
0663     double lumiSoFar = 0.0;
0664     while (std::getline(lumifile, line)) {
0665       std::istringstream iss(line);
0666       std::string a, b;
0667       if (!(iss >> a >> b)) {
0668         break;
0669       }  // error
0670       int run = std::stoi(a);
0671       double lumi = std::stod(b) / 1000.;
0672
0673       // check if the run is in the list
0674       if (std::find(intersection.begin(), intersection.end(), run) != intersection.end()) {
0675         lumiByRun.push_back(lumiSoFar + lumi);
0676         lumiMapByRun[run] = (lumiSoFar + lumi);
0677       } else {
0678         logWarning << " Run: " << run << " is not found in the intersection" << std::endl;
0679       }
0680       // fill the lumi so far irrespective if the run is in the intersection or not
0681       // this corrects the global scale
0682       lumiSoFar += lumi;
0683
0684       logInfo << run << " ====> lumi so far: " << lumiSoFar << std::endl;
0685     }
0686   }
0687
0688   // debug only
0689   for (UInt_t index = 0; index < intersection.size(); index++) {
0690     logInfo << index << " " << intersection[index] << std::endl;
0691   }
0692
0693   // book the vectors of values
0694   alignmentTrend dxyPhiMeans_;
0695   alignmentTrend dxyPhiChi2_;
0696   alignmentTrend dxyPhiHiErr_;
0697   alignmentTrend dxyPhiLoErr_;
0698   alignmentTrend dxyPhiKS_;
0699   alignmentTrend dxyPhiHi_;
0700   alignmentTrend dxyPhiLo_;
0701
0702   alignmentTrend dxyEtaMeans_;
0703   alignmentTrend dxyEtaChi2_;
0704   alignmentTrend dxyEtaHiErr_;
0705   alignmentTrend dxyEtaLoErr_;
0706   alignmentTrend dxyEtaKS_;
0707   alignmentTrend dxyEtaHi_;
0708   alignmentTrend dxyEtaLo_;
0709
0710   alignmentTrend dzPhiMeans_;
0711   alignmentTrend dzPhiChi2_;
0712   alignmentTrend dzPhiHiErr_;
0713   alignmentTrend dzPhiLoErr_;
0714   alignmentTrend dzPhiKS_;
0715   alignmentTrend dzPhiHi_;
0716   alignmentTrend dzPhiLo_;
0717
0718   alignmentTrend dzEtaMeans_;
0719   alignmentTrend dzEtaChi2_;
0720   alignmentTrend dzEtaHiErr_;
0721   alignmentTrend dzEtaLoErr_;
0722   alignmentTrend dzEtaKS_;
0723   alignmentTrend dzEtaHi_;
0724   alignmentTrend dzEtaLo_;
0725
0726   // unrolled histos
0727
0728   std::map<TString, std::vector<unrolledHisto> > dxyVect;
0729   std::map<TString, std::vector<unrolledHisto> > dzVect;
0730
0731   double lumiSoFar = 0.0;
0732
0733   // loop over the runs in the intersection
0734   //unsigned int last = (DEBUG==true) ? 50 : intersection.size();
0735
0736   //std::function<void()> f_processData = std::bind(processData,1,intersection,nDirs_,dirs,LegLabels,lumiSoFar,runs,lumiByRun,lumiMapByRun,useRMS,dxyPhiMeans_,dxyPhiHi_,dxyPhiLo_,dxyEtaMeans_,dxyEtaHi_,dxyEtaLo_,dzPhiMeans_,dzPhiHi_,dzPhiLo_,dzEtaMeans_,dzEtaHi_,dzEtaLo_,dxyVect,dzVect);
0737
0738   logInfo << " pre do-stuff: " << runs.size() << std::endl;
0739
0740   //we should use std::bind to create a functor and then pass it to the procPool
0741   auto f_processData = std::bind(processData, _1, intersection, nDirs_, dirs, LegLabels, useRMS);
0742
0743   //f_processData(0);
0744   //logInfo<<" post do-stuff: " <<  runs.size() << std::endl;
0745
0746   TProcPool procPool(std::min(nWorkers, intersection.size()));
0747   std::vector<size_t> range(std::min(nWorkers, intersection.size()));
0748   std::iota(range.begin(), range.end(), 0);
0749   //procPool.Map([&f_processData](size_t a) { f_processData(a); },{1,2,3});
0750   auto extracts = procPool.Map(f_processData, range);
0751
0752   // sort the extracts according to the global index
0753   std::sort(extracts.begin(), extracts.end(), [](const outTrends &a, const outTrends &b) -> bool {
0754     return a.m_index < b.m_index;
0755   });
0756
0757   // re-assemble everything together
0758   for (auto extractedTrend : extracts) {
0759     if (!useLumiByFile) {
0760       logInfo << "lumiSoFar: " << lumiSoFar << "/fb" << std::endl;
0761     }
0762
0763     runs.insert(std::end(runs), std::begin(extractedTrend.m_runs), std::end(extractedTrend.m_runs));
0764
0765     // luminosity needs a different treatment
0766     // we need to re-sum the luminosity so far
0767
0768     for (const auto &run : extractedTrend.m_runs) {
0769       if (!useLumiByFile) {
0770         logInfo << run << " ====> lumi so far: " << lumiSoFar + extractedTrend.m_lumiMapByRun[run] << std::endl;
0771         lumiByRun.push_back(lumiSoFar + extractedTrend.m_lumiMapByRun[run]);
0772         lumiMapByRun[run] = (lumiSoFar + extractedTrend.m_lumiMapByRun[run]);
0773       }
0774     }
0775
0776     lumiSoFar += (extractedTrend.m_lumiSoFar / 1000.);
0777
0778     /*
0779       lumiByRun.insert(std::end(lumiByRun), std::begin(extractedTrend.m_lumiByRun), std::end(extractedTrend.m_lumiByRun));
0780       lumiMapByRun.insert(extractedTrend.m_lumiMapByRun.begin(),extractedTrend.m_lumiMapByRun.end());
0781     */
0782
0783     for (const auto &label : LegLabels) {
0784       //******************************//
0785       dxyPhiMeans_[label].insert(std::end(dxyPhiMeans_[label]),
0786                                  std::begin(extractedTrend.m_dxyPhiMeans[label]),
0787                                  std::end(extractedTrend.m_dxyPhiMeans[label]));
0788       dxyPhiChi2_[label].insert(std::end(dxyPhiChi2_[label]),
0789                                 std::begin(extractedTrend.m_dxyPhiChi2[label]),
0790                                 std::end(extractedTrend.m_dxyPhiChi2[label]));
0791       dxyPhiKS_[label].insert(std::end(dxyPhiKS_[label]),
0792                               std::begin(extractedTrend.m_dxyPhiKS[label]),
0793                               std::end(extractedTrend.m_dxyPhiKS[label]));
0794
0795       dxyPhiHi_[label].insert(std::end(dxyPhiHi_[label]),
0796                               std::begin(extractedTrend.m_dxyPhiHi[label]),
0797                               std::end(extractedTrend.m_dxyPhiHi[label]));
0798       dxyPhiLo_[label].insert(std::end(dxyPhiLo_[label]),
0799                               std::begin(extractedTrend.m_dxyPhiLo[label]),
0800                               std::end(extractedTrend.m_dxyPhiLo[label]));
0801
0802       //******************************//
0803       dzPhiMeans_[label].insert(std::end(dzPhiMeans_[label]),
0804                                 std::begin(extractedTrend.m_dzPhiMeans[label]),
0805                                 std::end(extractedTrend.m_dzPhiMeans[label]));
0806       dzPhiChi2_[label].insert(std::end(dzPhiChi2_[label]),
0807                                std::begin(extractedTrend.m_dzPhiChi2[label]),
0808                                std::end(extractedTrend.m_dzPhiChi2[label]));
0809       dzPhiKS_[label].insert(std::end(dzPhiKS_[label]),
0810                              std::begin(extractedTrend.m_dzPhiKS[label]),
0811                              std::end(extractedTrend.m_dzPhiKS[label]));
0812
0813       dzPhiHi_[label].insert(std::end(dzPhiHi_[label]),
0814                              std::begin(extractedTrend.m_dzPhiHi[label]),
0815                              std::end(extractedTrend.m_dzPhiHi[label]));
0816       dzPhiLo_[label].insert(std::end(dzPhiLo_[label]),
0817                              std::begin(extractedTrend.m_dzPhiLo[label]),
0818                              std::end(extractedTrend.m_dzPhiLo[label]));
0819
0820       //******************************//
0821       dxyEtaMeans_[label].insert(std::end(dxyEtaMeans_[label]),
0822                                  std::begin(extractedTrend.m_dxyEtaMeans[label]),
0823                                  std::end(extractedTrend.m_dxyEtaMeans[label]));
0824       dxyEtaChi2_[label].insert(std::end(dxyEtaChi2_[label]),
0825                                 std::begin(extractedTrend.m_dxyEtaChi2[label]),
0826                                 std::end(extractedTrend.m_dxyEtaChi2[label]));
0827       dxyEtaKS_[label].insert(std::end(dxyEtaKS_[label]),
0828                               std::begin(extractedTrend.m_dxyEtaKS[label]),
0829                               std::end(extractedTrend.m_dxyEtaKS[label]));
0830
0831       dxyEtaHi_[label].insert(std::end(dxyEtaHi_[label]),
0832                               std::begin(extractedTrend.m_dxyEtaHi[label]),
0833                               std::end(extractedTrend.m_dxyEtaHi[label]));
0834       dxyEtaLo_[label].insert(std::end(dxyEtaLo_[label]),
0835                               std::begin(extractedTrend.m_dxyEtaLo[label]),
0836                               std::end(extractedTrend.m_dxyEtaLo[label]));
0837
0838       //******************************//
0839       dzEtaMeans_[label].insert(std::end(dzEtaMeans_[label]),
0840                                 std::begin(extractedTrend.m_dzEtaMeans[label]),
0841                                 std::end(extractedTrend.m_dzEtaMeans[label]));
0842       dzEtaChi2_[label].insert(std::end(dzEtaChi2_[label]),
0843                                std::begin(extractedTrend.m_dzEtaChi2[label]),
0844                                std::end(extractedTrend.m_dzEtaChi2[label]));
0845       dzEtaKS_[label].insert(std::end(dzEtaKS_[label]),
0846                              std::begin(extractedTrend.m_dzEtaKS[label]),
0847                              std::end(extractedTrend.m_dzEtaKS[label]));
0848
0849       dzEtaHi_[label].insert(std::end(dzEtaHi_[label]),
0850                              std::begin(extractedTrend.m_dzEtaHi[label]),
0851                              std::end(extractedTrend.m_dzEtaHi[label]));
0852       dzEtaLo_[label].insert(std::end(dzEtaLo_[label]),
0853                              std::begin(extractedTrend.m_dzEtaLo[label]),
0854                              std::end(extractedTrend.m_dzEtaLo[label]));
0855
0856       //******************************//
0857       dxyVect[label].insert(std::end(dxyVect[label]),
0858                             std::begin(extractedTrend.m_dxyVect[label]),
0859                             std::end(extractedTrend.m_dxyVect[label]));
0860       dzVect[label].insert(std::end(dzVect[label]),
0861                            std::begin(extractedTrend.m_dzVect[label]),
0862                            std::end(extractedTrend.m_dzVect[label]));
0863     }
0864   }
0865
0866   // extra vectors for low and high boundaries
0867
0868   for (const auto &label : LegLabels) {
0869     for (unsigned int it = 0; it < dxyPhiMeans_[label].size(); it++) {
0870       dxyPhiHiErr_[label].push_back(std::abs(dxyPhiHi_[label][it] - dxyPhiMeans_[label][it]));
0871       dxyPhiLoErr_[label].push_back(std::abs(dxyPhiLo_[label][it] - dxyPhiMeans_[label][it]));
0872       dxyEtaHiErr_[label].push_back(std::abs(dxyEtaHi_[label][it] - dxyEtaMeans_[label][it]));
0873       dxyEtaLoErr_[label].push_back(std::abs(dxyEtaLo_[label][it] - dxyEtaMeans_[label][it]));
0874
0875       if (VERBOSE) {
0876         logInfo << "label: " << label << " means:" << dxyEtaMeans_[label][it] << " low: " << dxyEtaLo_[label][it]
0877                 << " loErr: " << dxyEtaLoErr_[label][it] << std::endl;
0878       }
0879
0880       dzPhiHiErr_[label].push_back(std::abs(dzPhiHi_[label][it] - dzPhiMeans_[label][it]));
0881       dzPhiLoErr_[label].push_back(std::abs(dzPhiLo_[label][it] - dzPhiMeans_[label][it]));
0882       dzEtaHiErr_[label].push_back(std::abs(dzEtaHi_[label][it] - dzEtaMeans_[label][it]));
0883       dzEtaLoErr_[label].push_back(std::abs(dzEtaLo_[label][it] - dzEtaMeans_[label][it]));
0884     }
0885   }
0886
0887   // just check runs are ordered
0888   //for(const auto &run : runs){
0889   //  logInfo<<" "<< run ;
0890   // }
0891
0892   // main function call
0893   /*
0894     processData(0,intersection,nDirs_,dirs,LegLabels,lumiSoFar,runs,lumiByRun,lumiMapByRun,useRMS,
0895     dxyPhiMeans_,dxyPhiHi_,dxyPhiLo_,
0896     dxyEtaMeans_,dxyEtaHi_,dxyEtaLo_,
0897     dzPhiMeans_,dzPhiHi_,dzPhiLo_,
0898     dzEtaMeans_,dzEtaHi_,dzEtaLo_,
0899     dxyVect,dzVect
0900     );
0901   */
0902
0903   // do the trend-plotting!
0904
0905   TCanvas *c_dxy_phi_vs_run = new TCanvas("c_dxy_phi_vs_run", "dxy(#phi) bias vs run number", 2000, 800);
0906   TCanvas *c_dxy_eta_vs_run = new TCanvas("c_dxy_eta_vs_run", "dxy(#eta) bias vs run number", 2000, 800);
0907   TCanvas *c_dz_phi_vs_run = new TCanvas("c_dz_phi_vs_run", "dz(#phi) bias vs run number", 2000, 800);
0908   TCanvas *c_dz_eta_vs_run = new TCanvas("c_dz_eta_vs_run", "dz(#eta) bias vs run number", 2000, 800);
0909
0910   TCanvas *c_RMS_dxy_phi_vs_run = new TCanvas("c_RMS_dxy_phi_vs_run", "dxy(#phi) bias vs run number", 2000, 800);
0911   TCanvas *c_RMS_dxy_eta_vs_run = new TCanvas("c_RMS_dxy_eta_vs_run", "dxy(#eta) bias vs run number", 2000, 800);
0912   TCanvas *c_RMS_dz_phi_vs_run = new TCanvas("c_RMS_dz_phi_vs_run", "dxy(#phi) bias vs run number", 2000, 800);
0913   TCanvas *c_RMS_dz_eta_vs_run = new TCanvas("c_RMS_dz_eta_vs_run", "dxy(#eta) bias vs run number", 2000, 800);
0914
0915   TCanvas *c_mean_dxy_phi_vs_run = new TCanvas("c_mean_dxy_phi_vs_run", "dxy(#phi) bias vs run number", 2000, 800);
0916   TCanvas *c_mean_dxy_eta_vs_run = new TCanvas("c_mean_dxy_eta_vs_run", "dxy(#eta) bias vs run number", 2000, 800);
0917   TCanvas *c_mean_dz_phi_vs_run = new TCanvas("c_mean_dz_phi_vs_run", "dxy(#phi) bias vs run number", 2000, 800);
0918   TCanvas *c_mean_dz_eta_vs_run = new TCanvas("c_mean_dz_eta_vs_run", "dxy(#eta) bias vs run number", 2000, 800);
0919
0920   TCanvas *Scatter_dxy_vs_run = new TCanvas("Scatter_dxy_vs_run", "dxy bias vs run number", 2000, 800);
0921   Scatter_dxy_vs_run->Divide(1, nDirs_);
0922   TCanvas *Scatter_dz_vs_run = new TCanvas("Scatter_dz_vs_run", "dxy bias vs run number", 2000, 800);
0923   Scatter_dz_vs_run->Divide(1, nDirs_);
0924
0925   TCanvas *c_chisquare_vs_run = new TCanvas("c_chisquare_vs_run", "chi2 of pol0 fit vs run number", 2000, 800);
0926   c_chisquare_vs_run->Divide(2, 2);
0927
0928   TCanvas *c_KSScore_vs_run = new TCanvas("c_KSScore_vs_run", "KS score compatibility to 0 vs run number", 2000, 1000);
0929   c_KSScore_vs_run->Divide(2, 2);
0930
0931   // bias on the mean
0932
0933   TGraph *g_dxy_phi_vs_run[nDirs_];
0934   TGraphAsymmErrors *gerr_dxy_phi_vs_run[nDirs_];
0935   TGraph *g_chi2_dxy_phi_vs_run[nDirs_];
0936   TGraph *g_KS_dxy_phi_vs_run[nDirs_];
0937   TGraph *gprime_dxy_phi_vs_run[nDirs_];
0938   TGraph *g_dxy_phi_hi_vs_run[nDirs_];
0939   TGraph *g_dxy_phi_lo_vs_run[nDirs_];
0940
0941   TGraph *g_dxy_eta_vs_run[nDirs_];
0942   TGraphAsymmErrors *gerr_dxy_eta_vs_run[nDirs_];
0943   TGraph *g_chi2_dxy_eta_vs_run[nDirs_];
0944   TGraph *g_KS_dxy_eta_vs_run[nDirs_];
0945   TGraph *gprime_dxy_eta_vs_run[nDirs_];
0946   TGraph *g_dxy_eta_hi_vs_run[nDirs_];
0947   TGraph *g_dxy_eta_lo_vs_run[nDirs_];
0948
0949   TGraph *g_dz_phi_vs_run[nDirs_];
0950   TGraphAsymmErrors *gerr_dz_phi_vs_run[nDirs_];
0951   TGraph *g_chi2_dz_phi_vs_run[nDirs_];
0952   TGraph *g_KS_dz_phi_vs_run[nDirs_];
0953   TGraph *gprime_dz_phi_vs_run[nDirs_];
0954   TGraph *g_dz_phi_hi_vs_run[nDirs_];
0955   TGraph *g_dz_phi_lo_vs_run[nDirs_];
0956
0957   TGraph *g_dz_eta_vs_run[nDirs_];
0958   TGraphAsymmErrors *gerr_dz_eta_vs_run[nDirs_];
0959   TGraph *g_chi2_dz_eta_vs_run[nDirs_];
0960   TGraph *g_KS_dz_eta_vs_run[nDirs_];
0961   TGraph *gprime_dz_eta_vs_run[nDirs_];
0962   TGraph *g_dz_eta_hi_vs_run[nDirs_];
0963   TGraph *g_dz_eta_lo_vs_run[nDirs_];
0964
0965   // resolutions
0966
0967   TH1F *h_RMS_dxy_phi_vs_run[nDirs_];
0968   TH1F *h_RMS_dxy_eta_vs_run[nDirs_];
0969   TH1F *h_RMS_dz_phi_vs_run[nDirs_];
0970   TH1F *h_RMS_dz_eta_vs_run[nDirs_];
0971
0972   // scatters of integrated bias
0973
0974   TH2F *h2_scatter_dxy_vs_run[nDirs_];
0975   TH2F *h2_scatter_dz_vs_run[nDirs_];
0976
0977   // decide the type
0978
0979   TString theType = "";
0980   TString theTypeLabel = "";
0981   if (lumi_axis_format) {
0982     theType = "luminosity";
0983     theTypeLabel = "Processed luminosity [1/fb]";
0984     if (!useLumiByFile) {
0985       x_ticks = lumiByRun;
0986     } else {
0987       // x_ticks = lumiByRun; (generally it's wrong)
0988       // in case we are passing the luminosity per run by file, need to re-check the map
0989       // in order to fill the x ticks of the graph only for the run which actually pass the selection
0990       double lastLumi(0.);
0991       for (const auto &iRun : runs) {
0992         if (lumiMapByRun.find(iRun) != lumiMapByRun.end()) {
0993           x_ticks.push_back(lumiMapByRun[iRun]);
0994           lastLumi = x_ticks.back();
0995         } else {
0996           // if the run is not find in the lumiMapByRun we just use the lumi of
0997           // the last run analyzed (so we assume that run has lumi=0)
0998           x_ticks.push_back(lastLumi);
0999         }
1000       }
1001     }
1002   } else {
1003     if (!time_axis_format) {
1004       theType = "run number";
1005       theTypeLabel = "run number";
1006       x_ticks = runs;
1007     } else {
1008       theType = "date";
1009       theTypeLabel = "UTC date";
1010       for (const auto &run : runs) {
1011         x_ticks.push_back(times[run].Convert());
1012       }
1013     }
1014   }
1015
1016   //TLegend *my_lego = new TLegend(0.08,0.80,0.18,0.93);
1017   //my_lego-> SetNColumns(2);
1018
1019   TLegend *my_lego = new TLegend(0.75, 0.76, 0.92, 0.93);
1020   //TLegend *my_lego = new TLegend(0.75,0.26,0.92,0.43);
1021   //my_lego->SetHeader("2017 Data","C");
1022   TLine *line = new TLine(0, 0, 1, 0);
1023   line->SetLineColor(kBlue);
1024   line->SetLineStyle(9);
1025   line->SetLineWidth(1);
1026   my_lego->AddEntry(line, "pixel calibration update", "L");
1027   my_lego->SetFillColor(10);
1028   my_lego->SetTextSize(0.042);
1029   my_lego->SetTextFont(42);
1030   my_lego->SetFillColor(10);
1031   my_lego->SetLineColor(10);
1032   my_lego->SetShadowColor(10);
1033
1034   // arrays for storing RMS histograms
1035   TObjArray *arr_dxy_phi = new TObjArray();
1036   TObjArray *arr_dz_phi = new TObjArray();
1037   TObjArray *arr_dxy_eta = new TObjArray();
1038   TObjArray *arr_dz_eta = new TObjArray();
1039
1040   arr_dxy_phi->Expand(nDirs_);
1041   arr_dz_phi->Expand(nDirs_);
1042   arr_dxy_eta->Expand(nDirs_);
1043   arr_dz_eta->Expand(nDirs_);
1044
1045   int ccc = 0;
1046   for (const auto &tick : x_ticks) {
1047     outfile << "***************************************" << std::endl;
1048     outfile << tick << std::endl;
1049     for (Int_t j = 0; j < nDirs_; j++) {
1050       double RMSxyphi = std::abs(dxyPhiLo_[LegLabels[j]][ccc] - dxyPhiHi_[LegLabels[j]][ccc]);
1051       double RMSxyeta = std::abs(dxyEtaLo_[LegLabels[j]][ccc] - dxyEtaHi_[LegLabels[j]][ccc]);
1052       double RMSzphi = std::abs(dzPhiLo_[LegLabels[j]][ccc] - dzPhiHi_[LegLabels[j]][ccc]);
1053       double RMSzeta = std::abs(dzEtaLo_[LegLabels[j]][ccc] - dzEtaHi_[LegLabels[j]][ccc]);
1054
1055       if (RMSxyphi > 100 || RMSxyeta > 100 || RMSzphi > 100 || RMSzeta > 100) {
1056         outfile << LegLabels[j] << " dxy(phi) " << dxyPhiMeans_[LegLabels[j]][ccc] << " "
1057                 << dxyPhiLo_[LegLabels[j]][ccc] << " " << dxyPhiHi_[LegLabels[j]][ccc] << " "
1058                 << std::abs(dxyPhiLo_[LegLabels[j]][ccc] - dxyPhiHi_[LegLabels[j]][ccc]) << std::endl;
1059         outfile << LegLabels[j] << " dxy(eta) " << dxyEtaMeans_[LegLabels[j]][ccc] << " "
1060                 << dxyEtaLo_[LegLabels[j]][ccc] << " " << dxyEtaHi_[LegLabels[j]][ccc] << " "
1061                 << std::abs(dxyEtaLo_[LegLabels[j]][ccc] - dxyEtaHi_[LegLabels[j]][ccc]) << std::endl;
1062         outfile << LegLabels[j] << " dz (phi) " << dzPhiMeans_[LegLabels[j]][ccc] << " " << dzPhiLo_[LegLabels[j]][ccc]
1063                 << " " << dzPhiHi_[LegLabels[j]][ccc] << " "
1064                 << std::abs(dzPhiLo_[LegLabels[j]][ccc] - dzPhiHi_[LegLabels[j]][ccc]) << std::endl;
1065         outfile << LegLabels[j] << " dz (eta) " << dzEtaMeans_[LegLabels[j]][ccc] << " " << dzEtaLo_[LegLabels[j]][ccc]
1066                 << " " << dzEtaHi_[LegLabels[j]][ccc] << " "
1067                 << std::abs(dzEtaLo_[LegLabels[j]][ccc] - dzEtaHi_[LegLabels[j]][ccc]) << std::endl;
1068       }
1069     }
1070     ccc++;
1071   }
1072   outfile << std::endl;
1073
1074   pv::bundle theBundle =
1075       pv::bundle(nDirs_, theType, theTypeLabel, lumiMapByRun, times, lumi_axis_format, time_axis_format, useRMS);
1076   theBundle.printAll();
1077
1078   for (Int_t j = 0; j < nDirs_; j++) {
1079     // check on the sanity
1080     logInfo << "x_ticks.size()= " << x_ticks.size() << " dxyPhiMeans_[LegLabels[" << j
1081             << "]].size()=" << dxyPhiMeans_[LegLabels[j]].size() << std::endl;
1082
1083     // otherwise something very bad has happened
1084     assert(x_ticks.size() == dxyPhiMeans_[LegLabels[j]].size());
1085
1086     // *************************************
1087     // dxy vs phi
1088     // *************************************
1089
1090     auto dxyPhiInputs =
1091         pv::wrappedTrends(dxyPhiMeans_, dxyPhiLo_, dxyPhiHi_, dxyPhiLoErr_, dxyPhiHiErr_, dxyPhiChi2_, dxyPhiKS_);
1092
1093     outputGraphs(dxyPhiInputs,
1094                  x_ticks,
1095                  ex_ticks,
1096                  c_dxy_phi_vs_run,
1097                  c_mean_dxy_phi_vs_run,
1098                  c_RMS_dxy_phi_vs_run,
1099                  g_dxy_phi_vs_run[j],
1100                  g_chi2_dxy_phi_vs_run[j],
1101                  g_KS_dxy_phi_vs_run[j],
1102                  g_dxy_phi_lo_vs_run[j],
1103                  g_dxy_phi_hi_vs_run[j],
1104                  gerr_dxy_phi_vs_run[j],
1105                  h_RMS_dxy_phi_vs_run,
1106                  theBundle,
1107                  pv::dxyphi,
1108                  j,
1109                  arr_dxy_phi,
1110                  LegLabels[j],
1111                  my_lego);
1112
1113     // *************************************
1114     // dxy vs eta
1115     // *************************************
1116
1117     auto dxyEtaInputs =
1118         pv::wrappedTrends(dxyEtaMeans_, dxyEtaLo_, dxyEtaHi_, dxyEtaLoErr_, dxyEtaHiErr_, dxyEtaChi2_, dxyEtaKS_);
1119
1120     outputGraphs(dxyEtaInputs,
1121                  x_ticks,
1122                  ex_ticks,
1123                  c_dxy_eta_vs_run,
1124                  c_mean_dxy_eta_vs_run,
1125                  c_RMS_dxy_eta_vs_run,
1126                  g_dxy_eta_vs_run[j],
1127                  g_chi2_dxy_eta_vs_run[j],
1128                  g_KS_dxy_eta_vs_run[j],
1129                  g_dxy_eta_lo_vs_run[j],
1130                  g_dxy_eta_hi_vs_run[j],
1131                  gerr_dxy_eta_vs_run[j],
1132                  h_RMS_dxy_eta_vs_run,
1133                  theBundle,
1134                  pv::dxyeta,
1135                  j,
1136                  arr_dxy_eta,
1137                  LegLabels[j],
1138                  my_lego);
1139
1140     // *************************************
1141     // dz vs phi
1142     // *************************************
1143
1144     auto dzPhiInputs =
1145         pv::wrappedTrends(dzPhiMeans_, dzPhiLo_, dzPhiHi_, dzPhiLoErr_, dzPhiHiErr_, dzPhiChi2_, dzPhiKS_);
1146
1147     outputGraphs(dzPhiInputs,
1148                  x_ticks,
1149                  ex_ticks,
1150                  c_dz_phi_vs_run,
1151                  c_mean_dz_phi_vs_run,
1152                  c_RMS_dz_phi_vs_run,
1153                  g_dz_phi_vs_run[j],
1154                  g_chi2_dz_phi_vs_run[j],
1155                  g_KS_dz_phi_vs_run[j],
1156                  g_dz_phi_lo_vs_run[j],
1157                  g_dz_phi_hi_vs_run[j],
1158                  gerr_dz_phi_vs_run[j],
1159                  h_RMS_dz_phi_vs_run,
1160                  theBundle,
1161                  pv::dzphi,
1162                  j,
1163                  arr_dz_phi,
1164                  LegLabels[j],
1165                  my_lego);
1166
1167     // *************************************
1168     // dz vs eta
1169     // *************************************
1170
1171     auto dzEtaInputs =
1172         pv::wrappedTrends(dzEtaMeans_, dzEtaLo_, dzEtaHi_, dzEtaLoErr_, dzEtaHiErr_, dzEtaChi2_, dzEtaKS_);
1173
1174     outputGraphs(dzEtaInputs,
1175                  x_ticks,
1176                  ex_ticks,
1177                  c_dz_eta_vs_run,
1178                  c_mean_dz_eta_vs_run,
1179                  c_RMS_dz_eta_vs_run,
1180                  g_dz_eta_vs_run[j],
1181                  g_chi2_dz_eta_vs_run[j],
1182                  g_KS_dz_eta_vs_run[j],
1183                  g_dz_eta_lo_vs_run[j],
1184                  g_dz_eta_hi_vs_run[j],
1185                  gerr_dz_eta_vs_run[j],
1186                  h_RMS_dz_eta_vs_run,
1187                  theBundle,
1188                  pv::dzeta,
1189                  j,
1190                  arr_dz_eta,
1191                  LegLabels[j],
1192                  my_lego);
1193
1194     // *************************************
1195     // Integrated bias dxy scatter plots
1196     // *************************************
1197
1198     h2_scatter_dxy_vs_run[j] =
1199         new TH2F(Form("h2_scatter_dxy_%s", LegLabels[j].Data()),
1200                  Form("scatter of d_{xy} vs %s;%s;d_{xy} [cm]", theType.Data(), theTypeLabel.Data()),
1201                  x_ticks.size() - 1,
1202                  &(x_ticks[0]),
1203                  dxyVect[LegLabels[j]][0].get_n_bins(),
1204                  dxyVect[LegLabels[j]][0].get_y_min(),
1205                  dxyVect[LegLabels[j]][0].get_y_max());
1206     h2_scatter_dxy_vs_run[j]->SetStats(kFALSE);
1207
1208     for (unsigned int runindex = 0; runindex < x_ticks.size(); runindex++) {
1209       for (unsigned int binindex = 0; binindex < dxyVect[LegLabels[j]][runindex].get_n_bins(); binindex++) {
1210         h2_scatter_dxy_vs_run[j]->SetBinContent(runindex + 1,
1211                                                 binindex + 1,
1212                                                 dxyVect[LegLabels[j]][runindex].get_bin_contents().at(binindex) /
1213                                                     dxyVect[LegLabels[j]][runindex].get_integral());
1214       }
1215     }
1216
1217     //Scatter_dxy_vs_run->cd();
1218     h2_scatter_dxy_vs_run[j]->SetFillColorAlpha(pv::colors[j], 0.3);
1219     h2_scatter_dxy_vs_run[j]->SetMarkerColor(pv::colors[j]);
1220     h2_scatter_dxy_vs_run[j]->SetLineColor(pv::colors[j]);
1221     h2_scatter_dxy_vs_run[j]->SetMarkerStyle(pv::markers[j]);
1222
1223     auto h_dxypfx_tmp = (TProfile *)(((TH2F *)h2_scatter_dxy_vs_run[j])->ProfileX(Form("_apfx_%i", j), 1, -1, "o"));
1224     h_dxypfx_tmp->SetName(TString(h2_scatter_dxy_vs_run[j]->GetName()) + "_pfx");
1225     h_dxypfx_tmp->SetStats(kFALSE);
1226     h_dxypfx_tmp->SetMarkerColor(pv::colors[j]);
1227     h_dxypfx_tmp->SetLineColor(pv::colors[j]);
1228     h_dxypfx_tmp->SetLineWidth(2);
1229     h_dxypfx_tmp->SetMarkerSize(1);
1230     h_dxypfx_tmp->SetMarkerStyle(pv::markers[j]);
1231
1232     beautify(h2_scatter_dxy_vs_run[j]);
1233     beautify(h_dxypfx_tmp);
1234
1235     Scatter_dxy_vs_run->cd(j + 1);
1236     adjustmargins(Scatter_dxy_vs_run->cd(j + 1));
1237     //h_dxypfx_tmp->GetYaxis()->SetRangeUser(-0.01,0.01);
1238     //h2_scatter_dxy_vs_run[j]->GetYaxis()->SetRangeUser(-0.5,0.5);
1239     h2_scatter_dxy_vs_run[j]->Draw("colz");
1240     h_dxypfx_tmp->Draw("same");
1241
1242     // *************************************
1243     // Integrated bias dz scatter plots
1244     // *************************************
1245
1246     h2_scatter_dz_vs_run[j] =
1247         new TH2F(Form("h2_scatter_dz_%s", LegLabels[j].Data()),
1248                  Form("scatter of d_{z} vs %s;%s;d_{z} [cm]", theType.Data(), theTypeLabel.Data()),
1249                  x_ticks.size() - 1,
1250                  &(x_ticks[0]),
1251                  dzVect[LegLabels[j]][0].get_n_bins(),
1252                  dzVect[LegLabels[j]][0].get_y_min(),
1253                  dzVect[LegLabels[j]][0].get_y_max());
1254     h2_scatter_dz_vs_run[j]->SetStats(kFALSE);
1255
1256     for (unsigned int runindex = 0; runindex < x_ticks.size(); runindex++) {
1257       for (unsigned int binindex = 0; binindex < dzVect[LegLabels[j]][runindex].get_n_bins(); binindex++) {
1258         h2_scatter_dz_vs_run[j]->SetBinContent(runindex + 1,
1259                                                binindex + 1,
1260                                                dzVect[LegLabels[j]][runindex].get_bin_contents().at(binindex) /
1261                                                    dzVect[LegLabels[j]][runindex].get_integral());
1262       }
1263     }
1264
1265     //Scatter_dz_vs_run->cd();
1266     h2_scatter_dz_vs_run[j]->SetFillColorAlpha(pv::colors[j], 0.3);
1267     h2_scatter_dz_vs_run[j]->SetMarkerColor(pv::colors[j]);
1268     h2_scatter_dz_vs_run[j]->SetLineColor(pv::colors[j]);
1269     h2_scatter_dz_vs_run[j]->SetMarkerStyle(pv::markers[j]);
1270
1271     auto h_dzpfx_tmp = (TProfile *)(((TH2F *)h2_scatter_dz_vs_run[j])->ProfileX(Form("_apfx_%i", j), 1, -1, "o"));
1272     h_dzpfx_tmp->SetName(TString(h2_scatter_dz_vs_run[j]->GetName()) + "_pfx");
1273     h_dzpfx_tmp->SetStats(kFALSE);
1274     h_dzpfx_tmp->SetMarkerColor(pv::colors[j]);
1275     h_dzpfx_tmp->SetLineColor(pv::colors[j]);
1276     h_dzpfx_tmp->SetLineWidth(2);
1277     h_dzpfx_tmp->SetMarkerSize(1);
1278     h_dzpfx_tmp->SetMarkerStyle(pv::markers[j]);
1279
1280     beautify(h2_scatter_dz_vs_run[j]);
1281     beautify(h_dzpfx_tmp);
1282
1283     Scatter_dz_vs_run->cd(j + 1);
1284     adjustmargins(Scatter_dz_vs_run->cd(j + 1));
1285     //h_dzpfx_tmp->GetYaxis()->SetRangeUser(-0.01,0.01);
1286     //h2_scatter_dz_vs_run[j]->GetYaxis()->SetRangeUser(-0.5,0.5);
1287     h2_scatter_dz_vs_run[j]->Draw("colz");
1288     h_dzpfx_tmp->Draw("same");
1289
1290     // ****************************************
1291     // Canvas for chi2 goodness of pol0 fit
1292     // ****************************************
1293
1294     // 1st pad
1295     c_chisquare_vs_run->cd(1);
1296     adjustmargins(c_chisquare_vs_run->cd(1));
1297     g_chi2_dxy_phi_vs_run[j]->SetMarkerStyle(pv::markers[j]);
1298     g_chi2_dxy_phi_vs_run[j]->SetMarkerColor(pv::colors[j]);
1299     g_chi2_dxy_phi_vs_run[j]->SetLineColor(pv::colors[j]);
1300
1301     g_chi2_dxy_phi_vs_run[j]->SetName(Form("g_chi2_dxy_phi_%s", LegLabels[j].Data()));
1302     g_chi2_dxy_phi_vs_run[j]->SetTitle(Form("log_{10}(#chi2/ndf) of d_{xy}(#varphi) fit vs %s", theType.Data()));
1303     g_chi2_dxy_phi_vs_run[j]->GetXaxis()->SetTitle(theTypeLabel.Data());
1304     g_chi2_dxy_phi_vs_run[j]->GetYaxis()->SetTitle("log_{10}(#chi^{2}/ndf) of d_{xy}(#phi) pol0 fit");
1305     g_chi2_dxy_phi_vs_run[j]->GetYaxis()->SetRangeUser(-0.5, 4.5);
1306     if (lumi_axis_format) {
1307       g_chi2_dxy_phi_vs_run[j]->GetXaxis()->SetRangeUser(0., theBundle.getTotalLumi());
1308     }
1309     beautify(g_chi2_dxy_phi_vs_run[j]);
1310     //g_chi2_dxy_phi_vs_run[j]->GetYaxis()->SetTitleOffset(1.3);
1311
1312     if (j == 0) {
1313       g_chi2_dxy_phi_vs_run[j]->Draw("APL");
1314     } else {
1315       g_chi2_dxy_phi_vs_run[j]->Draw("PLsame");
1316     }
1317
1318     if (time_axis_format) {
1319       timify(g_chi2_dxy_phi_vs_run[j]);
1320     }
1321
1322     if (j == nDirs_ - 1) {
1323       my_lego->Draw("same");
1324     }
1325
1326     auto current_pad = static_cast<TCanvas *>(gPad);
1327     superImposeIOVBoundaries(current_pad, lumi_axis_format, time_axis_format, lumiMapByRun, times, false);
1328
1329     // 2nd pad
1330     c_chisquare_vs_run->cd(2);
1331     adjustmargins(c_chisquare_vs_run->cd(2));
1332     g_chi2_dxy_eta_vs_run[j]->SetMarkerStyle(pv::markers[j]);
1333     g_chi2_dxy_eta_vs_run[j]->SetMarkerColor(pv::colors[j]);
1334     g_chi2_dxy_eta_vs_run[j]->SetLineColor(pv::colors[j]);
1335
1336     g_chi2_dxy_eta_vs_run[j]->SetName(Form("g_chi2_dxy_eta_%s", LegLabels[j].Data()));
1337     g_chi2_dxy_eta_vs_run[j]->SetTitle(Form("log_{10}(#chi2/ndf) of d_{xy}(#eta) fit vs %s", theType.Data()));
1338     g_chi2_dxy_eta_vs_run[j]->GetXaxis()->SetTitle(theTypeLabel.Data());
1339     g_chi2_dxy_eta_vs_run[j]->GetYaxis()->SetTitle("log_{10}(#chi^{2}/ndf) of d_{xy}(#eta) pol0 fit");
1340     g_chi2_dxy_eta_vs_run[j]->GetYaxis()->SetRangeUser(-0.5, 4.5);
1341     if (lumi_axis_format) {
1342       g_chi2_dxy_eta_vs_run[j]->GetXaxis()->SetRangeUser(0., theBundle.getTotalLumi());
1343     }
1344     beautify(g_chi2_dxy_eta_vs_run[j]);
1345     //g_chi2_dxy_eta_vs_run[j]->GetYaxis()->SetTitleOffset(1.3);
1346
1347     if (j == 0) {
1348       g_chi2_dxy_eta_vs_run[j]->Draw("APL");
1349     } else {
1350       g_chi2_dxy_eta_vs_run[j]->Draw("PLsame");
1351     }
1352
1353     if (time_axis_format) {
1354       timify(g_chi2_dxy_eta_vs_run[j]);
1355     }
1356
1357     if (j == nDirs_ - 1) {
1358       my_lego->Draw("same");
1359     }
1360
1361     current_pad = static_cast<TCanvas *>(gPad);
1362     superImposeIOVBoundaries(current_pad, lumi_axis_format, time_axis_format, lumiMapByRun, times, false);
1363
1364     //3d pad
1365     c_chisquare_vs_run->cd(3);
1366     adjustmargins(c_chisquare_vs_run->cd(3));
1367     g_chi2_dz_phi_vs_run[j]->SetMarkerStyle(pv::markers[j]);
1368     g_chi2_dz_phi_vs_run[j]->SetMarkerColor(pv::colors[j]);
1369     g_chi2_dz_phi_vs_run[j]->SetLineColor(pv::colors[j]);
1370
1371     g_chi2_dz_phi_vs_run[j]->SetName(Form("g_chi2_dz_phi_%s", LegLabels[j].Data()));
1372     g_chi2_dz_phi_vs_run[j]->SetTitle(Form("log_{10}(#chi2/ndf) of d_{z}(#varphi) fit vs %s", theType.Data()));
1373     g_chi2_dz_phi_vs_run[j]->GetXaxis()->SetTitle(theTypeLabel.Data());
1374     g_chi2_dz_phi_vs_run[j]->GetYaxis()->SetTitle("log_{10}(#chi^{2}/ndf) of d_{z}(#phi) pol0 fit");
1375     g_chi2_dz_phi_vs_run[j]->GetYaxis()->SetRangeUser(-0.5, 4.5);
1376     if (lumi_axis_format) {
1377       g_chi2_dz_phi_vs_run[j]->GetXaxis()->SetRangeUser(0., theBundle.getTotalLumi());
1378     }
1379     beautify(g_chi2_dz_phi_vs_run[j]);
1380     //g_chi2_dz_phi_vs_run[j]->GetYaxis()->SetTitleOffset(1.3);
1381
1382     if (j == 0) {
1383       g_chi2_dz_phi_vs_run[j]->Draw("APL");
1384     } else {
1385       g_chi2_dz_phi_vs_run[j]->Draw("PLsame");
1386     }
1387
1388     if (time_axis_format) {
1389       timify(g_chi2_dz_phi_vs_run[j]);
1390     }
1391
1392     if (j == nDirs_ - 1) {
1393       my_lego->Draw("same");
1394     }
1395
1396     current_pad = static_cast<TCanvas *>(gPad);
1397     superImposeIOVBoundaries(current_pad, lumi_axis_format, time_axis_format, lumiMapByRun, times, false);
1398
1399     //4th pad
1400     c_chisquare_vs_run->cd(4);
1401     adjustmargins(c_chisquare_vs_run->cd(4));
1402     g_chi2_dz_eta_vs_run[j]->SetMarkerStyle(pv::markers[j]);
1403     g_chi2_dz_eta_vs_run[j]->SetMarkerColor(pv::colors[j]);
1404     g_chi2_dz_eta_vs_run[j]->SetLineColor(pv::colors[j]);
1405
1406     g_chi2_dz_eta_vs_run[j]->SetName(Form("g_chi2_dz_eta_%s", LegLabels[j].Data()));
1407     g_chi2_dz_eta_vs_run[j]->SetTitle(Form("log_{10}(#chi2/ndf) of d_{z}(#eta) fit vs %s", theType.Data()));
1408     g_chi2_dz_eta_vs_run[j]->GetXaxis()->SetTitle(theTypeLabel.Data());
1409     g_chi2_dz_eta_vs_run[j]->GetYaxis()->SetTitle("log_{10}(#chi^{2}/ndf) of d_{z}(#eta) pol0 fit");
1410     g_chi2_dz_eta_vs_run[j]->GetYaxis()->SetRangeUser(-0.5, 4.5);
1411     if (lumi_axis_format) {
1412       g_chi2_dz_eta_vs_run[j]->GetXaxis()->SetRangeUser(0., theBundle.getTotalLumi());
1413     }
1414     beautify(g_chi2_dz_eta_vs_run[j]);
1415     //g_chi2_dz_eta_vs_run[j]->GetYaxis()->SetTitleOffset(1.3);
1416
1417     if (j == 0) {
1418       g_chi2_dz_eta_vs_run[j]->Draw("APL");
1419     } else {
1420       g_chi2_dz_eta_vs_run[j]->Draw("PLsame");
1421     }
1422
1423     if (time_axis_format) {
1424       timify(g_chi2_dz_eta_vs_run[j]);
1425     }
1426
1427     if (j == nDirs_ - 1) {
1428       my_lego->Draw("same");
1429     }
1430
1431     current_pad = static_cast<TCanvas *>(gPad);
1432     superImposeIOVBoundaries(current_pad, lumi_axis_format, time_axis_format, lumiMapByRun, times, false);
1433
1434     // ****************************************
1435     // Canvas for Kolmogorov-Smirnov test
1436     // ****************************************
1437
1438     // 1st pad
1439     c_KSScore_vs_run->cd(1);
1440     adjustmargins(c_KSScore_vs_run->cd(1));
1441     g_KS_dxy_phi_vs_run[j]->SetMarkerStyle(pv::markers[j]);
1442     g_KS_dxy_phi_vs_run[j]->SetMarkerColor(pv::colors[j]);
1443     g_KS_dxy_phi_vs_run[j]->SetLineColor(pv::colors[j]);
1444
1445     g_KS_dxy_phi_vs_run[j]->SetName(Form("g_KS_dxy_phi_%s", LegLabels[j].Data()));
1446     g_KS_dxy_phi_vs_run[j]->SetTitle(Form("log_{10}(KS-score) of d_{xy}(#varphi) vs %s", theType.Data()));
1447     g_KS_dxy_phi_vs_run[j]->GetXaxis()->SetTitle(theTypeLabel.Data());
1448     g_KS_dxy_phi_vs_run[j]->GetYaxis()->SetTitle("log_{10}(KS-score) of d_{xy}(#phi) w.r.t 0");
1449     g_KS_dxy_phi_vs_run[j]->GetYaxis()->SetRangeUser(-20., 1.);
1450     beautify(g_KS_dxy_phi_vs_run[j]);
1451     //g_KS_dxy_phi_vs_run[j]->GetYaxis()->SetTitleOffset(1.3);
1452
1453     if (j == 0) {
1454       g_KS_dxy_phi_vs_run[j]->Draw("AP");
1455     } else {
1456       g_KS_dxy_phi_vs_run[j]->Draw("Psame");
1457     }
1458
1459     if (time_axis_format) {
1460       timify(g_KS_dxy_phi_vs_run[j]);
1461     }
1462
1463     if (j == nDirs_ - 1) {
1464       my_lego->Draw("same");
1465     }
1466
1467     // 2nd pad
1468     c_KSScore_vs_run->cd(2);
1469     adjustmargins(c_KSScore_vs_run->cd(2));
1470     g_KS_dxy_eta_vs_run[j]->SetMarkerStyle(pv::markers[j]);
1471     g_KS_dxy_eta_vs_run[j]->SetMarkerColor(pv::colors[j]);
1472     g_KS_dxy_eta_vs_run[j]->SetLineColor(pv::colors[j]);
1473
1474     g_KS_dxy_eta_vs_run[j]->SetName(Form("g_KS_dxy_eta_%s", LegLabels[j].Data()));
1475     g_KS_dxy_eta_vs_run[j]->SetTitle(Form("log_{10}(KS-score) of d_{xy}(#eta) vs %s", theType.Data()));
1476     g_KS_dxy_eta_vs_run[j]->GetXaxis()->SetTitle(theTypeLabel.Data());
1477     g_KS_dxy_eta_vs_run[j]->GetYaxis()->SetTitle("log_{10}(KS-score) of d_{xy}(#eta) w.r.t 0");
1478     g_KS_dxy_eta_vs_run[j]->GetYaxis()->SetRangeUser(-20., 1.);
1479     beautify(g_KS_dxy_eta_vs_run[j]);
1480     //g_KS_dxy_eta_vs_run[j]->GetYaxis()->SetTitleOffset(1.3);
1481
1482     if (j == 0) {
1483       g_KS_dxy_eta_vs_run[j]->Draw("AP");
1484     } else {
1485       g_KS_dxy_eta_vs_run[j]->Draw("Psame");
1486     }
1487
1488     if (time_axis_format) {
1489       timify(g_KS_dxy_eta_vs_run[j]);
1490     }
1491
1492     if (j == nDirs_ - 1) {
1493       my_lego->Draw("same");
1494     }
1495
1496     //3d pad
1497     c_KSScore_vs_run->cd(3);
1498     adjustmargins(c_KSScore_vs_run->cd(3));
1499     g_KS_dz_phi_vs_run[j]->SetMarkerStyle(pv::markers[j]);
1500     g_KS_dz_phi_vs_run[j]->SetMarkerColor(pv::colors[j]);
1501     g_KS_dz_phi_vs_run[j]->SetLineColor(pv::colors[j]);
1502
1503     g_KS_dz_phi_vs_run[j]->SetName(Form("g_KS_dz_phi_%s", LegLabels[j].Data()));
1504     g_KS_dz_phi_vs_run[j]->SetTitle(Form("log_{10}(KS-score) of d_{z}(#varphi) vs %s", theType.Data()));
1505     g_KS_dz_phi_vs_run[j]->GetXaxis()->SetTitle(theTypeLabel.Data());
1506     g_KS_dz_phi_vs_run[j]->GetYaxis()->SetTitle("log_{10}(KS-score) of d_{z}(#phi) w.r.t 0");
1507     g_KS_dz_phi_vs_run[j]->GetYaxis()->SetRangeUser(-20., 1.);
1508     beautify(g_KS_dz_phi_vs_run[j]);
1509     //g_KS_dz_phi_vs_run[j]->GetYaxis()->SetTitleOffset(1.3);
1510
1511     if (j == 0) {
1512       g_KS_dz_phi_vs_run[j]->Draw("AP");
1513     } else {
1514       g_KS_dz_phi_vs_run[j]->Draw("Psame");
1515     }
1516
1517     if (time_axis_format) {
1518       timify(g_KS_dz_phi_vs_run[j]);
1519     }
1520
1521     if (j == nDirs_ - 1) {
1522       my_lego->Draw("same");
1523     }
1524
1525     //4th pad
1526     c_KSScore_vs_run->cd(4);
1527     adjustmargins(c_KSScore_vs_run->cd(4));
1528     g_KS_dz_eta_vs_run[j]->SetMarkerStyle(pv::markers[j]);
1529     g_KS_dz_eta_vs_run[j]->SetMarkerColor(pv::colors[j]);
1530     g_KS_dz_eta_vs_run[j]->SetLineColor(pv::colors[j]);
1531
1532     g_KS_dz_eta_vs_run[j]->SetName(Form("g_KS_dz_eta_%s", LegLabels[j].Data()));
1533     g_KS_dz_eta_vs_run[j]->SetTitle(Form("log_{10}(KS-score) of d_{z}(#eta) vs %s", theType.Data()));
1534     g_KS_dz_eta_vs_run[j]->GetXaxis()->SetTitle(theTypeLabel.Data());
1535     g_KS_dz_eta_vs_run[j]->GetYaxis()->SetTitle("log_{10}(KS-score) of d_{z}(#eta) w.r.t 0");
1536     g_KS_dz_eta_vs_run[j]->GetYaxis()->SetRangeUser(-20., 1.);
1537     beautify(g_KS_dz_eta_vs_run[j]);
1538     //g_KS_dz_eta_vs_run[j]->GetYaxis()->SetTitleOffset(1.3);
1539
1540     if (j == 0) {
1541       g_KS_dz_eta_vs_run[j]->Draw("AP");
1542     } else {
1543       g_KS_dz_eta_vs_run[j]->Draw("Psame");
1544     }
1545
1546     if (time_axis_format) {
1547       timify(g_KS_dz_eta_vs_run[j]);
1548     }
1549
1550     if (j == nDirs_ - 1) {
1551       my_lego->Draw("same");
1552     }
1553   }
1554
1555   // delete the array for the maxima
1556   delete arr_dxy_phi;
1557   delete arr_dz_phi;
1558   delete arr_dxy_eta;
1559   delete arr_dz_eta;
1560
1561   TString append;
1562   if (lumi_axis_format) {
1563     append = "lumi";
1564   } else {
1565     if (time_axis_format) {
1566       append = "date";
1567     } else {
1568       append = "run";
1569     }
1570   }
1571
1572   c_dxy_phi_vs_run->SaveAs("dxy_phi_vs_" + append + ".pdf");
1573   c_dxy_phi_vs_run->SaveAs("dxy_phi_vs_" + append + ".png");
1574   c_dxy_phi_vs_run->SaveAs("dxy_phi_vs_" + append + ".eps");
1575
1576   c_dxy_eta_vs_run->SaveAs("dxy_eta_vs_" + append + ".pdf");
1577   c_dxy_eta_vs_run->SaveAs("dxy_eta_vs_" + append + ".png");
1578   c_dxy_eta_vs_run->SaveAs("dxy_eta_vs_" + append + ".eps");
1579
1580   c_dz_phi_vs_run->SaveAs("dz_phi_vs_" + append + ".pdf");
1581   c_dz_phi_vs_run->SaveAs("dz_phi_vs_" + append + ".png");
1582   c_dz_phi_vs_run->SaveAs("dz_phi_vs_" + append + ".eps");
1583
1584   c_dz_eta_vs_run->SaveAs("dz_eta_vs_" + append + ".pdf");
1585   c_dz_eta_vs_run->SaveAs("dz_eta_vs_" + append + ".png");
1586   c_dz_eta_vs_run->SaveAs("dz_eta_vs_" + append + ".eps");
1587
1588   TCanvas dummyC("dummyC", "dummyC", 2000, 800);
1589   dummyC.Print("combined.pdf[");
1590   c_dxy_phi_vs_run->Print("combined.pdf");
1591   c_dxy_eta_vs_run->Print("combined.pdf");
1592   c_dz_phi_vs_run->Print("combined.pdf");
1593   c_dz_eta_vs_run->Print("combined.pdf");
1594   dummyC.Print("combined.pdf]");
1595
1596   // mean
1597
1598   c_mean_dxy_phi_vs_run->SaveAs("mean_dxy_phi_vs_" + append + ".pdf");
1599   c_mean_dxy_phi_vs_run->SaveAs("mean_dxy_phi_vs_" + append + ".png");
1600
1601   c_mean_dxy_eta_vs_run->SaveAs("mean_dxy_eta_vs_" + append + ".pdf");
1602   c_mean_dxy_eta_vs_run->SaveAs("mean_dxy_eta_vs_" + append + ".png");
1603
1604   c_mean_dz_phi_vs_run->SaveAs("mean_dz_phi_vs_" + append + ".pdf");
1605   c_mean_dz_phi_vs_run->SaveAs("mean_dz_phi_vs_" + append + ".png");
1606
1607   c_mean_dz_eta_vs_run->SaveAs("mean_dz_eta_vs_" + append + ".pdf");
1608   c_mean_dz_eta_vs_run->SaveAs("mean_dz_eta_vs_" + append + ".png");
1609
1610   TCanvas dummyC2("dummyC2", "dummyC2", 2000, 800);
1611   dummyC2.Print("means.pdf[");
1612   c_mean_dxy_phi_vs_run->Print("means.pdf");
1613   c_mean_dxy_eta_vs_run->Print("means.pdf");
1614   c_mean_dz_phi_vs_run->Print("means.pdf");
1615   c_mean_dz_eta_vs_run->Print("means.pdf");
1616   dummyC2.Print("means.pdf]");
1617
1618   // RMS
1619
1620   c_RMS_dxy_phi_vs_run->SaveAs("RMS_dxy_phi_vs_" + append + ".pdf");
1621   c_RMS_dxy_phi_vs_run->SaveAs("RMS_dxy_phi_vs_" + append + ".png");
1622
1623   c_RMS_dxy_eta_vs_run->SaveAs("RMS_dxy_eta_vs_" + append + ".pdf");
1624   c_RMS_dxy_eta_vs_run->SaveAs("RMS_dxy_eta_vs_" + append + ".png");
1625
1626   c_RMS_dz_phi_vs_run->SaveAs("RMS_dz_phi_vs_" + append + ".pdf");
1627   c_RMS_dz_phi_vs_run->SaveAs("RMS_dz_phi_vs_" + append + ".png");
1628
1629   c_RMS_dz_eta_vs_run->SaveAs("RMS_dz_eta_vs_" + append + ".pdf");
1630   c_RMS_dz_eta_vs_run->SaveAs("RMS_dz_eta_vs_" + append + ".png");
1631
1632   TCanvas dummyC3("dummyC3", "dummyC3", 2000, 800);
1633   dummyC3.Print("RMSs.pdf[");
1634   c_RMS_dxy_phi_vs_run->Print("RMSs.pdf");
1635   c_RMS_dxy_eta_vs_run->Print("RMSs.pdf");
1636   c_RMS_dz_phi_vs_run->Print("RMSs.pdf");
1637   c_RMS_dz_eta_vs_run->Print("RMSs.pdf");
1638   dummyC3.Print("RMSs.pdf]");
1639
1640   // scatter
1641
1642   Scatter_dxy_vs_run->SaveAs("Scatter_dxy_vs_" + append + ".pdf");
1643   Scatter_dxy_vs_run->SaveAs("Scatter_dxy_vs_" + append + ".png");
1644
1645   Scatter_dz_vs_run->SaveAs("Scatter_dz_vs_" + append + ".pdf");
1646   Scatter_dz_vs_run->SaveAs("Scatter_dz_vs_" + append + ".png");
1647
1648   // chi2
1649
1650   c_chisquare_vs_run->SaveAs("chi2pol0fit_vs_" + append + ".pdf");
1651   c_chisquare_vs_run->SaveAs("chi2pol0fit_vs_" + append + ".png");
1652
1653   // KS score
1654
1655   c_KSScore_vs_run->SaveAs("KSScore_vs_" + append + ".pdf");
1656   c_KSScore_vs_run->SaveAs("KSScore_vs_" + append + ".png");
1657
1658   // do all the deletes
1659
1660   for (int iDir = 0; iDir < nDirs_; iDir++) {
1661     delete g_dxy_phi_vs_run[iDir];
1662     delete g_chi2_dxy_phi_vs_run[iDir];
1663     delete g_KS_dxy_phi_vs_run[iDir];
1664     delete g_dxy_phi_hi_vs_run[iDir];
1665     delete g_dxy_phi_lo_vs_run[iDir];
1666
1667     delete g_dxy_eta_vs_run[iDir];
1668     delete g_chi2_dxy_eta_vs_run[iDir];
1669     delete g_KS_dxy_eta_vs_run[iDir];
1670     delete g_dxy_eta_hi_vs_run[iDir];
1671     delete g_dxy_eta_lo_vs_run[iDir];
1672
1673     delete g_dz_phi_vs_run[iDir];
1674     delete g_chi2_dz_phi_vs_run[iDir];
1675     delete g_KS_dz_phi_vs_run[iDir];
1676     delete g_dz_phi_hi_vs_run[iDir];
1677     delete g_dz_phi_lo_vs_run[iDir];
1678
1679     delete g_dz_eta_vs_run[iDir];
1680     delete g_chi2_dz_eta_vs_run[iDir];
1681     delete g_KS_dz_eta_vs_run[iDir];
1682     delete g_dz_eta_hi_vs_run[iDir];
1683     delete g_dz_eta_lo_vs_run[iDir];
1684
1685     delete h_RMS_dxy_phi_vs_run[iDir];
1686     delete h_RMS_dxy_eta_vs_run[iDir];
1687     delete h_RMS_dz_phi_vs_run[iDir];
1688     delete h_RMS_dz_eta_vs_run[iDir];
1689   }
1690
1691   // mv the run-by-run plots into the folders
1692
1693   gSystem->mkdir("Biases");
1694   TString processline = ".! mv Bias*.p* ./Biases/";
1695   logInfo << "Executing: \n" << processline << "\n" << std::endl;
1696   gROOT->ProcessLine(processline.Data());
1697   gSystem->Sleep(100);
1698   processline.Clear();
1699
1700   gSystem->mkdir("ResolutionsVsPt");
1701   processline = ".! mv ResolutionsVsPt*.p* ./ResolutionsVsPt/";
1702   logInfo << "Executing: \n" << processline << "\n" << std::endl;
1703   gROOT->ProcessLine(processline.Data());
1704   gSystem->Sleep(100);
1705   processline.Clear();
1706
1707   gSystem->mkdir("Resolutions");
1708   processline = ".! mv Resolutions*.p* ./Resolutions/";
1709   logInfo << "Executing: \n" << processline << "\n" << std::endl;
1710   gROOT->ProcessLine(processline.Data());
1711   gSystem->Sleep(100);
1712   processline.Clear();
1713
1714   gSystem->mkdir("Pulls");
1715   processline = ".! mv Pulls*.p* ./Pulls/";
1716   logInfo << "Executing: \n" << processline << "\n" << std::endl;
1717   gROOT->ProcessLine(processline.Data());
1718   gSystem->Sleep(100);
1719   processline.Clear();
1720
1721   timer.Stop();
1722   timer.Print();
1723 }
1724
1725 /*! \fn outputGraphs
1726  *  \brief function to build the output graphs
1727  */
1728
1729 /*--------------------------------------------------------------------*/
1730 void outputGraphs(const pv::wrappedTrends &allInputs,
1731                   const std::vector<double> &ticks,
1732                   const std::vector<double> &ex_ticks,
1733                   TCanvas *&canv,
1734                   TCanvas *&mean_canv,
1735                   TCanvas *&rms_canv,
1736                   TGraph *&g_mean,
1737                   TGraph *&g_chi2,
1738                   TGraph *&g_KS,
1739                   TGraph *&g_low,
1740                   TGraph *&g_high,
1741                   TGraphAsymmErrors *&g_asym,
1742                   TH1F *h_RMS[],
1743                   const pv::bundle &mybundle,
1744                   const pv::view &theView,
1745                   const int index,
1746                   TObjArray *&array,
1747                   const TString &label,
1748                   TLegend *&legend)
1749 /*--------------------------------------------------------------------*/
1750 {
1751   g_mean = new TGraph(ticks.size(), &(ticks[0]), &((allInputs.getMean()[label])[0]));
1752   g_chi2 = new TGraph(ticks.size(), &(ticks[0]), &((allInputs.getChi2()[label])[0]));
1753   g_KS = new TGraph(ticks.size(), &(ticks[0]), &((allInputs.getKS()[label])[0]));
1754   g_high = new TGraph(ticks.size(), &(ticks[0]), &((allInputs.getHigh()[label])[0]));
1755   g_low = new TGraph(ticks.size(), &(ticks[0]), &((allInputs.getLow()[label])[0]));
1756
1757   g_asym = new TGraphAsymmErrors(ticks.size(),
1758                                  &(ticks[0]),
1759                                  &((allInputs.getMean()[label])[0]),
1760                                  &(ex_ticks[0]),
1761                                  &(ex_ticks[0]),
1762                                  &((allInputs.getLowErr()[label])[0]),
1763                                  &((allInputs.getHighErr()[label])[0]));
1764
1765   adjustmargins(canv);
1766   canv->cd();
1767   g_asym->SetFillStyle(3005);
1768   g_asym->SetFillColor(pv::colors[index]);
1769   g_mean->SetMarkerStyle(pv::markers[index]);
1770   g_mean->SetMarkerColor(pv::colors[index]);
1771   g_mean->SetLineColor(pv::colors[index]);
1772   g_mean->SetMarkerSize(1.);
1773   g_high->SetLineColor(pv::colors[index]);
1774   g_low->SetLineColor(pv::colors[index]);
1775   beautify(g_mean);
1776   beautify(g_asym);
1777
1778   if (theView == pv::dxyphi) {
1779     legend->AddEntry(g_mean, label, "PL");
1780   }
1781
1782   const char *coord;
1783   const char *kin;
1784   float ampl;
1785
1786   switch (theView) {
1787     case pv::dxyphi:
1788       coord = "xy";
1789       kin = "phi";
1790       ampl = 40;
1791       break;
1792     case pv::dzphi:
1793       coord = "z";
1794       kin = "phi";
1795       ampl = 40;
1796       break;
1797     case pv::dxyeta:
1798       coord = "xy";
1799       kin = "eta";
1800       ampl = 40;
1801       break;
1802     case pv::dzeta:
1803       coord = "z";
1804       kin = "eta";
1805       ampl = 60;
1806       break;
1807     default:
1808       coord = "unknown";
1809       kin = "unknown";
1810       break;
1811   }
1812
1813   g_mean->SetName(Form("g_bias_d%s_%s_%s", coord, kin, label.Data()));
1814   g_mean->SetTitle(Form("Bias of d_{%s}(#%s) vs %s", coord, kin, mybundle.getDataType()));
1815   g_mean->GetXaxis()->SetTitle(mybundle.getDataTypeLabel());
1816   g_mean->GetYaxis()->SetTitle(Form("#LT d_{%s}(#%s) #GT [#mum]", coord, kin));
1817   g_mean->GetYaxis()->SetRangeUser(-ampl, ampl);
1818
1819   logInfo << "===================================================================================================="
1820           << std::endl;
1821   logInfo << "Total Luminosity: " << mybundle.getTotalLumi() << std::endl;
1822   logInfo << "===================================================================================================="
1823           << std::endl;
1824
1825   g_asym->SetName(Form("gerr_bias_d%s_%s_%s", coord, kin, label.Data()));
1826   g_asym->SetTitle(Form("Bias of d_{%s}(#%s) vs %s", coord, kin, mybundle.getDataType()));
1827   g_asym->GetXaxis()->SetTitle(mybundle.getDataTypeLabel());
1828   g_asym->GetYaxis()->SetTitle(Form("#LT d_{%s}(#%s) #GT [#mum]", coord, kin));
1829   g_asym->GetYaxis()->SetRangeUser(-ampl, ampl);
1830
1831   //g_mean->GetXaxis()->UnZoom();
1832   //g_asym->GetXaxis()->UnZoom();
1833
1834   if (index == 0) {
1835     g_asym->GetXaxis()->SetRangeUser(0., mybundle.getTotalLumi());
1836     g_mean->Draw("AP");
1837     g_asym->Draw("P3same");
1838   } else {
1839     g_mean->Draw("Psame");
1840     g_asym->Draw("3same");
1841   }
1842   g_high->Draw("Lsame");
1843   g_low->Draw("Lsame");
1844
1845   if (mybundle.isTimeAxis()) {
1846     timify(g_asym);
1847     timify(g_mean);
1848     timify(g_high);
1849     timify(g_low);
1850   }
1851
1852   if (mybundle.isLumiAxis()) {
1853     g_asym->GetXaxis()->SetRangeUser(0., mybundle.getTotalLumi());
1854     g_mean->GetXaxis()->SetRangeUser(0., mybundle.getTotalLumi());
1855     g_high->GetXaxis()->SetRangeUser(0., mybundle.getTotalLumi());
1856     g_low->GetXaxis()->SetRangeUser(0., mybundle.getTotalLumi());
1857   }
1858
1859   if (index == (mybundle.getNObjects() - 1)) {
1860     legend->Draw("same");
1861     TH1F *theZero = DrawConstantGraph(g_mean, 1, 0.);
1862     theZero->Draw("E1same][");
1863   }
1864
1865   auto current_pad = static_cast<TPad *>(gPad);
1866   cmsPrel(current_pad);
1867
1868   superImposeIOVBoundaries(
1869       canv, mybundle.isLumiAxis(), mybundle.isTimeAxis(), mybundle.getLumiMapByRun(), mybundle.getTimes());
1870
1871   // mean only
1872   adjustmargins(mean_canv);
1873   mean_canv->cd();
1874   auto gprime = (TGraph *)g_mean->Clone();
1875   if (index == 0) {
1876     gprime->GetYaxis()->SetRangeUser(-10., 10.);
1877     if (mybundle.isLumiAxis())
1878       gprime->GetXaxis()->SetRangeUser(0., mybundle.getTotalLumi());
1879     gprime->Draw("AP");
1880   } else {
1881     gprime->Draw("Psame");
1882   }
1883
1884   if (index == (mybundle.getNObjects() - 1)) {
1885     legend->Draw("same");
1886   }
1887
1888   if (index == 0) {
1889     TH1F *theZero = DrawConstantGraph(gprime, 2, 0.);
1890     theZero->Draw("E1same][");
1891
1892     auto current_pad = static_cast<TPad *>(gPad);
1893     cmsPrel(current_pad);
1894
1895     superImposeIOVBoundaries(
1896         mean_canv, mybundle.isLumiAxis(), mybundle.isTimeAxis(), mybundle.getLumiMapByRun(), mybundle.getTimes());
1897   }
1898
1899   // scatter or RMS TH1
1900   h_RMS[index] = new TH1F(Form("h_RMS_dz_eta_%s", label.Data()),
1901                           Form("scatter of d_{%s}(#%s) vs %s;%s;%s of d_{%s}(#%s) [#mum]",
1902                                coord,
1903                                kin,
1904                                mybundle.getDataType(),
1905                                mybundle.getDataTypeLabel(),
1906                                (mybundle.isUsingRMS() ? "RMS" : "peak-to-peak deviation"),
1907                                coord,
1908                                kin),
1909                           ticks.size() - 1,
1910                           &(ticks[0]));
1911   h_RMS[index]->SetStats(kFALSE);
1912
1913   int bincounter = 0;
1914   for (const auto &tick : ticks) {
1915     bincounter++;
1916     h_RMS[index]->SetBinContent(
1917         bincounter, std::abs(allInputs.getHigh()[label][bincounter - 1] - allInputs.getLow()[label][bincounter - 1]));
1918     h_RMS[index]->SetBinError(bincounter, 0.01);
1919   }
1920
1921   h_RMS[index]->SetLineColor(pv::colors[index]);
1922   h_RMS[index]->SetLineWidth(2);
1923   h_RMS[index]->SetMarkerSize(1.);
1924   h_RMS[index]->SetMarkerStyle(pv::markers[index]);
1925   h_RMS[index]->SetMarkerColor(pv::colors[index]);
1926   adjustmargins(rms_canv);
1927   rms_canv->cd();
1928   beautify(h_RMS[index]);
1929
1930   if (mybundle.isTimeAxis()) {
1931     timify(h_RMS[index]);
1932   }
1933
1934   array->Add(h_RMS[index]);
1935
1936   // at the last file re-loop
1937   if (index == (mybundle.getNObjects() - 1)) {
1938     auto theMax = getMaximumFromArray(array);
1939     logInfo << "the max for d" << coord << "(" << kin << ") RMS is " << theMax << std::endl;
1940
1941     for (Int_t k = 0; k < mybundle.getNObjects(); k++) {
1942       //h_RMS[k]->GetYaxis()->SetRangeUser(-theMax*0.45,theMax*1.40);
1943       h_RMS[k]->GetYaxis()->SetRangeUser(0., theMax * 1.80);
1944       if (k == 0) {
1945         h_RMS[k]->Draw("L");
1946       } else {
1947         h_RMS[k]->Draw("Lsame");
1948       }
1949     }
1950     legend->Draw("same");
1951     TH1F *theConst = DrawConstant(h_RMS[index], 1, 0.);
1952     //theConst->Draw("same][");
1953
1954     current_pad = static_cast<TPad *>(gPad);
1955     cmsPrel(current_pad);
1956
1957     superImposeIOVBoundaries(
1958         rms_canv, mybundle.isLumiAxis(), mybundle.isTimeAxis(), mybundle.getLumiMapByRun(), mybundle.getTimes());
1959   }
1960 }
1961
1962 /*! \fn list_files
1963  *  \brief utility function to list of filles in a directory
1964  */
1965
1966 /*--------------------------------------------------------------------*/
1967 std::vector<int> list_files(const char *dirname, const char *ext)
1968 /*--------------------------------------------------------------------*/
1969 {
1970   std::vector<int> theRunNumbers;
1971
1972   TSystemDirectory dir(dirname, dirname);
1973   TList *files = dir.GetListOfFiles();
1974   if (files) {
1975     TSystemFile *file;
1976     TString fname;
1977     TIter next(files);
1978     while ((file = (TSystemFile *)next())) {
1979       fname = file->GetName();
1980       if (!file->IsDirectory() && fname.EndsWith(ext) && fname.BeginsWith("PVValidation")) {
1981         //logInfo << fname.Data() << std::endl;
1982         TObjArray *bits = fname.Tokenize("_");
1983         TString theRun = bits->At(2)->GetName();
1984         //logInfo << theRun << std::endl;
1985         TString formatRun = (theRun.ReplaceAll(".root", "")).ReplaceAll("_", "");
1986         //logInfo << dirname << " "<< formatRun.Atoi() << std::endl;
1987         theRunNumbers.push_back(formatRun.Atoi());
1988       }
1989     }
1990   }
1991   return theRunNumbers;
1992 }
1993
1994 /*! \fn arrangeOutCanvas
1995  *  \brief utility function to arrange the plots per run nicely in a TCanvas
1996  */
1997
1998 /*--------------------------------------------------------------------*/
1999 void arrangeOutCanvas(TCanvas *canv,
2000                       TH1F *m_11Trend[100],
2001                       TH1F *m_12Trend[100],
2002                       TH1F *m_21Trend[100],
2003                       TH1F *m_22Trend[100],
2004                       Int_t nDirs,
2005                       TString LegLabels[10],
2006                       unsigned int theRun) {
2007   /*--------------------------------------------------------------------*/
2008
2009   TLegend *lego = new TLegend(0.19, 0.80, 0.79, 0.93);
2010   //lego-> SetNColumns(2);
2011   lego->SetFillColor(10);
2012   lego->SetTextSize(0.042);
2013   lego->SetTextFont(42);
2014   lego->SetFillColor(10);
2015   lego->SetLineColor(10);
2016   lego->SetShadowColor(10);
2017
2018   TPaveText *ptDate = new TPaveText(0.19, 0.95, 0.45, 0.99, "blNDC");
2019   ptDate->SetFillColor(kYellow);
2020   //ptDate->SetFillColor(10);
2021   ptDate->SetBorderSize(1);
2022   ptDate->SetLineColor(kBlue);
2023   ptDate->SetLineWidth(1);
2024   ptDate->SetTextFont(42);
2025   TText *textDate = ptDate->AddText(Form("Run: %i", theRun));
2026   textDate->SetTextSize(0.04);
2027   textDate->SetTextColor(kBlue);
2028   textDate->SetTextAlign(22);
2029
2030   canv->SetFillColor(10);
2031   canv->Divide(2, 2);
2032
2033   TH1F *dBiasTrend[4][nDirs];
2034
2035   for (Int_t i = 0; i < nDirs; i++) {
2036     dBiasTrend[0][i] = m_11Trend[i];
2037     dBiasTrend[1][i] = m_12Trend[i];
2038     dBiasTrend[2][i] = m_21Trend[i];
2039     dBiasTrend[3][i] = m_22Trend[i];
2040   }
2041
2042   Double_t absmin[4] = {999., 999., 999., 999.};
2043   Double_t absmax[4] = {-999., -999. - 999., -999.};
2044
2045   for (Int_t k = 0; k < 4; k++) {
2046     canv->cd(k + 1)->SetBottomMargin(0.14);
2047     canv->cd(k + 1)->SetLeftMargin(0.18);
2048     canv->cd(k + 1)->SetRightMargin(0.01);
2049     canv->cd(k + 1)->SetTopMargin(0.06);
2050
2051     canv->cd(k + 1);
2052
2053     for (Int_t i = 0; i < nDirs; i++) {
2054       if (dBiasTrend[k][i]->GetMaximum() > absmax[k])
2055         absmax[k] = dBiasTrend[k][i]->GetMaximum();
2056       if (dBiasTrend[k][i]->GetMinimum() < absmin[k])
2057         absmin[k] = dBiasTrend[k][i]->GetMinimum();
2058     }
2059
2060     Double_t safeDelta = (absmax[k] - absmin[k]) / 8.;
2061     Double_t theExtreme = std::max(absmax[k], TMath::Abs(absmin[k]));
2062
2063     for (Int_t i = 0; i < nDirs; i++) {
2064       if (i == 0) {
2065         TString theTitle = dBiasTrend[k][i]->GetName();
2066
2067         if (theTitle.Contains("norm")) {
2068           //dBiasTrend[k][i]->GetYaxis()->SetRangeUser(std::min(-0.48,absmin[k]-safeDelta/2.),std::max(0.48,absmax[k]+safeDelta/2.));
2069           dBiasTrend[k][i]->GetYaxis()->SetRangeUser(0., 1.8);
2070         } else {
2071           if (!theTitle.Contains("width")) {
2072             if (theTitle.Contains("ladder") || theTitle.Contains("modZ")) {
2073               dBiasTrend[k][i]->GetYaxis()->SetRangeUser(std::min(-40., -theExtreme - (safeDelta / 2.)),
2074                                                          std::max(40., theExtreme + (safeDelta / 2.)));
2075             } else {
2076               dBiasTrend[k][i]->GetYaxis()->SetRangeUser(-40., 40.);
2077             }
2078           } else {
2079             // dBiasTrend[k][i]->GetYaxis()->SetRangeUser(0.,theExtreme+(safeDelta/2.));
2080             // if(theTitle.Contains("eta")) {
2081             //   dBiasTrend[k][i]->GetYaxis()->SetRangeUser(0.,500.);
2082             // } else {
2083             //   dBiasTrend[k][i]->GetYaxis()->SetRangeUser(0.,200.);
2084             // }
2085             auto my_view = checkTheView(theTitle);
2086
2087             //logInfo<<" ----------------------------------> " << theTitle << " view: " << my_view <<  std::endl;
2088
2089             switch (my_view) {
2090               case pv::dxyphi:
2091                 dBiasTrend[k][i]->GetYaxis()->SetRangeUser(0., 200.);
2092                 break;
2093               case pv::dzphi:
2094                 dBiasTrend[k][i]->GetYaxis()->SetRangeUser(0., 400.);
2095                 break;
2096               case pv::dxyeta:
2097                 dBiasTrend[k][i]->GetYaxis()->SetRangeUser(0., 300.);
2098                 break;
2099               case pv::dzeta:
2100                 dBiasTrend[k][i]->GetYaxis()->SetRangeUser(0., 1.e3);
2101                 break;
2102               case pv::pT:
2103                 dBiasTrend[k][i]->GetYaxis()->SetRangeUser(0., 200.);
2104                 break;
2105               case pv::generic:
2106                 dBiasTrend[k][i]->GetYaxis()->SetRangeUser(0., 350.);
2107                 break;
2108               default:
2109                 dBiasTrend[k][i]->GetYaxis()->SetRangeUser(0., 300.);
2110             }
2111           }
2112         }
2113
2114         dBiasTrend[k][i]->Draw("Le1");
2115         makeNewXAxis(dBiasTrend[k][i]);
2116
2117         Double_t theC = -1.;
2118
2119         if (theTitle.Contains("width")) {
2120           if (theTitle.Contains("norm")) {
2121             theC = 1.;
2122           } else {
2123             theC = -1.;
2124           }
2125         } else {
2126           theC = 0.;
2127         }
2128
2129         TH1F *theConst = DrawConstant(dBiasTrend[k][i], 1, theC);
2130         theConst->Draw("PLsame][");
2131
2132       } else {
2133         dBiasTrend[k][i]->Draw("Le1sames");
2134         makeNewXAxis(dBiasTrend[k][i]);
2135       }
2136       TPad *current_pad = static_cast<TPad *>(canv->GetPad(k + 1));
2137       cmsPrel(current_pad, 2);
2138       ptDate->Draw("same");
2139
2140       if (k == 0) {
2141         lego->AddEntry(dBiasTrend[k][i], LegLabels[i]);
2142       }
2143     }
2144
2145     lego->Draw();
2146   }
2147 }
2148
2149 /*! \fn MakeNiceTrendPlotStype
2150  *  \brief utility function to embellish trend plot style
2151  */
2152
2153 /*--------------------------------------------------------------------*/
2154 void MakeNiceTrendPlotStyle(TH1 *hist, Int_t color)
2155 /*--------------------------------------------------------------------*/
2156 {
2157   hist->SetStats(kFALSE);
2158   hist->SetLineWidth(2);
2159   hist->GetXaxis()->CenterTitle(true);
2160   hist->GetYaxis()->CenterTitle(true);
2161   hist->GetXaxis()->SetTitleFont(42);
2162   hist->GetYaxis()->SetTitleFont(42);
2163   hist->GetXaxis()->SetTitleSize(0.065);
2164   hist->GetYaxis()->SetTitleSize(0.065);
2165   hist->GetXaxis()->SetTitleOffset(1.0);
2166   hist->GetYaxis()->SetTitleOffset(1.2);
2167   hist->GetXaxis()->SetLabelFont(42);
2168   hist->GetYaxis()->SetLabelFont(42);
2169   hist->GetYaxis()->SetLabelSize(.05);
2170   hist->GetXaxis()->SetLabelSize(.07);
2171   //hist->GetXaxis()->SetNdivisions(505);
2172   if (color != 8) {
2173     hist->SetMarkerSize(1.5);
2174   } else {
2175     hist->SetLineWidth(3);
2176     hist->SetMarkerSize(0.0);
2177   }
2178   hist->SetMarkerStyle(pv::markers[color]);
2179   hist->SetLineColor(pv::colors[color]);
2180   hist->SetMarkerColor(pv::colors[color]);
2181 }
2182
2183 /*! \fn maxkeNewAxis
2184  *  \brief utility function to re-make x-axis with correct binning
2185  */
2186
2187 /*--------------------------------------------------------------------*/
2188 void makeNewXAxis(TH1 *h)
2189 /*--------------------------------------------------------------------*/
2190 {
2191   TString myTitle = h->GetName();
2192   float axmin = -999;
2193   float axmax = 999.;
2194   int ndiv = 510;
2195   if (myTitle.Contains("eta")) {
2196     axmin = -2.7;
2197     axmax = 2.7;
2198     ndiv = 505;
2199   } else if (myTitle.Contains("phi")) {
2200     axmin = -TMath::Pi();
2201     axmax = TMath::Pi();
2202     ndiv = 510;
2203   } else if (myTitle.Contains("pT")) {
2204     axmin = 0;
2205     axmax = 19.99;
2206     ndiv = 510;
2207   } else if (myTitle.Contains("ladder")) {
2208     axmin = 0;
2209     axmax = 12;
2210     ndiv = 510;
2211   } else if (myTitle.Contains("modZ")) {
2212     axmin = 0;
2213     axmax = 8;
2214     ndiv = 510;
2215   } else {
2216     logError << "unrecognized variable" << std::endl;
2217   }
2218
2219   // Remove the current axis
2220   h->GetXaxis()->SetLabelOffset(999);
2221   h->GetXaxis()->SetTickLength(0);
2222
2223   if (myTitle.Contains("phi")) {
2224     h->GetXaxis()->SetTitle("#varphi (sector) [rad]");
2225   }
2226
2227   // Redraw the new axis
2228   gPad->Update();
2229
2230   TGaxis *newaxis =
2231       new TGaxis(gPad->GetUxmin(), gPad->GetUymin(), gPad->GetUxmax(), gPad->GetUymin(), axmin, axmax, ndiv, "SDH");
2232
2233   TGaxis *newaxisup =
2234       new TGaxis(gPad->GetUxmin(), gPad->GetUymax(), gPad->GetUxmax(), gPad->GetUymax(), axmin, axmax, ndiv, "-SDH");
2235
2236   newaxis->SetLabelOffset(0.02);
2237   newaxis->SetLabelFont(42);
2238   newaxis->SetLabelSize(0.05);
2239
2240   newaxisup->SetLabelOffset(-0.02);
2241   newaxisup->SetLabelFont(42);
2242   newaxisup->SetLabelSize(0);
2243
2244   newaxis->Draw();
2245   newaxisup->Draw();
2246 }
2247
2248 /*! \fn setStyle
2249  *  \brief main function to set the plotting style
2250  */
2251
2252 /*--------------------------------------------------------------------*/
2253 void setStyle() {
2254   /*--------------------------------------------------------------------*/
2255
2256   TGaxis::SetMaxDigits(6);
2257
2258   TH1::StatOverflows(kTRUE);
2259   gStyle->SetOptTitle(0);
2260   gStyle->SetOptStat("e");
2261   //gStyle->SetPadTopMargin(0.05);
2262   //gStyle->SetPadBottomMargin(0.15);
2263   //gStyle->SetPadLeftMargin(0.17);
2264   //gStyle->SetPadRightMargin(0.02);
2265   gStyle->SetPadBorderMode(0);
2266   gStyle->SetTitleFillColor(10);
2267   gStyle->SetTitleFont(42);
2268   gStyle->SetTitleColor(1);
2269   gStyle->SetTitleTextColor(1);
2270   gStyle->SetTitleFontSize(0.06);
2271   gStyle->SetTitleBorderSize(0);
2272   gStyle->SetStatColor(kWhite);
2273   gStyle->SetStatFont(42);
2274   gStyle->SetStatFontSize(0.05);  ///---> gStyle->SetStatFontSize(0.025);
2275   gStyle->SetStatTextColor(1);
2276   gStyle->SetStatFormat("6.4g");
2277   gStyle->SetStatBorderSize(1);
2278   gStyle->SetPadTickX(1);  // To get tick marks on the opposite side of the frame
2279   gStyle->SetPadTickY(1);
2280   gStyle->SetPadBorderMode(0);
2281   gStyle->SetOptFit(1);
2282   gStyle->SetNdivisions(510);
2283
2284   //gStyle->SetPalette(kInvertedDarkBodyRadiator);
2285
2286   const Int_t NRGBs = 5;
2287   const Int_t NCont = 255;
2288
2289   /*
2290   Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
2291   Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
2292   Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
2293   Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
2294   */
2295
2296   Double_t stops[NRGBs] = {0.00, 0.01, 0.05, 0.09, 0.1};
2297   Double_t red[NRGBs] = {1.00, 0.84, 0.61, 0.34, 0.00};
2298   Double_t green[NRGBs] = {1.00, 0.84, 0.61, 0.34, 0.00};
2299   Double_t blue[NRGBs] = {1.00, 0.84, 0.61, 0.34, 0.00};
2300
2301   TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
2302   gStyle->SetNumberContours(NCont);
2303 }
2304
2305 // /*--------------------------------------------------------------------*/
2306 // void cmsPrel(TPad* pad) {
2307 // /*--------------------------------------------------------------------*/
2308
2309 //   float H = pad->GetWh();
2310 //   float W = pad->GetWw();
2311 //   float l = pad->GetLeftMargin();
2312 //   float t = pad->GetTopMargin();
2313 //   float r = pad->GetRightMargin();
2314 //   float b = pad->GetBottomMargin();
2315 //   float relPosX = 0.009;
2316 //   float relPosY = 0.045;
2317 //   float lumiTextOffset = 0.8;
2318
2319 //   TLatex *latex = new TLatex();
2320 //   latex->SetNDC();
2321 //   latex->SetTextSize(0.045);
2322
2323 //   float posX_    = 1-r - relPosX*(1-l-r);
2324 //   float posXCMS_ = posX_- 0.125;
2325 //   float posY_ =  1-t + 0.05; /// - relPosY*(1-t-b);
2326 //   float factor = 1./0.86;
2327
2328 //   latex->SetTextAlign(33);
2329 //   latex->SetTextFont(61);
2330 //   latex->SetTextSize(0.045*factor);
2331 //   latex->DrawLatex(posXCMS_,posY_,"CMS");
2332 //   latex->SetTextSize(0.045);
2333 //   latex->SetTextFont(42); //22
2334 //   latex->DrawLatex(posX_,posY_,"Internal (13 TeV)");
2335 //   //latex->DrawLatex(posX_,posY_,"CMS Preliminary (13 TeV)");
2336 //   //latex->DrawLatex(posX_,posY_,"CMS 2017 Work in progress (13 TeV)");
2337
2338 // }
2339
2340 /*! \fn cmsPrel
2341  *  \brief utility function to put the CMS paraphernalia on canvas
2342  */
2343
2344 /*--------------------------------------------------------------------*/
2345 void cmsPrel(TPad *pad, size_t ipads) {
2346   /*--------------------------------------------------------------------*/
2347
2348   float H = pad->GetWh();
2349   float W = pad->GetWw();
2350   float l = pad->GetLeftMargin();
2351   float t = pad->GetTopMargin();
2352   float r = pad->GetRightMargin();
2353   float b = pad->GetBottomMargin();
2354   float relPosX = 0.009;
2355   float relPosY = 0.045;
2356   float lumiTextOffset = 0.8;
2357
2358   TLatex *latex = new TLatex();
2359   latex->SetNDC();
2360   latex->SetTextSize(0.045);
2361
2362   float posX_ = 1 - (r / ipads);
2363   float posY_ = 1 - t + 0.05;  /// - relPosY*(1-t-b);
2364   float factor = 1. / 0.82;
2365
2366   latex->SetTextAlign(33);
2367   latex->SetTextSize(0.045);
2368   latex->SetTextFont(42);  //22
2369   latex->DrawLatex(posX_, posY_, "Internal (13 TeV)");
2370
2371   UInt_t w;
2372   UInt_t h;
2373   latex->GetTextExtent(w, h, "Internal (13 TeV)");
2374   float size = w / (W / ipads);
2375   //logInfo<<w<<" "<<" "<<W<<" "<<size<<std::endl;
2376   float posXCMS_ = posX_ - size * (1 + 0.025 * ipads);
2377
2378   latex->SetTextAlign(33);
2379   latex->SetTextFont(61);
2380   latex->SetTextSize(0.045 * factor);
2381   latex->DrawLatex(posXCMS_, posY_ + 0.004, "CMS");
2382
2383   //latex->DrawLatex(posX_,posY_,"CMS Preliminary (13 TeV)");
2384   //latex->DrawLatex(posX_,posY_,"CMS 2017 Work in progress (13 TeV)");
2385 }
2386
2387 /*! \fn DrawConstant
2388  *  \brief utility function to draw a constant histogram
2389  */
2390
2391 /*--------------------------------------------------------------------*/
2392 TH1F *DrawConstant(TH1F *hist, Int_t iter, Double_t theConst)
2393 /*--------------------------------------------------------------------*/
2394 {
2395   Int_t nbins = hist->GetNbinsX();
2396   Double_t lowedge = hist->GetBinLowEdge(1);
2397   Double_t highedge = hist->GetBinLowEdge(nbins + 1);
2398
2399   TH1F *hzero = new TH1F(Form("hconst_%s_%i", hist->GetName(), iter),
2400                          Form("hconst_%s_%i", hist->GetName(), iter),
2401                          nbins,
2402                          lowedge,
2403                          highedge);
2404   for (Int_t i = 0; i <= hzero->GetNbinsX(); i++) {
2405     hzero->SetBinContent(i, theConst);
2406     hzero->SetBinError(i, 0.);
2407   }
2408   hzero->SetLineWidth(2);
2409   hzero->SetLineStyle(9);
2410   hzero->SetLineColor(kMagenta);
2411
2412   return hzero;
2413 }
2414
2415 /*! \fn DrawConstant
2416  *  \brief utility function to draw a constant histogram with erros !=0
2417  */
2418
2419 /*--------------------------------------------------------------------*/
2420 TH1F *DrawConstantWithErr(TH1F *hist, Int_t iter, Double_t theConst)
2421 /*--------------------------------------------------------------------*/
2422 {
2423   Int_t nbins = hist->GetNbinsX();
2424   Double_t lowedge = hist->GetBinLowEdge(1);
2425   Double_t highedge = hist->GetBinLowEdge(nbins + 1);
2426
2427   TH1F *hzero = new TH1F(Form("hconst_%s_%i", hist->GetName(), iter),
2428                          Form("hconst_%s_%i", hist->GetName(), iter),
2429                          nbins,
2430                          lowedge,
2431                          highedge);
2432   for (Int_t i = 0; i <= hzero->GetNbinsX(); i++) {
2433     hzero->SetBinContent(i, theConst);
2434     hzero->SetBinError(i, hist->GetBinError(i));
2435   }
2436   hzero->SetLineWidth(2);
2437   hzero->SetLineStyle(9);
2438   hzero->SetLineColor(kMagenta);
2439
2440   return hzero;
2441 }
2442
2443 /*! \fn DrawConstantGraph
2444  *  \brief utility function to draw a constant TGraph
2445  */
2446
2447 /*--------------------------------------------------------------------*/
2448 TH1F *DrawConstantGraph(TGraph *graph, Int_t iter, Double_t theConst)
2449 /*--------------------------------------------------------------------*/
2450 {
2451   Double_t xmin = graph->GetXaxis()->GetXmin();  //TMath::MinElement(graph->GetN(),graph->GetX());
2452   Double_t xmax = graph->GetXaxis()->GetXmax();  //TMath::MaxElement(graph->GetN(),graph->GetX());
2453
2454   //logInfo<<xmin<<" : "<<xmax<<std::endl;
2455
2456   TH1F *hzero = new TH1F(Form("hconst_%s_%i", graph->GetName(), iter),
2457                          Form("hconst_%s_%i", graph->GetName(), iter),
2458                          graph->GetN(),
2459                          xmin,
2460                          xmax);
2461   for (Int_t i = 0; i <= hzero->GetNbinsX(); i++) {
2462     hzero->SetBinContent(i, theConst);
2463     hzero->SetBinError(i, 0.);
2464   }
2465
2466   hzero->SetLineWidth(2);
2467   hzero->SetLineStyle(9);
2468   hzero->SetLineColor(kMagenta);
2469
2470   return hzero;
2471 }
2472
2473 /*! \fn getUnrolledHisto
2474  *  \brief utility function to tranform a TH1 into a vector of floats
2475  */
2476
2477 /*--------------------------------------------------------------------*/
2478 unrolledHisto getUnrolledHisto(TH1F *hist)
2479 /*--------------------------------------------------------------------*/
2480 {
2481   /*
2482     Double_t y_min = hist->GetBinLowEdge(1);
2483     Double_t y_max = hist->GetBinLowEdge(hist->GetNbinsX()+1);
2484   */
2485
2486   Double_t y_min = -0.1;
2487   Double_t y_max = 0.1;
2488
2489   std::vector<Double_t> contents;
2490   for (int j = 0; j < hist->GetNbinsX(); j++) {
2491     if (std::abs(hist->GetXaxis()->GetBinCenter(j)) <= 0.1)
2492       contents.push_back(hist->GetBinContent(j + 1));
2493   }
2494
2495   auto ret = unrolledHisto(y_min, y_max, contents.size(), contents);
2496   return ret;
2497 }
2498
2499 /*! \fn getBiases
2500  *  \brief utility function to extract characterization of the PV bias plot
2501  */
2502
2503 /*--------------------------------------------------------------------*/
2504 pv::biases getBiases(TH1F *hist)
2505 /*--------------------------------------------------------------------*/
2506 {
2507   int nbins = hist->GetNbinsX();
2508
2509   //extract median from histogram
2510   double *y = new double[nbins];
2511   double *err = new double[nbins];
2512
2513   // remember for weight means <x> = sum_i (x_i* w_i) / sum_i w_i ; where w_i = 1/sigma^2_i
2514
2515   for (int j = 0; j < nbins; j++) {
2516     y[j] = hist->GetBinContent(j + 1);
2517     if (hist->GetBinError(j + 1) != 0.) {
2518       err[j] = 1. / (hist->GetBinError(j + 1) * hist->GetBinError(j + 1));
2519     } else {
2520       err[j] = 0.;
2521     }
2522   }
2523
2524   Double_t w_mean = TMath::Mean(nbins, y, err);
2525   Double_t w_rms = TMath::RMS(nbins, y, err);
2526
2527   Double_t mean = TMath::Mean(nbins, y);
2528   Double_t rms = TMath::RMS(nbins, y);
2529
2530   Double_t max = hist->GetMaximum();
2531   Double_t min = hist->GetMinimum();
2532
2533   // in case one would like to use a pol0 fit
2534   hist->Fit("pol0", "Q0+");
2535   TF1 *f = (TF1 *)hist->FindObject("pol0");
2536   //f->SetLineColor(hist->GetLineColor());
2537   //f->SetLineStyle(hist->GetLineStyle());
2538   Double_t chi2 = f->GetChisquare();
2539   Int_t ndf = f->GetNDF();
2540
2541   TH1F *theZero = DrawConstantWithErr(hist, 1, 1.);
2542   TH1F *displaced = (TH1F *)hist->Clone("displaced");
2543   displaced->Add(theZero);
2544   Double_t ksScore = std::max(-20., TMath::Log10(displaced->KolmogorovTest(theZero)));
2545   Double_t chi2Score = displaced->Chi2Test(theZero);
2546
2547   /*
2548     std::pair<std::pair<Double_t,Double_t>, Double_t> result;
2549     std::pair<Double_t,Double_t> resultBounds;
2550     resultBounds = useRMS_ ? std::make_pair(mean-rms,mean+rms) :  std::make_pair(min,max)  ;
2551     result = make_pair(resultBounds,mean);
2552   */
2553
2554   pv::biases result(mean, rms, w_mean, w_rms, min, max, chi2, ndf, ksScore);
2555
2556   delete theZero;
2557   delete displaced;
2558   return result;
2559 }
2560
2561 /*! \fn beautify
2562  *  \brief utility function to beautify a TGraph
2563  */
2564
2565 /*--------------------------------------------------------------------*/
2566 void beautify(TGraph *g) {
2567   /*--------------------------------------------------------------------*/
2568   g->GetXaxis()->SetLabelFont(42);
2569   g->GetYaxis()->SetLabelFont(42);
2570   g->GetYaxis()->SetLabelSize(.055);
2571   g->GetXaxis()->SetLabelSize(.055);
2572   g->GetYaxis()->SetTitleSize(.055);
2573   g->GetXaxis()->SetTitleSize(.055);
2574   g->GetXaxis()->SetTitleOffset(1.1);
2575   g->GetYaxis()->SetTitleOffset(0.6);
2576   g->GetXaxis()->SetTitleFont(42);
2577   g->GetYaxis()->SetTitleFont(42);
2578   g->GetXaxis()->CenterTitle(true);
2579   g->GetYaxis()->CenterTitle(true);
2580   g->GetXaxis()->SetNdivisions(505);
2581 }
2582
2583 /*! \fn beautify
2584  *  \brief utility function to beautify a TH1
2585  */
2586
2587 /*--------------------------------------------------------------------*/
2588 void beautify(TH1 *h) {
2589   /*--------------------------------------------------------------------*/
2590   h->SetMinimum(0.);
2591   h->GetXaxis()->SetLabelFont(42);
2592   h->GetYaxis()->SetLabelFont(42);
2593   h->GetYaxis()->SetLabelSize(.055);
2594   h->GetXaxis()->SetLabelSize(.055);
2595   h->GetYaxis()->SetTitleSize(.055);
2596   h->GetXaxis()->SetTitleSize(.055);
2597   h->GetXaxis()->SetTitleOffset(1.1);
2598   h->GetYaxis()->SetTitleOffset(0.6);
2599   h->GetXaxis()->SetTitleFont(42);
2600   h->GetYaxis()->SetTitleFont(42);
2601   h->GetXaxis()->CenterTitle(true);
2602   h->GetYaxis()->CenterTitle(true);
2603   h->GetXaxis()->SetNdivisions(505);
2604 }
2605
2606 /*! \fn adjustmargins
2607  *  \brief utility function to adjust margins of a TCanvas
2608  */
2609
2610 /*--------------------------------------------------------------------*/
2611 void adjustmargins(TCanvas *canv) {
2612   /*--------------------------------------------------------------------*/
2613   canv->cd()->SetBottomMargin(0.14);
2614   canv->cd()->SetLeftMargin(0.07);
2615   canv->cd()->SetRightMargin(0.03);
2616   canv->cd()->SetTopMargin(0.06);
2617 }
2618
2619 /*! \fn adjustmargins
2620  *  \brief utility function to adjust margins of a TVirtualPad
2621  */
2622
2623 /*--------------------------------------------------------------------*/
2624 void adjustmargins(TVirtualPad *canv) {
2625   /*--------------------------------------------------------------------*/
2626   canv->SetBottomMargin(0.12);
2627   canv->SetLeftMargin(0.07);
2628   canv->SetRightMargin(0.01);
2629   canv->SetTopMargin(0.02);
2630 }
2631
2632 /*! \fn checkTH1AndReturn
2633  *  \brief utility function to check if a TH1 exists in a TFile before returning it
2634  */
2635
2636 /*--------------------------------------------------------------------*/
2637 TH1F *checkTH1AndReturn(TFile *f, TString address) {
2638   /*--------------------------------------------------------------------*/
2639   TH1F *h = NULL;
2640   if (f->GetListOfKeys()->Contains(address)) {
2641     h = (TH1F *)f->Get(address);
2642   }
2643   return h;
2644 }
2645
2646 /*! \fn checkThView
2647  *  \brief utility function to return the pv::view corresponding to a given PV bias plot
2648  */
2649
2650 /*--------------------------------------------------------------------*/
2651 pv::view checkTheView(const TString &toCheck) {
2652   /*--------------------------------------------------------------------*/
2653   if (toCheck.Contains("dxy")) {
2654     if (toCheck.Contains("phi") || toCheck.Contains("ladder")) {
2655       return pv::dxyphi;
2656     } else if (toCheck.Contains("eta") || toCheck.Contains("modZ")) {
2657       return pv::dxyeta;
2658     } else {
2659       return pv::pT;
2660     }
2661   } else if (toCheck.Contains("dz")) {
2662     if (toCheck.Contains("phi") || toCheck.Contains("ladder")) {
2663       return pv::dzphi;
2664     } else if (toCheck.Contains("eta") || toCheck.Contains("modZ")) {
2665       return pv::dzeta;
2666     } else {
2667       return pv::pT;
2668     }
2669   } else {
2670     return pv::generic;
2671   }
2672 }
2673
2674 /*! \fn timify
2675  *  \brief utility function to make a histogram x-axis of the time type
2676  */
2677
2678 /*--------------------------------------------------------------------*/
2679 template <typename T>
2680 void timify(T *mgr)
2681 /*--------------------------------------------------------------------*/
2682 {
2683   mgr->GetXaxis()->SetTimeDisplay(1);
2684   mgr->GetXaxis()->SetNdivisions(510);
2685   mgr->GetXaxis()->SetTimeFormat("%Y-%m-%d");
2686   mgr->GetXaxis()->SetTimeOffset(0, "gmt");
2687   mgr->GetXaxis()->SetLabelSize(.035);
2688 }
2689
2690 struct increase {
2691   template <class T>
2692   bool operator()(T const &a, T const &b) const {
2693     return a > b;
2694   }
2695 };
2696
2697 /*! \fn getMaximumFromArray
2698  *  \brief utility function to extract the maximum out of an array of histograms
2699  */
2700
2701 /*--------------------------------------------------------------------*/
2702 Double_t getMaximumFromArray(TObjArray *array)
2703 /*--------------------------------------------------------------------*/
2704 {
2705   Double_t theMaximum = -999.;  //(static_cast<TH1*>(array->At(0)))->GetMaximum();
2706
2707   for (Int_t i = 0; i < array->GetSize(); i++) {
2708     double theMaxForThisHist;
2709     auto hist = static_cast<TH1 *>(array->At(i));
2710     std::vector<double> maxima;
2711     for (int j = 0; j < hist->GetNbinsX(); j++)
2712       maxima.push_back(hist->GetBinContent(j));
2713     std::sort(std::begin(maxima), std::end(maxima));  //,increase());
2714     double rms_maxima = TMath::RMS(hist->GetNbinsX(), &(maxima[0]));
2715     double mean_maxima = TMath::Mean(hist->GetNbinsX(), &(maxima[0]));
2716
2717     const Int_t nq = 100;
2718     Double_t xq[nq];  // position where to compute the quantiles in [0,1]
2719     Double_t yq[nq];  // array to contain the quantiles
2720     for (Int_t i = 0; i < nq; i++)
2721       xq[i] = 0.9 + (Float_t(i + 1) / nq) * 0.10;
2722     TMath::Quantiles(maxima.size(), nq, &(maxima[0]), yq, xq);
2723
2724     //for(int q=0;q<nq;q++){
2725     //  logInfo<<q<<" "<<xq[q]<<" "<<yq[q]<<std::endl;
2726     //}
2727
2728     //for (const auto &element : maxima){
2729     //  if(element<1.5*mean_maxima){
2730     //  theMaxForThisHist=element;
2731     //  break;
2732     //  }
2733     //}
2734
2735     theMaxForThisHist = yq[80];
2736
2737     logInfo << "rms_maxima[" << i << "]: " << rms_maxima << " mean maxima[" << mean_maxima
2738             << "] purged maximum:" << theMaxForThisHist << std::endl;
2739
2740     if (theMaxForThisHist > theMaximum)
2741       theMaximum = theMaxForThisHist;
2742
2743     /*
2744     if( (static_cast<TH1*>(array->At(i)))->GetMaximum() > theMaximum){
2745       theMaximum = (static_cast<TH1*>(array->At(i)))->GetMaximum();
2746       //logInfo<<"i= "<<i<<" theMaximum="<<theMaximum<<endl;
2747     }
2748     */
2749   }
2750
2751   return theMaximum;
2752 }
2753
2754 /*! \fn superImposeIOVBoundaries
2755  *  \brief utility function to superimpose the IOV boundaries on the trend plot
2756  */
2757
2758 /*--------------------------------------------------------------------*/
2759 void superImposeIOVBoundaries(TCanvas *c,
2760                               bool lumi_axis_format,
2761                               bool time_axis_format,
2762                               const std::map<int, double> &lumiMapByRun,
2763                               const std::map<int, TDatime> &timeMap,
2764                               bool drawText)
2765 /*--------------------------------------------------------------------*/
2766 {
2767   // get the vector of runs in the lumiMap
2768   std::vector<int> vruns;
2769   for (auto const &imap : lumiMapByRun) {
2770     vruns.push_back(imap.first);
2771     //logInfo<<" run:" << imap.first << " lumi: "<< imap.second << std::endl;
2772   }
2773
2774   //std::vector<vint> truns;
2775   for (auto const &imap : timeMap) {
2776     logInfo << " run:" << imap.first << " time: " << imap.second.Convert() << std::endl;
2777   }
2778
2779   /* Run-2 Ultra-Legacy ReReco IOVs (from tag SiPixelTemplateDBObject_38T_v16_offline)
2780      1      271866   2016-04-28   2cc9ecb98e8ba900b26701d6adf052ba59c1f5e1  SiPixelTemplateDBObject 2016
2781      2      276315   2016-07-04   28a6077528012fe0abb03df9ef52e5976137c324  SiPixelTemplateDBObject
2782      3      278271   2016-08-05   43744865be66299a3a37599961a9faf5a0d2bd78  SiPixelTemplateDBObject
2783      4      280928   2016-09-16   7533fd5b60c10a9c55cddc6d2677d3fae6e8bb14  SiPixelTemplateDBObject
2784      5      290543   2017-03-30   05803622f45cf4ee2c07af2bb97875bd1010bfea  SiPixelTemplateDBObject 2017
2785      6      297281   2017-06-22   0ed971165dda7ae1db50c6636dab049087cad9a1  SiPixelTemplateDBObject
2786      7      298653   2017-07-09   e3af935c8f4eded04455a88959501d7408895501  SiPixelTemplateDBObject
2787      8      299443   2017-07-19   e2203ad6babf1885d754ae1392e4841a9b20c02e  SiPixelTemplateDBObject
2788      9      300389   2017-08-03   d512d75856a89b7602c143046f9e030e112c81e5  SiPixelTemplateDBObject
2789      10     301046   2017-08-12   c624d964356bf70df4a33761c32f8976a0d89257  SiPixelTemplateDBObject
2790      11     302131   2017-08-31   ca6f56bb82434b69bd951d5ca89cdce841473ce0  SiPixelTemplateDBObject
2791      12     303790   2017-09-23   cbb7b82c563a21f06bf7e57bec3a490625c98d18  SiPixelTemplateDBObject
2792      13     303998   2017-09-27   7aec56e15aed26f7b686535ecaff59911947f802  SiPixelTemplateDBObject
2793      14     304911   2017-10-13   d5bc4c312735d2d97e617ae3d45c735639a1cd82  SiPixelTemplateDBObject
2794      15     313041   2018-03-28   f3752b4a1fac4aa65c6439ab75d216dcc3ed27a9  SiPixelTemplateDBObject 2018
2795      16     314881   2018-04-22   bc8784a015d78068c51c9a581328b064299a31b4  SiPixelTemplateDBObject
2796      17     316758   2018-05-22   3a8abe693d1d3df829212e1049330e68f3392282  SiPixelTemplateDBObject
2797      18     317475   2018-06-05   83e09d34c122040bca0f5daf4b89a0435663c230  SiPixelTemplateDBObject
2798      19     317485   2018-06-05   3a8abe693d1d3df829212e1049330e68f3392282  SiPixelTemplateDBObject
2799      20     317527   2018-06-06   356b03fcd7e869ef8f28153318f0cd32c27b1f40  SiPixelTemplateDBObject
2800      21     317661   2018-06-10   1cbb2fc3edf448c50d00fac3aa47c8cf3b19ac6f  SiPixelTemplateDBObject
2801      22     317664   2018-06-11   356b03fcd7e869ef8f28153318f0cd32c27b1f40  SiPixelTemplateDBObject
2802      23     318227   2018-06-21   1cbb2fc3edf448c50d00fac3aa47c8cf3b19ac6f  SiPixelTemplateDBObject
2803      24     320377   2018-07-26   5e9cd265167ec543aac49685cf706a58014c08f8  SiPixelTemplateDBObject
2804      25     321831   2018-08-27   e47d453934b47f06e53a323bd9ae16b38ec1bbd1  SiPixelTemplateDBObject
2805      26     322510   2018-09-09   df783de5f30a99bc036d8973e48da78513206aba  SiPixelTemplateDBObject
2806      27     322603   2018-09-09   e47d453934b47f06e53a323bd9ae16b38ec1bbd1  SiPixelTemplateDBObject
2807      28     323232   2018-09-21   9e660fd65e392f01f69cde77a52f59e934d7737d  SiPixelTemplateDBObject
2808      29     324245   2018-10-08   6caa50b30aa80ede330585888cdc5e804853d1da  SiPixelTemplateDBObject
2809   */
2810
2811   static const int nIOVs =
2812       29;  //   1       2       3       4       5       6       7       8       9      10      11      12      13      14      15
2813   int IOVboundaries[nIOVs] = {
2814       271866,
2815       276315,
2816       278271,
2817       280928,
2818       290543,
2819       297281,
2820       298653,
2821       299443,
2822       300389,
2823       301046,
2824       302131,
2825       303790,
2826       303998,
2827       304911,
2828       313041,
2829       //  16      17      18      19      20      21      22      23      24      25      26      27      28      29
2830       314881,
2831       316758,
2832       317475,
2833       317485,
2834       317527,
2835       317661,
2836       317664,
2837       318227,
2838       320377,
2839       321831,
2840       322510,
2841       322603,
2842       323232,
2843       324245};
2844
2845   int benchmarkruns[nIOVs] = {271866, 276315, 278271, 280928, 290543, 297281, 298653, 299443, 300389, 301046,
2846                               302131, 303790, 303998, 304911, 313041, 314881, 316758, 317475, 317485, 317527,
2847                               317661, 317664, 318227, 320377, 321831, 322510, 322603, 323232, 324245};
2848
2849   std::string dates[nIOVs] = {"2016-04-28", "2016-07-04", "2016-08-05", "2016-09-16", "2017-03-30", "2017-06-22",
2850                               "2017-07-09", "2017-07-19", "2017-08-03", "2017-08-12", "2017-08-31", "2017-09-23",
2851                               "2017-09-27", "2017-10-13", "2018-03-28", "2018-04-22", "2018-05-22", "2018-06-05",
2852                               "2018-06-05", "2018-06-06", "2018-06-10", "2018-06-11", "2018-06-21", "2018-07-26",
2853                               "2018-08-27", "2018-09-09", "2018-09-09", "2018-09-21", "2018-10-08"};
2854
2855   TArrow *IOV_lines[nIOVs];
2856   c->cd();
2857   c->Update();
2858
2859   TArrow *a_lines[nIOVs];
2860   TArrow *b_lines[nIOVs];
2861   for (Int_t IOV = 0; IOV < nIOVs; IOV++) {
2862     // check we are not in the RMS histogram to avoid first line
2863     if (IOVboundaries[IOV] < vruns.front())
2864       continue;
2865     //&& ((TString)c->GetName()).Contains("RMS")) continue;
2866     int closestrun = pv::closest(vruns, IOVboundaries[IOV]);
2867     int closestbenchmark = pv::closest(vruns, benchmarkruns[IOV]);
2868
2869     if (lumi_axis_format) {
2870       if (closestrun < 0)
2871         continue;
2872       //logInfo<< "natural boundary: " << IOVboundaries[IOV] << " closest:" << closestrun << std::endl;
2873
2874       a_lines[IOV] = new TArrow(
2875           lumiMapByRun.at(closestrun), (c->GetUymin()), lumiMapByRun.at(closestrun), c->GetUymax(), 0.5, "|>");
2876
2877       if (closestbenchmark < 0)
2878         continue;
2879       b_lines[IOV] = new TArrow(lumiMapByRun.at(closestbenchmark),
2880                                 (c->GetUymin()),
2881                                 lumiMapByRun.at(closestbenchmark),
2882                                 c->GetUymax(),
2883                                 0.5,
2884                                 "|>");
2885
2886     } else if (time_axis_format) {
2887       if (closestrun < 0)
2888         continue;
2889       logInfo << "natural boundary: " << IOVboundaries[IOV] << " closest:" << closestrun << std::endl;
2890       a_lines[IOV] = new TArrow(
2891           timeMap.at(closestrun).Convert(), (c->GetUymin()), timeMap.at(closestrun).Convert(), c->GetUymax(), 0.5, "|>");
2892
2893       if (closestbenchmark < 0)
2894         continue;
2895       b_lines[IOV] = new TArrow(timeMap.at(closestbenchmark).Convert(),
2896                                 (c->GetUymin()),
2897                                 timeMap.at(closestbenchmark).Convert(),
2898                                 c->GetUymax(),
2899                                 0.5,
2900                                 "|>");
2901
2902     } else {
2903       a_lines[IOV] = new TArrow(IOVboundaries[IOV],
2904                                 (c->GetUymin()),
2905                                 IOVboundaries[IOV],
2906                                 0.65 * c->GetUymax(),
2907                                 0.5,
2908                                 "|>");  //(c->GetUymin()+0.2*(c->GetUymax()-c->GetUymin()) ),0.5,"|>");
2909       b_lines[IOV] = new TArrow(benchmarkruns[IOV],
2910                                 (c->GetUymin()),
2911                                 benchmarkruns[IOV],
2912                                 0.65 * c->GetUymax(),
2913                                 0.5,
2914                                 "|>");  //(c->GetUymin()+0.2*(c->GetUymax()-c->GetUymin()) ),0.5,"|>");
2915     }
2916     a_lines[IOV]->SetLineColor(kBlue);
2917     a_lines[IOV]->SetLineStyle(9);
2918     a_lines[IOV]->SetLineWidth(1);
2919     a_lines[IOV]->Draw("same");
2920
2921     b_lines[IOV]->SetLineColor(kGray);
2922     b_lines[IOV]->SetLineStyle(1);
2923     b_lines[IOV]->SetLineWidth(2);
2924     //b_lines[IOV]->Draw("same");
2925   }
2926
2927   TPaveText *runnumbers[nIOVs];
2928
2929   for (Int_t IOV = 0; IOV < nIOVs; IOV++) {
2930     if (IOVboundaries[IOV] < vruns.front())
2931       continue;
2932     //&& ((TString)c->GetName()).Contains("RMS")) continue;
2933     int closestrun = pv::closest(vruns, IOVboundaries[IOV]);
2934
2935     Int_t ix1;
2936     Int_t ix2;
2937     Int_t iw = gPad->GetWw();
2938     Int_t ih = gPad->GetWh();
2939     Double_t x1p, y1p, x2p, y2p;
2940     gPad->GetPadPar(x1p, y1p, x2p, y2p);
2941     ix1 = (Int_t)(iw * x1p);
2942     ix2 = (Int_t)(iw * x2p);
2943     Double_t wndc = TMath::Min(1., (Double_t)iw / (Double_t)ih);
2944     Double_t rw = wndc / (Double_t)iw;
2945     Double_t x1ndc = (Double_t)ix1 * rw;
2946     Double_t x2ndc = (Double_t)ix2 * rw;
2947     Double_t rx1, ry1, rx2, ry2;
2948     gPad->GetRange(rx1, ry1, rx2, ry2);
2949     Double_t rx = (x2ndc - x1ndc) / (rx2 - rx1);
2950     Double_t _sx;
2951     if (lumi_axis_format) {
2952       if (closestrun < 0)
2953         break;
2954       _sx = rx * (lumiMapByRun.at(closestrun) - rx1) + x1ndc;  //-0.05;
2955     } else if (time_axis_format) {
2956       if (closestrun < 0)
2957         break;
2958       _sx = rx * (timeMap.at(closestrun).Convert() - rx1) + x1ndc;  //-0.05;
2959     } else {
2960       _sx = rx * (IOVboundaries[IOV] - rx1) + x1ndc;  //-0.05
2961     }
2962     Double_t _dx = _sx + 0.03;
2963
2964     Int_t index = IOV % 5;
2965     // if(IOV<5)
2966     //   index=IOV;
2967     // else{
2968     //   index=IOV-5;
2969     // }
2970
2971     //runnumbers[IOV] = new TPaveText(_sx+0.001,0.14+(0.03*index),_dx,(0.17+0.03*index),"blNDC");
2972     runnumbers[IOV] = new TPaveText(_sx + 0.001, 0.89 - (0.05 * index), _dx + 0.024, (0.92 - 0.05 * index), "blNDC");
2973     //runnumbers[IOV]->SetTextAlign(11);
2974     TText *textRun = runnumbers[IOV]->AddText(Form("%i", int(IOVboundaries[IOV])));
2975     //TText *textRun = runnumbers[IOV]->AddText(Form("%s",dates[IOV].c_str()));
2976     textRun->SetTextSize(0.038);
2977     textRun->SetTextColor(kBlue);
2978     runnumbers[IOV]->SetFillColor(10);
2979     runnumbers[IOV]->SetLineColor(kBlue);
2980     runnumbers[IOV]->SetBorderSize(1);
2981     runnumbers[IOV]->SetLineWidth(1);
2982     runnumbers[IOV]->SetTextColor(kBlue);
2983     runnumbers[IOV]->SetTextFont(42);
2984     if (drawText)
2985       runnumbers[IOV]->Draw("same");
2986   }
2987 }
2988
2989 /*! \fn processData
2990  *  \brief function where the magic happens, take the raw inputs and creates the output trends
2991  */
2992
2993 /*--------------------------------------------------------------------*/
2994 outTrends processData(size_t iter,
2995                       std::vector<int> intersection,
2996                       const Int_t nDirs_,
2997                       const char *dirs[10],
2998                       TString LegLabels[10],
2999                       bool useRMS)
3000 /*--------------------------------------------------------------------*/
3001 {
3002   outTrends ret;
3003   ret.init();
3004
3005   unsigned int effSize = std::min(nWorkers, intersection.size());
3006
3007   unsigned int pitch = std::floor(intersection.size() / effSize);
3008   unsigned int first = iter * pitch;
3009   unsigned int last = (iter == (effSize - 1)) ? intersection.size() : ((iter + 1) * pitch);
3010
3011   logInfo << "iter:" << iter << "| pitch: " << pitch << " [" << first << "-" << last << ")" << std::endl;
3012
3013   ret.m_index = iter;
3014
3015   for (unsigned int n = first; n < last; n++) {
3016     //in case of debug, use only 50
3017     //for(unsigned int n=0; n<50;n++){
3018
3019     //if(intersection.at(n)!=283946)
3020     //  continue;
3021
3022     if (VERBOSE) {
3023       logInfo << "iter: " << iter << " " << n << " " << intersection.at(n) << std::endl;
3024     }
3025
3026     TFile *fins[nDirs_];
3027
3028     TH1F *dxyPhiMeanTrend[nDirs_];
3029     TH1F *dxyPhiWidthTrend[nDirs_];
3030     TH1F *dzPhiMeanTrend[nDirs_];
3031     TH1F *dzPhiWidthTrend[nDirs_];
3032
3033     TH1F *dxyLadderMeanTrend[nDirs_];
3034     TH1F *dxyLadderWidthTrend[nDirs_];
3035     TH1F *dzLadderWidthTrend[nDirs_];
3036     TH1F *dzLadderMeanTrend[nDirs_];
3037
3038     TH1F *dxyModZMeanTrend[nDirs_];
3039     TH1F *dxyModZWidthTrend[nDirs_];
3040     TH1F *dzModZMeanTrend[nDirs_];
3041     TH1F *dzModZWidthTrend[nDirs_];
3042
3043     TH1F *dxyEtaMeanTrend[nDirs_];
3044     TH1F *dxyEtaWidthTrend[nDirs_];
3045     TH1F *dzEtaMeanTrend[nDirs_];
3046     TH1F *dzEtaWidthTrend[nDirs_];
3047
3048     TH1F *dxyNormPhiWidthTrend[nDirs_];
3049     TH1F *dxyNormEtaWidthTrend[nDirs_];
3050     TH1F *dzNormPhiWidthTrend[nDirs_];
3051     TH1F *dzNormEtaWidthTrend[nDirs_];
3052
3053     TH1F *dxyNormPtWidthTrend[nDirs_];
3054     TH1F *dzNormPtWidthTrend[nDirs_];
3055     TH1F *dxyPtWidthTrend[nDirs_];
3056     TH1F *dzPtWidthTrend[nDirs_];
3057
3058     TH1F *dxyIntegralTrend[nDirs_];
3059     TH1F *dzIntegralTrend[nDirs_];
3060
3061     bool areAllFilesOK = true;
3062     Int_t lastOpen = 0;
3063
3064     // loop over the objects
3065     for (Int_t j = 0; j < nDirs_; j++) {
3066       //fins[j] = TFile::Open(Form("%s/PVValidation_%s_%i.root",dirs[j],dirs[j],intersection[n]));
3067
3068       size_t position = std::string(dirs[j]).find("/");
3069       string stem = std::string(dirs[j]).substr(position + 1);  // get from position to the end
3070
3071       fins[j] = new TFile(Form("%s/PVValidation_%s_%i.root", dirs[j], stem.c_str(), intersection[n]));
3072       if (fins[j]->IsZombie()) {
3073         logError << Form("%s/PVValidation_%s_%i.root", dirs[j], stem.c_str(), intersection[n])
3074                  << " is a Zombie! cannot combine" << std::endl;
3075         areAllFilesOK = false;
3076         lastOpen = j;
3077         break;
3078       }
3079
3080       if (VERBOSE) {
3081         logInfo << Form("%s/PVValidation_%s_%i.root", dirs[j], stem.c_str(), intersection[n])
3082                 << " has size: " << fins[j]->GetSize() << " b ";
3083       }
3084
3085       // sanity check
3086       TH1F *h_tracks = (TH1F *)fins[j]->Get("PVValidation/EventFeatures/h_nTracks");
3087       if (j == 0) {
3088         TH1F *h_lumi = (TH1F *)fins[j]->Get("PVValidation/EventFeatures/h_lumiFromConfig");
3089         double lumi = h_lumi->GetBinContent(1);
3090         ret.m_lumiSoFar += lumi;
3091         //logInfo<<"lumi: "<<lumi
3092         //   <<" ,lumi so far: "<<lumiSoFar<<std::endl;
3093
3094         //outfile<<"run "<<intersection[n]<<" lumi: "<<lumi
3095         //     <<" ,lumi so far: "<<lumiSoFar<<std::endl;
3096       }
3097
3098       Double_t numEvents = h_tracks->GetEntries();
3099
3100       if (numEvents < 2500) {
3101         logWarning << "excluding run " << intersection[n] << " because it has less than 2.5k events" << std::endl;
3102         areAllFilesOK = false;
3103         lastOpen = j;
3104         break;
3105       }
3106
3107       dxyPhiMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dxy_phi");
3108       //dxyPhiMeanTrend[j]     = checkTH1AndReturn(fins[j],"PVValidation/MeanTrends/means_dxy_phi");
3109       dxyPhiWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_phi");
3110       dzPhiMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dz_phi");
3111       dzPhiWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_phi");
3112
3113       dxyLadderMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dxy_ladder");
3114       dxyLadderWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_ladder");
3115       dzLadderMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dz_ladder");
3116       dzLadderWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_ladder");
3117
3118       dxyEtaMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dxy_eta");
3119       dxyEtaWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_eta");
3120       dzEtaMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dz_eta");
3121       dzEtaWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_eta");
3122
3123       dxyModZMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dxy_modZ");
3124       dxyModZWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_modZ");
3125       dzModZMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dz_modZ");
3126       dzModZWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_modZ");
3127
3128       dxyNormPhiWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dxy_phi");
3129       dxyNormEtaWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dxy_eta");
3130       dzNormPhiWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dz_phi");
3131       dzNormEtaWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dz_eta");
3132
3133       dxyNormPtWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dxy_pTCentral");
3134       dzNormPtWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dz_pTCentral");
3135       dxyPtWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_pTCentral");
3136       dzPtWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_pTCentral");
3137
3138       dxyIntegralTrend[j] = (TH1F *)fins[j]->Get("PVValidation/ProbeTrackFeatures/h_probedxyRefitV");
3139       dzIntegralTrend[j] = (TH1F *)fins[j]->Get("PVValidation/ProbeTrackFeatures/h_probedzRefitV");
3140
3141       // fill the vectors of biases
3142
3143       auto dxyPhiBiases = getBiases(dxyPhiMeanTrend[j]);
3144
3145       //logInfo<<"\n" <<j<<" "<< LegLabels[j] << " dxy(phi) mean: "<< dxyPhiBiases.getWeightedMean()
3146       //       <<" dxy(phi) max: "<< dxyPhiBiases.getMax()
3147       //       <<" dxy(phi) min: "<< dxyPhiBiases.getMin()
3148       //       << std::endl;
3149
3150       ret.m_dxyPhiMeans[LegLabels[j]].push_back(dxyPhiBiases.getWeightedMean());
3151       ret.m_dxyPhiChi2[LegLabels[j]].push_back(TMath::Log10(dxyPhiBiases.getNormChi2()));
3152       ret.m_dxyPhiKS[LegLabels[j]].push_back(dxyPhiBiases.getKSScore());
3153
3154       //logInfo<<"\n" <<j<<" "<< LegLabels[j] << " dxy(phi) ks score: "<< dxyPhiBiases.getKSScore() << std::endl;
3155
3156       useRMS
3157           ? ret.m_dxyPhiLo[LegLabels[j]].push_back(dxyPhiBiases.getWeightedMean() - 2 * dxyPhiBiases.getWeightedRMS())
3158           : ret.m_dxyPhiLo[LegLabels[j]].push_back(dxyPhiBiases.getMin());
3159       useRMS
3160           ? ret.m_dxyPhiHi[LegLabels[j]].push_back(dxyPhiBiases.getWeightedMean() + 2 * dxyPhiBiases.getWeightedRMS())
3161           : ret.m_dxyPhiHi[LegLabels[j]].push_back(dxyPhiBiases.getMax());
3162
3163       auto dxyEtaBiases = getBiases(dxyEtaMeanTrend[j]);
3164       ret.m_dxyEtaMeans[LegLabels[j]].push_back(dxyEtaBiases.getWeightedMean());
3165       ret.m_dxyEtaChi2[LegLabels[j]].push_back(TMath::Log10(dxyEtaBiases.getNormChi2()));
3166       ret.m_dxyEtaKS[LegLabels[j]].push_back(dxyEtaBiases.getKSScore());
3167       useRMS
3168           ? ret.m_dxyEtaLo[LegLabels[j]].push_back(dxyEtaBiases.getWeightedMean() - 2 * dxyEtaBiases.getWeightedRMS())
3169           : ret.m_dxyEtaLo[LegLabels[j]].push_back(dxyEtaBiases.getMin());
3170       useRMS
3171           ? ret.m_dxyEtaHi[LegLabels[j]].push_back(dxyEtaBiases.getWeightedMean() + 2 * dxyEtaBiases.getWeightedRMS())
3172           : ret.m_dxyEtaHi[LegLabels[j]].push_back(dxyEtaBiases.getMax());
3173
3174       auto dzPhiBiases = getBiases(dzPhiMeanTrend[j]);
3175       ret.m_dzPhiMeans[LegLabels[j]].push_back(dzPhiBiases.getWeightedMean());
3176       ret.m_dzPhiChi2[LegLabels[j]].push_back(TMath::Log10(dzPhiBiases.getNormChi2()));
3177       ret.m_dzPhiKS[LegLabels[j]].push_back(dzPhiBiases.getKSScore());
3178       useRMS ? ret.m_dzPhiLo[LegLabels[j]].push_back(dzPhiBiases.getWeightedMean() - 2 * dzPhiBiases.getWeightedRMS())
3179              : ret.m_dzPhiLo[LegLabels[j]].push_back(dzPhiBiases.getMin());
3180       useRMS ? ret.m_dzPhiHi[LegLabels[j]].push_back(dzPhiBiases.getWeightedMean() + 2 * dzPhiBiases.getWeightedRMS())
3181              : ret.m_dzPhiHi[LegLabels[j]].push_back(dzPhiBiases.getMax());
3182
3183       auto dzEtaBiases = getBiases(dzEtaMeanTrend[j]);
3184       ret.m_dzEtaMeans[LegLabels[j]].push_back(dzEtaBiases.getWeightedMean());
3185       ret.m_dzEtaChi2[LegLabels[j]].push_back(TMath::Log10(dzEtaBiases.getNormChi2()));
3186       ret.m_dzEtaKS[LegLabels[j]].push_back(dzEtaBiases.getKSScore());
3187       useRMS ? ret.m_dzEtaLo[LegLabels[j]].push_back(dzEtaBiases.getWeightedMean() - 2 * dzEtaBiases.getWeightedRMS())
3188              : ret.m_dzEtaLo[LegLabels[j]].push_back(dzEtaBiases.getMin());
3189       useRMS ? ret.m_dzEtaHi[LegLabels[j]].push_back(dzEtaBiases.getWeightedMean() + 2 * dzEtaBiases.getWeightedRMS())
3190              : ret.m_dzEtaHi[LegLabels[j]].push_back(dzEtaBiases.getMax());
3191
3192       // unrolled histograms
3193       ret.m_dxyVect[LegLabels[j]].push_back(getUnrolledHisto(dxyIntegralTrend[j]));
3194       ret.m_dzVect[LegLabels[j]].push_back(getUnrolledHisto(dzIntegralTrend[j]));
3195
3196       //logInfo<<std::endl;
3197       //logInfo<<" n. bins: "<< dxyVect[LegLabels[j]].back().get_n_bins()
3198       //       <<" y-min:   "<< dxyVect[LegLabels[j]].back().get_y_min()
3199       //       <<" y-max:   "<< dxyVect[LegLabels[j]].back().get_y_max() << std::endl;
3200
3201       // beautify the histograms
3202       MakeNiceTrendPlotStyle(dxyPhiMeanTrend[j], j);
3203       MakeNiceTrendPlotStyle(dxyPhiWidthTrend[j], j);
3204       MakeNiceTrendPlotStyle(dzPhiMeanTrend[j], j);
3205       MakeNiceTrendPlotStyle(dzPhiWidthTrend[j], j);
3206
3207       MakeNiceTrendPlotStyle(dxyLadderMeanTrend[j], j);
3208       MakeNiceTrendPlotStyle(dxyLadderWidthTrend[j], j);
3209       MakeNiceTrendPlotStyle(dzLadderMeanTrend[j], j);
3210       MakeNiceTrendPlotStyle(dzLadderWidthTrend[j], j);
3211
3212       MakeNiceTrendPlotStyle(dxyEtaMeanTrend[j], j);
3213       MakeNiceTrendPlotStyle(dxyEtaWidthTrend[j], j);
3214       MakeNiceTrendPlotStyle(dzEtaMeanTrend[j], j);
3215       MakeNiceTrendPlotStyle(dzEtaWidthTrend[j], j);
3216
3217       MakeNiceTrendPlotStyle(dxyModZMeanTrend[j], j);
3218       MakeNiceTrendPlotStyle(dxyModZWidthTrend[j], j);
3219       MakeNiceTrendPlotStyle(dzModZMeanTrend[j], j);
3220       MakeNiceTrendPlotStyle(dzModZWidthTrend[j], j);
3221
3222       MakeNiceTrendPlotStyle(dxyNormPhiWidthTrend[j], j);
3223       MakeNiceTrendPlotStyle(dxyNormEtaWidthTrend[j], j);
3224       MakeNiceTrendPlotStyle(dzNormPhiWidthTrend[j], j);
3225       MakeNiceTrendPlotStyle(dzNormEtaWidthTrend[j], j);
3226
3227       MakeNiceTrendPlotStyle(dxyNormPtWidthTrend[j], j);
3228       MakeNiceTrendPlotStyle(dzNormPtWidthTrend[j], j);
3229       MakeNiceTrendPlotStyle(dxyPtWidthTrend[j], j);
3230       MakeNiceTrendPlotStyle(dzPtWidthTrend[j], j);
3231     }
3232
3233     if (!areAllFilesOK) {
3234       // do all the necessary deletions
3235       logWarning << "====> not all files are OK" << std::endl;
3236
3237       for (int i = 0; i < lastOpen; i++) {
3238         fins[i]->Close();
3239       }
3240       continue;
3241     } else {
3242       ret.m_runs.push_back(intersection.at(n));
3243       // push back the vector of lumi (in fb at this point)
3244       ret.m_lumiByRun.push_back(ret.m_lumiSoFar / 1000.);
3245       ret.m_lumiMapByRun[intersection.at(n)] = ret.m_lumiSoFar / 1000.;
3246     }
3247
3248     if (VERBOSE) {
3249       logInfo << "I am still here - runs.size(): " << ret.m_runs.size() << std::endl;
3250     }
3251
3252     // Bias plots
3253
3254     TCanvas *BiasesCanvas = new TCanvas(Form("Biases_%i", intersection.at(n)), "Biases", 1200, 1200);
3255     arrangeOutCanvas(BiasesCanvas,
3256                      dxyPhiMeanTrend,
3257                      dzPhiMeanTrend,
3258                      dxyEtaMeanTrend,
3259                      dzEtaMeanTrend,
3260                      nDirs_,
3261                      LegLabels,
3262                      intersection.at(n));
3263
3264     BiasesCanvas->SaveAs(Form("Biases_%i.pdf", intersection.at(n)));
3265     BiasesCanvas->SaveAs(Form("Biases_%i.png", intersection.at(n)));
3266
3267     // Bias vs L1 modules position
3268
3269     TCanvas *BiasesL1Canvas = new TCanvas(Form("BiasesL1_%i", intersection.at(n)), "BiasesL1", 1200, 1200);
3270     arrangeOutCanvas(BiasesL1Canvas,
3271                      dxyLadderMeanTrend,
3272                      dzLadderMeanTrend,
3273                      dxyModZMeanTrend,
3274                      dzModZMeanTrend,
3275                      nDirs_,
3276                      LegLabels,
3277                      intersection.at(n));
3278
3279     BiasesL1Canvas->SaveAs(Form("BiasesL1_%i.pdf", intersection.at(n)));
3280     BiasesL1Canvas->SaveAs(Form("BiasesL1_%i.png", intersection.at(n)));
3281
3282     // Resolution plots
3283
3284     TCanvas *ResolutionsCanvas = new TCanvas(Form("Resolutions_%i", intersection.at(n)), "Resolutions", 1200, 1200);
3285     arrangeOutCanvas(ResolutionsCanvas,
3286                      dxyPhiWidthTrend,
3287                      dzPhiWidthTrend,
3288                      dxyEtaWidthTrend,
3289                      dzEtaWidthTrend,
3290                      nDirs_,
3291                      LegLabels,
3292                      intersection.at(n));
3293
3294     ResolutionsCanvas->SaveAs(Form("Resolutions_%i.pdf", intersection.at(n)));
3295     ResolutionsCanvas->SaveAs(Form("Resolutions_%i.png", intersection.at(n)));
3296
3297     // Resolution plots vs L1 modules position
3298
3299     TCanvas *ResolutionsL1Canvas = new TCanvas(Form("ResolutionsL1_%i", intersection.at(n)), "Resolutions", 1200, 1200);
3300     arrangeOutCanvas(ResolutionsL1Canvas,
3301                      dxyLadderWidthTrend,
3302                      dzLadderWidthTrend,
3303                      dxyModZWidthTrend,
3304                      dzModZWidthTrend,
3305                      nDirs_,
3306                      LegLabels,
3307                      intersection.at(n));
3308
3309     ResolutionsL1Canvas->SaveAs(Form("ResolutionsL1_%i.pdf", intersection.at(n)));
3310     ResolutionsL1Canvas->SaveAs(Form("ResolutionsL1_%i.png", intersection.at(n)));
3311
3312     // Pull plots
3313
3314     TCanvas *PullsCanvas = new TCanvas(Form("Pulls_%i", intersection.at(n)), "Pulls", 1200, 1200);
3315     arrangeOutCanvas(PullsCanvas,
3316                      dxyNormPhiWidthTrend,
3317                      dzNormPhiWidthTrend,
3318                      dxyNormEtaWidthTrend,
3319                      dzNormEtaWidthTrend,
3320                      nDirs_,
3321                      LegLabels,
3322                      intersection.at(n));
3323
3324     PullsCanvas->SaveAs(Form("Pulls_%i.pdf", intersection.at(n)));
3325     PullsCanvas->SaveAs(Form("Pulls_%i.png", intersection.at(n)));
3326
3327     // pT plots
3328
3329     TCanvas *ResolutionsVsPt =
3330         new TCanvas(Form("ResolutionsVsPT_%i", intersection.at(n)), "ResolutionsVsPt", 1200, 1200);
3331     arrangeOutCanvas(ResolutionsVsPt,
3332                      dxyPtWidthTrend,
3333                      dzPtWidthTrend,
3334                      dxyNormPtWidthTrend,
3335                      dzNormPtWidthTrend,
3336                      nDirs_,
3337                      LegLabels,
3338                      intersection.at(n));
3339
3340     ResolutionsVsPt->SaveAs(Form("ResolutionsVsPt_%i.pdf", intersection.at(n)));
3341     ResolutionsVsPt->SaveAs(Form("ResolutionsVsPt_%i.png", intersection.at(n)));
3342
3343     // do all the necessary deletions
3344
3345     for (int i = 0; i < nDirs_; i++) {
3346       delete dxyPhiMeanTrend[i];
3347       delete dzPhiMeanTrend[i];
3348       delete dxyEtaMeanTrend[i];
3349       delete dzEtaMeanTrend[i];
3350
3351       delete dxyPhiWidthTrend[i];
3352       delete dzPhiWidthTrend[i];
3353       delete dxyEtaWidthTrend[i];
3354       delete dzEtaWidthTrend[i];
3355
3356       delete dxyNormPhiWidthTrend[i];
3357       delete dxyNormEtaWidthTrend[i];
3358       delete dzNormPhiWidthTrend[i];
3359       delete dzNormEtaWidthTrend[i];
3360
3361       delete dxyNormPtWidthTrend[i];
3362       delete dzNormPtWidthTrend[i];
3363       delete dxyPtWidthTrend[i];
3364       delete dzPtWidthTrend[i];
3365
3366       fins[i]->Close();
3367     }
3368
3369     delete BiasesCanvas;
3370     delete BiasesL1Canvas;
3371     delete ResolutionsCanvas;
3372     delete ResolutionsL1Canvas;
3373     delete PullsCanvas;
3374     delete ResolutionsVsPt;
3375
3376     if (VERBOSE) {
3377       logInfo << std::endl;
3378     }
3379   }
3380
3381   return ret;
3382 }