Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-11-25 03:29:09 

0001 #include <string>
0002 #include <sstream>
0003 #include <vector>
0004 #include <iomanip>
0005 #include <stdio.h>
0006 #include <iostream>
0007 #include <fstream>
0008 #include <cassert>
0009 #include "TFile.h"
0010 #include "TPaveStats.h"
0011 #include "TROOT.h"
0012 #include "TList.h"
0013 #include "TNtuple.h"
0014 #include "TTree.h"
0015 #include "TError.h"
0016 #include "TH1.h"
0017 #include "TArrow.h"
0018 #include "TH2.h"
0019 #include "THStack.h"
0020 #include "TStyle.h"
0021 #include "TLegendEntry.h"
0022 #include "TPaveText.h"
0023 #include "TCut.h"
0024 #include "TLegend.h"
0025 #include "TGraphErrors.h"
0026 #include "TF1.h"
0027 #include "TMath.h"
0028 #include "TVectorD.h"
0029 #include "TCanvas.h"
0030 #include "TColor.h"
0031 #include "TAxis.h"
0032 #include "TGaxis.h"
0033 #include "TROOT.h"
0034 #include "TObjArray.h"
0035 #include "TGraphErrors.h"
0036 #include "TF1.h"
0037 #include "TMinuit.h"
0038 #include "TString.h"
0039 #include "TMath.h"
0040 #include <TDatime.h>
0041 #include <TSpectrum.h>
0042 #include <TSystem.h>
0043 #include <TTimeStamp.h>
0044 #include <TStopwatch.h>
0045 #include <TObjString.h>
0046 //#include "Alignment/OfflineValidation/macros/TkAlStyle.cc"
0047 #include "Alignment/OfflineValidation/macros/CMS_lumi.h"
0048 #define PLOTTING_MACRO  // to remove message logger
0049 #include "Alignment/OfflineValidation/interface/PVValidationHelpers.h"
0050 
0051 /* 
0052    This is an auxilliary class to store the list of files
0053    to be used to plot
0054 */
0055 
0056 class PVValidationVariables {
0057 public:
0058   PVValidationVariables(TString fileName, TString baseDir, TString legName = "", int color = 1, int style = 1);
0059   int getLineColor() { return lineColor; }
0060   int getMarkerStyle() { return markerStyle; }
0061   int getLineStyle() { return lineStyle; }
0062   TString getName() { return legendName; }
0063   TFile *getFile() { return file; }
0064   TString getFileName() { return fname; }
0065 
0066 private:
0067   TFile *file;
0068   int lineColor;
0069   int lineStyle;
0070   int markerStyle;
0071   TString legendName;
0072   TString fname;
0073 };
0074 
0075 PVValidationVariables::PVValidationVariables(
0076     TString fileName, TString baseDir, TString legName, int lColor, int lStyle) {
0077   fname = fileName;
0078   lineColor = lColor;
0079 
0080   int ndigits = 1 + std::floor(std::log10(lStyle));
0081   if (ndigits == 4) {
0082     lineStyle = lStyle % 100;
0083     markerStyle = lStyle / 100;
0084   } else {
0085     lineStyle = 1;
0086     markerStyle = lStyle;
0087   }
0088 
0089   if (legName == "") {
0090     std::string s_fileName = fileName.Data();
0091     int start = 0;
0092     if (s_fileName.find('/'))
0093       start = s_fileName.find_last_of('/') + 1;
0094     int stop = s_fileName.find_last_of('.');
0095     legendName = s_fileName.substr(start, stop - start);
0096   } else {
0097     legendName = legName;
0098   }
0099 
0100   // check if the base dir exists
0101   file = TFile::Open(fileName.Data(), "READ");
0102 
0103   if (!file) {
0104     std::cout << "ERROR! file " << fileName.Data() << " does not exist!" << std::endl;
0105     assert(false);
0106   }
0107 
0108   if (file->Get(baseDir.Data())) {
0109     std::cout << "found base directory: " << baseDir.Data() << std::endl;
0110   } else {
0111     std::cout << "no directory named: " << baseDir.Data() << std::endl;
0112     assert(false);
0113   }
0114 }
0115 
0116 /*
0117   This is an auxilliary enum and typedef used to handle the fit parameter types
0118 */
0119 
0120 namespace params {
0121 
0122   typedef std::pair<Double_t, Double_t> measurement;
0123 
0124   enum estimator { MEAN = 1, WIDTH = 2, MEDIAN = 3, MAD = 4, UNKWN = -1 };
0125 }  // namespace params
0126 
0127 /*
0128   This is an auxilliary struct used to handle the plot limits
0129 */
0130 
0131 struct Limits {
0132   // initializers list
0133 
0134   Limits()
0135       : _m_dxyPhiMax(80.),
0136         _m_dzPhiMax(80.),
0137         _m_dxyEtaMax(80.),
0138         _m_dzEtaMax(80.),
0139         _m_dxyPhiNormMax(0.5),
0140         _m_dzPhiNormMax(0.5),
0141         _m_dxyEtaNormMax(0.5),
0142         _m_dzEtaNormMax(0.5),
0143         _w_dxyPhiMax(120.),
0144         _w_dzPhiMax(180.),
0145         _w_dxyEtaMax(120.),
0146         _w_dzEtaMax(1000.),
0147         _w_dxyPhiNormMax(2.0),
0148         _w_dzPhiNormMax(2.0),
0149         _w_dxyEtaNormMax(2.0),
0150         _w_dzEtaNormMax(2.0) {}
0151 
0152   // getter methods
0153 
0154   std::pair<float, float> get_dxyPhiMax() const {
0155     std::pair<float, float> res(_m_dxyPhiMax, _w_dxyPhiMax);
0156     return res;
0157   }
0158 
0159   std::pair<float, float> get_dzPhiMax() const {
0160     std::pair<float, float> res(_m_dzPhiMax, _w_dzPhiMax);
0161     return res;
0162   }
0163 
0164   std::pair<float, float> get_dxyEtaMax() const {
0165     std::pair<float, float> res(_m_dxyEtaMax, _w_dxyEtaMax);
0166     return res;
0167   }
0168 
0169   std::pair<float, float> get_dzEtaMax() const {
0170     std::pair<float, float> res(_m_dzEtaMax, _w_dzEtaMax);
0171     return res;
0172   }
0173 
0174   std::pair<float, float> get_dxyPhiNormMax() const {
0175     std::pair<float, float> res(_m_dxyPhiNormMax, _w_dxyPhiNormMax);
0176     return res;
0177   }
0178 
0179   std::pair<float, float> get_dzPhiNormMax() const {
0180     std::pair<float, float> res(_m_dzPhiNormMax, _w_dzPhiNormMax);
0181     return res;
0182   }
0183 
0184   std::pair<float, float> get_dxyEtaNormMax() const {
0185     std::pair<float, float> res(_m_dxyEtaNormMax, _w_dxyEtaNormMax);
0186     return res;
0187   }
0188 
0189   std::pair<float, float> get_dzEtaNormMax() const {
0190     std::pair<float, float> res(_m_dzEtaNormMax, _w_dzEtaNormMax);
0191     return res;
0192   }
0193 
0194   // initializes to different values, if needed
0195 
0196   void init(float m_dxyPhiMax,
0197             float m_dzPhiMax,
0198             float m_dxyEtaMax,
0199             float m_dzEtaMax,
0200             float m_dxyPhiNormMax,
0201             float m_dzPhiNormMax,
0202             float m_dxyEtaNormMax,
0203             float m_dzEtaNormMax,
0204             float w_dxyPhiMax,
0205             float w_dzPhiMax,
0206             float w_dxyEtaMax,
0207             float w_dzEtaMax,
0208             float w_dxyPhiNormMax,
0209             float w_dzPhiNormMax,
0210             float w_dxyEtaNormMax,
0211             float w_dzEtaNormMax) {
0212     _m_dxyPhiMax = m_dxyPhiMax;
0213     _m_dzPhiMax = m_dzPhiMax;
0214     _m_dxyEtaMax = m_dxyEtaMax;
0215     _m_dzEtaMax = m_dzEtaMax;
0216     _m_dxyPhiNormMax = m_dxyPhiNormMax;
0217     _m_dzPhiNormMax = m_dzPhiNormMax;
0218     _m_dxyEtaNormMax = m_dxyEtaNormMax;
0219     _m_dzEtaNormMax = m_dzEtaNormMax;
0220     _w_dxyPhiMax = w_dxyPhiMax;
0221     _w_dzPhiMax = w_dzPhiMax;
0222     _w_dxyEtaMax = w_dxyEtaMax;
0223     _w_dzEtaMax = w_dzEtaMax;
0224     _w_dxyPhiNormMax = w_dxyPhiNormMax;
0225     _w_dzPhiNormMax = w_dzPhiNormMax;
0226     _w_dxyEtaNormMax = w_dxyEtaNormMax;
0227     _w_dzEtaNormMax = w_dzEtaNormMax;
0228   }
0229 
0230   void printAll() {
0231     std::cout << "======================================================" << std::endl;
0232     std::cout << "  The y-axis ranges on the plots will be the following:      " << std::endl;
0233 
0234     std::cout << "  mean of dxy vs Phi:         " << _m_dxyPhiMax << std::endl;
0235     std::cout << "  mean of dz  vs Phi:         " << _m_dzPhiMax << std::endl;
0236     std::cout << "  mean of dxy vs Eta:         " << _m_dxyEtaMax << std::endl;
0237     std::cout << "  mean of dz  vs Eta:         " << _m_dzEtaMax << std::endl;
0238 
0239     std::cout << "  mean of dxy vs Phi (norm):  " << _m_dxyPhiNormMax << std::endl;
0240     std::cout << "  mean of dz  vs Phi (norm):  " << _m_dzPhiNormMax << std::endl;
0241     std::cout << "  mean of dxy vs Eta (norm):  " << _m_dxyEtaNormMax << std::endl;
0242     std::cout << "  mean of dz  vs Eta (norm):  " << _m_dzEtaNormMax << std::endl;
0243 
0244     std::cout << "  width of dxy vs Phi:        " << _w_dxyPhiMax << std::endl;
0245     std::cout << "  width of dz  vs Phi:        " << _w_dzPhiMax << std::endl;
0246     std::cout << "  width of dxy vs Eta:        " << _w_dxyEtaMax << std::endl;
0247     std::cout << "  width of dz  vs Eta:        " << _w_dzEtaMax << std::endl;
0248 
0249     std::cout << "  width of dxy vs Phi (norm): " << _w_dxyPhiNormMax << std::endl;
0250     std::cout << "  width of dz  vs Phi (norm): " << _w_dzPhiNormMax << std::endl;
0251     std::cout << "  width of dxy vs Eta (norm): " << _w_dxyEtaNormMax << std::endl;
0252     std::cout << "  width of dz  vs Eta (norm): " << _w_dzEtaNormMax << std::endl;
0253 
0254     std::cout << "======================================================" << std::endl;
0255   }
0256 
0257 private:
0258   float _m_dxyPhiMax;
0259   float _m_dzPhiMax;
0260   float _m_dxyEtaMax;
0261   float _m_dzEtaMax;
0262   float _m_dxyPhiNormMax;
0263   float _m_dzPhiNormMax;
0264   float _m_dxyEtaNormMax;
0265   float _m_dzEtaNormMax;
0266 
0267   float _w_dxyPhiMax;
0268   float _w_dzPhiMax;
0269   float _w_dxyEtaMax;
0270   float _w_dzEtaMax;
0271   float _w_dxyPhiNormMax;
0272   float _w_dzPhiNormMax;
0273   float _w_dxyEtaNormMax;
0274   float _w_dzEtaNormMax;
0275 };
0276 
0277 Limits *thePlotLimits = new Limits();
0278 std::vector<PVValidationVariables *> sourceList;
0279 
0280 #define ARRAY_SIZE(array) (sizeof((array)) / sizeof((array[0])))
0281 
0282 void arrangeCanvas(TCanvas *canv,
0283                    TH1F *meanplots[100],
0284                    TH1F *widthplots[100],
0285                    Int_t nFiles,
0286                    TString LegLabels[10],
0287                    TString theDate = "bogus",
0288                    bool onlyBias = false,
0289                    bool setAutoLimits = true);
0290 void arrangeCanvas2D(TCanvas *canv,
0291                      TH2F *meanmaps[100],
0292                      TH2F *widthmaps[100],
0293                      Int_t nFiles,
0294                      TString LegLabels[10],
0295                      TString theDate = "bogus");
0296 void arrangeFitCanvas(
0297     TCanvas *canv, TH1F *meanplots[100], Int_t nFiles, TString LegLabels[10], TString theDate = "bogus");
0298 
0299 void arrangeBiasCanvas(TCanvas *canv,
0300                        TH1F *dxyPhiMeanTrend[100],
0301                        TH1F *dzPhiMeanTrend[100],
0302                        TH1F *dxyEtaMeanTrend[100],
0303                        TH1F *dzEtaMeanTrend[100],
0304                        Int_t nFiles,
0305                        TString LegLabels[10],
0306                        TString theDate = "bogus",
0307                        bool setAutoLimits = true);
0308 
0309 params::measurement getMedian(TH1F *histo);
0310 params::measurement getMAD(TH1F *histo);
0311 
0312 std::pair<params::measurement, params::measurement> fitResiduals(TH1 *hist, bool singleTime = false);
0313 
0314 Double_t DoubleSidedCB(double *x, double *par);
0315 std::pair<params::measurement, params::measurement> fitResidualsCB(TH1 *hist);
0316 
0317 Double_t tp0Fit(Double_t *x, Double_t *par5);
0318 std::pair<params::measurement, params::measurement> fitStudentTResiduals(TH1 *hist);
0319 
0320 void FillTrendPlot(TH1F *trendPlot, TH1F *residualsPlot[100], params::estimator firPar_, TString var_, Int_t nbins);
0321 
0322 // global (here for the moment)
0323 Int_t nBins_ = 48;
0324 void FillMap(TH2F *trendMap,
0325              std::vector<std::vector<TH1F *> > residualsMapPlot,
0326              params::estimator fitPar_,
0327              const int nBinsX = nBins_,
0328              const int nBinsY = nBins_);
0329 
0330 std::pair<TH2F *, TH2F *> trimTheMap(TH2 *hist);
0331 
0332 void MakeNiceTrendPlotStyle(TH1 *hist, Int_t color, Int_t style);
0333 void MakeNicePlotStyle(TH1 *hist);
0334 void MakeNiceMapStyle(TH2 *hist);
0335 void MakeNiceTF1Style(TF1 *f1, Int_t color);
0336 
0337 void FitPVResiduals(TString namesandlabels,
0338                     bool stdres = true,
0339                     bool do2DMaps = false,
0340                     TString theDate = "bogus",
0341                     bool setAutoLimits = true);
0342 TH1F *DrawZero(TH1F *hist, Int_t nbins, Double_t lowedge, Double_t highedge, Int_t iter);
0343 TH1F *DrawConstant(TH1F *hist, Int_t nbins, Double_t lowedge, Double_t highedge, Int_t iter, Double_t theConst);
0344 void makeNewXAxis(TH1F *h);
0345 void makeNewPairOfAxes(TH2F *h);
0346 
0347 // ancillary fitting functions
0348 Double_t fULine(Double_t *x, Double_t *par);
0349 Double_t fDLine(Double_t *x, Double_t *par);
0350 void FitULine(TH1 *hist);
0351 void FitDLine(TH1 *hist);
0352 
0353 params::measurement getTheRangeUser(TH1F *thePlot, Limits *thePlotLimits, bool tag = false);
0354 
0355 void setStyle();
0356 
0357 // global variables
0358 
0359 std::ofstream outfile("FittedDeltaZ.txt");
0360 
0361 // use the maximum of the three supported phases
0362 Int_t nLadders_ = 20;
0363 Int_t nModZ_ = 9;
0364 
0365 const Int_t nPtBins_ = 48;
0366 Float_t _boundMin = -0.5;
0367 Float_t _boundSx = (nBins_ / 4.) - 0.5;
0368 Float_t _boundDx = 3 * (nBins_ / 4.) - 0.5;
0369 Float_t _boundMax = nBins_ - 0.5;
0370 Float_t etaRange = 2.5;
0371 Float_t minPt_ = 1.;
0372 Float_t maxPt_ = 20.;
0373 bool isDebugMode = false;
0374 
0375 // pT binning as in paragraph 3.2 of CMS-PAS-TRK-10-005 (https://cds.cern.ch/record/1279383/files/TRK-10-005-pas.pdf)
0376 // this is the default
0377 
0378 std::array<float, nPtBins_ + 1> mypT_bins = {{0.5,  0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6,  1.7,
0379                                               1.8,  1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9,  3.0,
0380                                               3.1,  3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.25, 4.5,
0381                                               4.75, 5.0, 5.5, 6.,  7.,  8.,  9.,  11., 14., 20.}};
0382 
0383 // inline function
0384 int check(const double a[], int n) {
0385   //while (--n > 0 && a[n] == a[0])    // exact match
0386   while (--n > 0 && (a[n] - a[0]) < 0.01)  // merged input files, protection agains numerical precision
0387     ;
0388   return n != 0;
0389 }
0390 
0391 // fill the list of files
0392 //*************************************************************
0393 void loadFileList(const char *inputFile, TString baseDir, TString legendName, int lineColor, int lineStyle)
0394 //*************************************************************
0395 {
0396   gErrorIgnoreLevel = kFatal;
0397   sourceList.push_back(new PVValidationVariables(inputFile, baseDir, legendName, lineColor, lineStyle));
0398 }
0399 
0400 //*************************************************************
0401 void FitPVResiduals(TString namesandlabels, bool stdres, bool do2DMaps, TString theDate, bool setAutoLimits)
0402 //*************************************************************
0403 {
0404   // only for fatal errors (useful in debugging)
0405   gErrorIgnoreLevel = kFatal;
0406 
0407   TH1::AddDirectory(kFALSE);
0408   bool fromLoader = false;
0409 
0410   TTimeStamp start_time;
0411   TStopwatch timer;
0412   timer.Start();
0413 
0414   if (!setAutoLimits) {
0415     std::cout << "FitPVResiduals::FitPVResiduals(): Overriding autolimits!" << std::endl;
0416     thePlotLimits->printAll();
0417   } else {
0418     std::cout << "FitPVResiduals::FitPVResiduals(): plot axis range will be automatically adjusted" << std::endl;
0419   }
0420 
0421   //TkAlStyle::set(INTERNAL);   // set publication status
0422 
0423   Int_t def_markers[9] = {kFullSquare,
0424                           kFullCircle,
0425                           kFullTriangleDown,
0426                           kOpenSquare,
0427                           kDot,
0428                           kOpenCircle,
0429                           kFullTriangleDown,
0430                           kFullTriangleUp,
0431                           kOpenTriangleDown};
0432   Int_t def_colors[9] = {kBlack, kRed, kBlue, kMagenta, kGreen, kCyan, kViolet, kOrange, kGreen + 2};
0433 
0434   Int_t markers[9];
0435   Int_t colors[9];
0436 
0437   setStyle();
0438 
0439   // check if the loader is empty
0440   if (sourceList.size() != 0) {
0441     fromLoader = true;
0442   }
0443 
0444   // if enters here, whatever is passed from command line is neglected
0445   if (fromLoader) {
0446     std::cout << "FitPVResiduals::FitPVResiduals(): file list specified from loader" << std::endl;
0447     std::cout << "======================================================" << std::endl;
0448     std::cout << "!!    arguments passed from CLI will be neglected   !!" << std::endl;
0449     std::cout << "======================================================" << std::endl;
0450     for (std::vector<PVValidationVariables *>::iterator it = sourceList.begin(); it != sourceList.end(); ++it) {
0451       std::cout << "name:  " << std::setw(20) << (*it)->getName() << " |file:  " << std::setw(15) << (*it)->getFile()
0452                 << " |color: " << std::setw(5) << (*it)->getLineColor() << " |style: " << std::setw(5)
0453                 << (*it)->getLineStyle() << " |marker:" << std::setw(5) << (*it)->getMarkerStyle() << std::endl;
0454     }
0455     std::cout << "======================================================" << std::endl;
0456   }
0457 
0458   Int_t theFileCount = 0;
0459   TList *FileList = new TList();
0460   TList *LabelList = new TList();
0461 
0462   if (!fromLoader) {
0463     namesandlabels.Remove(TString::kTrailing, ',');
0464     TObjArray *nameandlabelpairs = namesandlabels.Tokenize(",");
0465     for (Int_t i = 0; i < nameandlabelpairs->GetEntries(); ++i) {
0466       TObjArray *aFileLegPair = TString(nameandlabelpairs->At(i)->GetName()).Tokenize("=");
0467 
0468       if (aFileLegPair->GetEntries() == 2) {
0469         FileList->Add(TFile::Open(aFileLegPair->At(0)->GetName(), "READ"));  // 2
0470         LabelList->Add(aFileLegPair->At(1));
0471       } else {
0472         std::cout << "Please give file name and legend entry in the following form:\n"
0473                   << " filename1=legendentry1,filename2=legendentry2\n";
0474         exit(EXIT_FAILURE);
0475       }
0476     }
0477     theFileCount = FileList->GetSize();
0478   } else {
0479     for (std::vector<PVValidationVariables *>::iterator it = sourceList.begin(); it != sourceList.end(); ++it) {
0480       //FileList->Add((*it)->getFile()); // was extremely slow
0481       FileList->Add(TFile::Open((*it)->getFileName(), "READ"));
0482     }
0483     theFileCount = sourceList.size();
0484   }
0485 
0486   if (theFileCount == 0) {
0487     std::cout << "FitPVResiduals::FitPVResiduals(): empty input file list has been passed." << std::endl;
0488     exit(EXIT_FAILURE);
0489   }
0490 
0491   //  const Int_t nFiles_ = FileList->GetSize();
0492   const Int_t nFiles_ = theFileCount;
0493   TString LegLabels[10];
0494   TFile *fins[nFiles_];
0495 
0496   // already used in the global variables
0497   //const Int_t nBins_ =24;
0498 
0499   for (Int_t j = 0; j < nFiles_; j++) {
0500     // Retrieve files
0501     fins[j] = (TFile *)FileList->At(j);
0502 
0503     // Retrieve labels
0504     if (!fromLoader) {
0505       TObjString *legend = (TObjString *)LabelList->At(j);
0506       LegLabels[j] = legend->String();
0507       markers[j] = def_markers[j];
0508       colors[j] = def_colors[j];
0509     } else {
0510       LegLabels[j] = sourceList[j]->getName();
0511       markers[j] = sourceList[j]->getMarkerStyle();
0512       colors[j] = sourceList[j]->getLineColor();
0513     }
0514     LegLabels[j].ReplaceAll("_", " ");
0515     std::cout << "FitPVResiduals::FitPVResiduals(): label[" << j << "] " << LegLabels[j] << std::endl;
0516   }
0517 
0518   //
0519   // initialize all the histograms to be taken from file
0520   //
0521 
0522   // integrated residuals
0523   TH1F *dxyRefit[nFiles_];
0524   TH1F *dzRefit[nFiles_];
0525 
0526   TH1F *dxySigRefit[nFiles_];
0527   TH1F *dzSigRefit[nFiles_];
0528 
0529   // dca absolute residuals
0530   TH1F *dxyPhiResiduals[nFiles_][nBins_];
0531   TH1F *dxyEtaResiduals[nFiles_][nBins_];
0532 
0533   TH1F *dzPhiResiduals[nFiles_][nBins_];
0534   TH1F *dzEtaResiduals[nFiles_][nBins_];
0535 
0536   // dca transvers residuals
0537   TH1F *dxPhiResiduals[nFiles_][nBins_];
0538   TH1F *dxEtaResiduals[nFiles_][nBins_];
0539 
0540   TH1F *dyPhiResiduals[nFiles_][nBins_];
0541   TH1F *dyEtaResiduals[nFiles_][nBins_];
0542 
0543   // dca normalized residuals
0544   TH1F *dxyNormPhiResiduals[nFiles_][nBins_];
0545   TH1F *dxyNormEtaResiduals[nFiles_][nBins_];
0546 
0547   TH1F *dzNormPhiResiduals[nFiles_][nBins_];
0548   TH1F *dzNormEtaResiduals[nFiles_][nBins_];
0549 
0550   // double-differential residuals
0551   TH1F *dxyMapResiduals[nFiles_][nBins_][nBins_];
0552   TH1F *dzMapResiduals[nFiles_][nBins_][nBins_];
0553 
0554   TH1F *dxyNormMapResiduals[nFiles_][nBins_][nBins_];
0555   TH1F *dzNormMapResiduals[nFiles_][nBins_][nBins_];
0556 
0557   // double-differential residuals L1
0558   TH1F *dxyL1MapResiduals[nFiles_][nLadders_][nModZ_];
0559   TH1F *dzL1MapResiduals[nFiles_][nLadders_][nModZ_];
0560 
0561   TH1F *dxyL1NormMapResiduals[nFiles_][nLadders_][nModZ_];
0562   TH1F *dzL1NormMapResiduals[nFiles_][nLadders_][nModZ_];
0563 
0564   // dca residuals vs pT
0565   TH1F *dzNormPtResiduals[nFiles_][nPtBins_];
0566   TH1F *dxyNormPtResiduals[nFiles_][nPtBins_];
0567 
0568   TH1F *dzPtResiduals[nFiles_][nPtBins_];
0569   TH1F *dxyPtResiduals[nFiles_][nPtBins_];
0570 
0571   // dca residuals vs module number/ladder
0572   TH1F *dzNormLadderResiduals[nFiles_][nLadders_];
0573   TH1F *dxyNormLadderResiduals[nFiles_][nLadders_];
0574 
0575   TH1F *dzLadderResiduals[nFiles_][nLadders_];
0576   TH1F *dxyLadderResiduals[nFiles_][nLadders_];
0577 
0578   TH1F *dzNormModZResiduals[nFiles_][nModZ_];
0579   TH1F *dxyNormModZResiduals[nFiles_][nModZ_];
0580 
0581   TH1F *dzModZResiduals[nFiles_][nModZ_];
0582   TH1F *dxyModZResiduals[nFiles_][nModZ_];
0583 
0584   // for sanity checks
0585   TH1F *theEtaHistos[nFiles_];
0586   TH1F *thebinsHistos[nFiles_];
0587   TH1F *theLaddersHistos[nFiles_];
0588   TH1F *theModZHistos[nFiles_];
0589   TH1F *thePtInfoHistos[nFiles_];
0590 
0591   double theEtaMax_[nFiles_];
0592   double theNBINS[nFiles_];
0593   double theLadders[nFiles_];
0594   double theModZ[nFiles_];
0595 
0596   double thePtMax[nFiles_];
0597   double thePtMin[nFiles_];
0598   double thePTBINS[nFiles_];
0599 
0600   TTimeStamp initialization_done;
0601 
0602   if (isDebugMode) {
0603     timer.Stop();
0604     std::cout << "check point 1: " << timer.CpuTime() << " " << timer.RealTime() << std::endl;
0605     timer.Continue();
0606   }
0607 
0608   for (Int_t i = 0; i < nFiles_; i++) {
0609     fins[i]->cd("PVValidation/EventFeatures/");
0610 
0611     if (gDirectory->GetListOfKeys()->Contains("etaMax")) {
0612       gDirectory->GetObject("etaMax", theEtaHistos[i]);
0613       theEtaMax_[i] = theEtaHistos[i]->GetBinContent(1) / theEtaHistos[i]->GetEntries();
0614       std::cout << "File n. " << i << " has theEtaMax[" << i << "] = " << theEtaMax_[i] << std::endl;
0615     } else {
0616       theEtaMax_[i] = 2.5;
0617       std::cout << "File n. " << i << " getting the default pseudo-rapidity range: " << theEtaMax_[i] << std::endl;
0618     }
0619 
0620     if (gDirectory->GetListOfKeys()->Contains("nbins")) {
0621       gDirectory->GetObject("nbins", thebinsHistos[i]);
0622       theNBINS[i] = thebinsHistos[i]->GetBinContent(1) / thebinsHistos[i]->GetEntries();
0623       std::cout << "File n. " << i << " has theNBINS[" << i << "] = " << theNBINS[i] << std::endl;
0624     } else {
0625       theNBINS[i] = 48.;
0626       std::cout << "File n. " << i << " getting the default n. of bins: " << theNBINS[i] << std::endl;
0627     }
0628 
0629     if (gDirectory->GetListOfKeys()->Contains("nladders")) {
0630       gDirectory->GetObject("nladders", theLaddersHistos[i]);
0631       theLadders[i] = theLaddersHistos[i]->GetBinContent(1) / theLaddersHistos[i]->GetEntries();
0632       std::cout << "File n. " << i << " has theNLadders[" << i << "] = " << theLadders[i] << std::endl;
0633     } else {
0634       theLadders[i] = -1.;
0635       std::cout << "File n. " << i << " getting the default n. ladders: " << theLadders[i] << std::endl;
0636     }
0637 
0638     if (gDirectory->GetListOfKeys()->Contains("nModZ")) {
0639       gDirectory->GetObject("nModZ", theModZHistos[i]);
0640       theModZ[i] = theModZHistos[i]->GetBinContent(1) / theModZHistos[i]->GetEntries();
0641       std::cout << "File n. " << i << " has theNModZ[" << i << "] = " << theModZ[i] << std::endl;
0642     } else {
0643       theModZ[i] = -1.;
0644       std::cout << "File n. " << i << " getting the default n. modules along Z: " << theModZ[i] << std::endl;
0645     }
0646 
0647     if (gDirectory->GetListOfKeys()->Contains("pTinfo")) {
0648       gDirectory->GetObject("pTinfo", thePtInfoHistos[i]);
0649       thePTBINS[i] = thePtInfoHistos[i]->GetBinContent(1) * 3. / thePtInfoHistos[i]->GetEntries();
0650       ;
0651       thePtMin[i] = thePtInfoHistos[i]->GetBinContent(2) * 3. / thePtInfoHistos[i]->GetEntries();
0652       thePtMax[i] = thePtInfoHistos[i]->GetBinContent(3) * 3. / thePtInfoHistos[i]->GetEntries();
0653       std::cout << "File n. " << i << " has thePTBINS[" << i << "] = " << thePTBINS[i] << " pT min:  " << thePtMin[i]
0654                 << " pT max: " << thePtMax[i] << std::endl;
0655     } else {
0656       // if there is no histogram then set the extremes to 0, but the n. of bins should be still the default
0657       // this protects running against the old file format.
0658       thePTBINS[i] = nPtBins_;
0659       thePtMin[i] = 0.;
0660       thePtMax[i] = 0.;
0661       std::cout << "File n. " << i << " getting the default pT binning: ";
0662       for (const auto &bin : mypT_bins) {
0663         std::cout << bin << " ";
0664       }
0665       std::cout << std::endl;
0666     }
0667 
0668     // get the non-differential residuals plots
0669 
0670     fins[i]->cd("PVValidation/ProbeTrackFeatures");
0671 
0672     gDirectory->GetObject("h_probedxyRefitV", dxyRefit[i]);
0673     gDirectory->GetObject("h_probedzRefitV", dzRefit[i]);
0674 
0675     gDirectory->GetObject("h_probeRefitVSigXY", dxySigRefit[i]);
0676     gDirectory->GetObject("h_probeRefitVSigZ", dzSigRefit[i]);
0677 
0678     for (Int_t j = 0; j < theNBINS[i]; j++) {
0679       if (stdres) {
0680         // DCA absolute residuals
0681 
0682         fins[i]->cd("PVValidation/Abs_Transv_Phi_Residuals/");
0683         gDirectory->GetObject(Form("histo_dxy_phi_plot%i", j), dxyPhiResiduals[i][j]);
0684         gDirectory->GetObject(Form("histo_dx_phi_plot%i", j), dxPhiResiduals[i][j]);
0685         gDirectory->GetObject(Form("histo_dy_phi_plot%i", j), dyPhiResiduals[i][j]);
0686 
0687         fins[i]->cd("PVValidation/Abs_Transv_Eta_Residuals/");
0688         gDirectory->GetObject(Form("histo_dxy_eta_plot%i", j), dxyEtaResiduals[i][j]);
0689         gDirectory->GetObject(Form("histo_dx_eta_plot%i", j), dxEtaResiduals[i][j]);
0690         gDirectory->GetObject(Form("histo_dy_eta_plot%i", j), dyEtaResiduals[i][j]);
0691 
0692         dzPhiResiduals[i][j] = (TH1F *)fins[i]->Get(Form("PVValidation/Abs_Long_Phi_Residuals/histo_dz_phi_plot%i", j));
0693         dzEtaResiduals[i][j] = (TH1F *)fins[i]->Get(Form("PVValidation/Abs_Long_Eta_Residuals/histo_dz_eta_plot%i", j));
0694 
0695         // DCA normalized residuals
0696         dxyNormPhiResiduals[i][j] =
0697             (TH1F *)fins[i]->Get(Form("PVValidation/Norm_Transv_Phi_Residuals/histo_norm_dxy_phi_plot%i", j));
0698         dxyNormEtaResiduals[i][j] =
0699             (TH1F *)fins[i]->Get(Form("PVValidation/Norm_Transv_Eta_Residuals/histo_norm_dxy_eta_plot%i", j));
0700         dzNormPhiResiduals[i][j] =
0701             (TH1F *)fins[i]->Get(Form("PVValidation/Norm_Long_Phi_Residuals/histo_norm_dz_phi_plot%i", j));
0702         dzNormEtaResiduals[i][j] =
0703             (TH1F *)fins[i]->Get(Form("PVValidation/Norm_Long_Eta_Residuals/histo_norm_dz_eta_plot%i", j));
0704 
0705         // double differential residuals
0706 
0707         if (do2DMaps) {
0708           for (Int_t k = 0; k < theNBINS[i]; k++) {
0709             // absolute residuals
0710             fins[i]->cd("PVValidation/Abs_DoubleDiffResiduals/");
0711             gDirectory->GetObject(Form("histo_dxy_eta_plot%i_phi_plot%i", j, k), dxyMapResiduals[i][j][k]);
0712             gDirectory->GetObject(Form("histo_dz_eta_plot%i_phi_plot%i", j, k), dzMapResiduals[i][j][k]);
0713 
0714             // normalized residuals
0715             fins[i]->cd("PVValidation/Norm_DoubleDiffResiduals/");
0716             gDirectory->GetObject(Form("histo_norm_dxy_eta_plot%i_phi_plot%i", j, k), dxyNormMapResiduals[i][j][k]);
0717             gDirectory->GetObject(Form("histo_norm_dz_eta_plot%i_phi_plot%i", j, k), dzNormMapResiduals[i][j][k]);
0718           }
0719         }
0720       } else {
0721         // DCA absolute residuals
0722         dxyPhiResiduals[i][j] =
0723             (TH1F *)fins[i]->Get(Form("PVValidation/Abs_Transv_Phi_Residuals/histo_IP2D_phi_plot%i", j));
0724         dxyEtaResiduals[i][j] =
0725             (TH1F *)fins[i]->Get(Form("PVValidation/Abs_Transv_Eta_Residuals/histo_IP2D_eta_plot%i", j));
0726         dzPhiResiduals[i][j] =
0727             (TH1F *)fins[i]->Get(Form("PVValidation/Abs_Long_Phi_Residuals/histo_resz_phi_plot%i", j));
0728         dzEtaResiduals[i][j] =
0729             (TH1F *)fins[i]->Get(Form("PVValidation/Abs_Long_Eta_Residuals/histo_resz_eta_plot%i", j));
0730 
0731         // DCA normalized residuals
0732         dxyNormPhiResiduals[i][j] =
0733             (TH1F *)fins[i]->Get(Form("PVValidation/Norm_Transv_Phi_Residuals/histo_norm_IP2D_phi_plot%i", j));
0734         dxyNormEtaResiduals[i][j] =
0735             (TH1F *)fins[i]->Get(Form("PVValidation/Norm_Transv_Eta_Residuals/histo_norm_IP2D_eta_plot%i", j));
0736         dzNormPhiResiduals[i][j] =
0737             (TH1F *)fins[i]->Get(Form("PVValidation/Norm_Long_Phi_Residuals/histo_norm_resz_phi_plot%i", j));
0738         dzNormEtaResiduals[i][j] =
0739             (TH1F *)fins[i]->Get(Form("PVValidation/Norm_Long_Eta_Residuals/histo_norm_resz_eta_plot%i", j));
0740 
0741         // double differential residuals
0742         if (do2DMaps) {
0743           for (Int_t k = 0; k < theNBINS[i]; k++) {
0744             // absolute residuals
0745             fins[i]->cd("PVValidation/Abs_DoubleDiffResiduals");
0746             gDirectory->GetObject(Form("PVValidation/Abs_DoubleDiffResiduals/histo_dxy_eta_plot%i_phi_plot%i", j, k),
0747                                   dxyMapResiduals[i][j][k]);
0748             gDirectory->GetObject(Form("PVValidation/Abs_DoubleDiffResiduals/histo_dz_eta_plot%i_phi_plot%i", j, k),
0749                                   dzMapResiduals[i][j][k]);
0750 
0751             // normalized residuals
0752             fins[i]->cd("PVValidation/Norm_DoubleDiffResiduals");
0753             gDirectory->GetObject(
0754                 Form("PVValidation/Norm_DoubleDiffResiduals/histo_norm_dxy_eta_plot%i_phi_plot%i", j, k),
0755                 dxyNormMapResiduals[i][j][k]);
0756             gDirectory->GetObject(
0757                 Form("PVValidation/Norm_DoubleDiffResiduals/histo_norm_dz_eta_plot%i_phi_plot%i", j, k),
0758                 dzNormMapResiduals[i][j][k]);
0759           }
0760         }  // if do2DMaps
0761       }
0762     }
0763 
0764     // residuals vs pT
0765 
0766     for (Int_t l = 0; l < thePTBINS[i] - 1; l++) {
0767       dxyPtResiduals[i][l] = (TH1F *)fins[i]->Get(Form("PVValidation/Abs_Transv_pT_Residuals/histo_dxy_pT_plot%i", l));
0768       dzPtResiduals[i][l] = (TH1F *)fins[i]->Get(Form("PVValidation/Abs_Long_pT_Residuals/histo_dz_pT_plot%i", l));
0769 
0770       dxyNormPtResiduals[i][l] =
0771           (TH1F *)fins[i]->Get(Form("PVValidation/Norm_Transv_pT_Residuals/histo_norm_dxy_pT_plot%i", l));
0772       dzNormPtResiduals[i][l] =
0773           (TH1F *)fins[i]->Get(Form("PVValidation/Norm_Long_pT_Residuals/histo_norm_dz_pT_plot%i", l));
0774     }
0775 
0776     // residuals vs module number / ladder
0777 
0778     if (theLadders[i] > 0 && theModZ[i] > 0) {
0779       for (Int_t iLadder = 0; iLadder < theLadders[i]; iLadder++) {
0780         dzNormLadderResiduals[i][iLadder] =
0781             (TH1F *)fins[i]->Get(Form("PVValidation/Norm_Long_ladder_Residuals/histo_norm_dz_ladder_plot%i", iLadder));
0782         dxyNormLadderResiduals[i][iLadder] = (TH1F *)fins[i]->Get(
0783             Form("PVValidation/Norm_Transv_ladder_Residuals/histo_norm_dxy_ladder_plot%i", iLadder));
0784 
0785         dzLadderResiduals[i][iLadder] =
0786             (TH1F *)fins[i]->Get(Form("PVValidation/Abs_Long_ladder_Residuals/histo_dz_ladder_plot%i", iLadder));
0787         dxyLadderResiduals[i][iLadder] = (TH1F *)fins[i]->Get(
0788             Form("PVValidation/Abs_Transv_ladderNoOverlap_Residuals/histo_dxy_ladder_plot%i", iLadder));
0789 
0790         if (do2DMaps) {
0791           for (Int_t iMod = 0; iMod < theModZ[i]; iMod++) {
0792             dxyL1MapResiduals[i][iLadder][iMod] =
0793                 (TH1F *)fins[i]->Get(Form("PVValidation/Abs_L1Residuals/histo_dxy_ladder%i_module%i", iLadder, iMod));
0794             dzL1MapResiduals[i][iLadder][iMod] =
0795                 (TH1F *)fins[i]->Get(Form("PVValidation/Abs_L1Residuals/histo_dz_ladder%i_module%i", iLadder, iMod));
0796 
0797             dxyL1NormMapResiduals[i][iLadder][iMod] = (TH1F *)fins[i]->Get(
0798                 Form("PVValidation/Norm_L1Residuals/histo_norm_dxy_ladder%i_module%i", iLadder, iMod));
0799             dzL1NormMapResiduals[i][iLadder][iMod] = (TH1F *)fins[i]->Get(
0800                 Form("PVValidation/Norm_L1Residuals/histo_norm_dz_ladder%i_module%i", iLadder, iMod));
0801           }
0802         }
0803       }
0804     }
0805 
0806     if (theModZ[i] > 0) {
0807       for (Int_t iMod = 0; iMod < theModZ[i]; iMod++) {
0808         dzNormModZResiduals[i][iMod] =
0809             (TH1F *)fins[i]->Get(Form("PVValidation/Norm_Long_modZ_Residuals/histo_norm_dz_modZ_plot%i", iMod));
0810         dxyNormModZResiduals[i][iMod] =
0811             (TH1F *)fins[i]->Get(Form("PVValidation/Norm_Transv_modZ_Residuals/histo_norm_dxy_modZ_plot%i", iMod));
0812 
0813         dzModZResiduals[i][iMod] =
0814             (TH1F *)fins[i]->Get(Form("PVValidation/Abs_Long_modZ_Residuals/histo_dz_modZ_plot%i", iMod));
0815         dxyModZResiduals[i][iMod] =
0816             (TH1F *)fins[i]->Get(Form("PVValidation/Abs_Transv_modZ_Residuals/histo_dxy_modZ_plot%i", iMod));
0817       }
0818     }
0819 
0820     // close the files after retrieving them
0821     fins[i]->Close();
0822   }
0823 
0824   TTimeStamp caching_done;
0825 
0826   if (isDebugMode) {
0827     timer.Stop();
0828     std::cout << "check point 2: " << timer.CpuTime() << " " << timer.RealTime() << std::endl;
0829     timer.Continue();
0830   }
0831 
0832   // checks if all pseudo-rapidity ranges coincide
0833   // if not, exits
0834   if (check(theEtaMax_, nFiles_)) {
0835     std::cout << "======================================================" << std::endl;
0836     std::cout << "FitPVResiduals::FitPVResiduals(): the eta range is different" << std::endl;
0837     std::cout << "exiting..." << std::endl;
0838     exit(EXIT_FAILURE);
0839   } else {
0840     etaRange = theEtaMax_[0];
0841     std::cout << "======================================================" << std::endl;
0842     std::cout << "FitPVResiduals::FitPVResiduals(): the eta range is [" << -etaRange << " ; " << etaRange << "]"
0843               << std::endl;
0844     std::cout << "======================================================" << std::endl;
0845   }
0846 
0847   // checks if all nbins ranges coincide
0848   // if not, exits
0849   if (check(theNBINS, nFiles_)) {
0850     std::cout << "======================================================" << std::endl;
0851     std::cout << "FitPVResiduals::FitPVResiduals(): the number of bins is different" << std::endl;
0852     std::cout << "exiting..." << std::endl;
0853     exit(EXIT_FAILURE);
0854   } else {
0855     nBins_ = theNBINS[0];
0856 
0857     // adjust also the limit for the fit
0858     _boundSx = (nBins_ / 4.) - 0.5;
0859     _boundDx = 3 * (nBins_ / 4.) - 0.5;
0860     _boundMax = nBins_ - 0.5;
0861 
0862     std::cout << "======================================================" << std::endl;
0863     std::cout << "FitPVResiduals::FitPVResiduals(): the number of bins is: " << nBins_ << std::endl;
0864     std::cout << "======================================================" << std::endl;
0865   }
0866 
0867   // checks if the geometries are consistent to produce the ladder plots
0868   if (check(theLadders, nFiles_)) {
0869     std::cout << "======================================================" << std::endl;
0870     std::cout << "FitPVResiduals::FitPVResiduals(): the number of ladders is different" << std::endl;
0871     std::cout << "won't do the ladder analysis..." << std::endl;
0872     std::cout << "======================================================" << std::endl;
0873     nLadders_ = -1;
0874   } else {
0875     nLadders_ = theLadders[0];
0876     std::cout << "======================================================" << std::endl;
0877     std::cout << "FitPVResiduals::FitPVResiduals(): the number of ladders is: " << nLadders_ << std::endl;
0878     std::cout << "======================================================" << std::endl;
0879   }
0880 
0881   // checks if the geometries are consistent to produce the moduleZ plots
0882   if (check(theModZ, nFiles_)) {
0883     std::cout << "======================================================" << std::endl;
0884     std::cout << "FitPVResiduals::FitPVResiduals(): the number of modules in Z is different" << std::endl;
0885     std::cout << "won't do the ladder analysis..." << std::endl;
0886     std::cout << "======================================================" << std::endl;
0887     nModZ_ = -1;
0888   } else {
0889     nModZ_ = theModZ[0];
0890     std::cout << "======================================================" << std::endl;
0891     std::cout << "FitPVResiduals::FitPVResiduals(): the number of modules in Z is: " << nModZ_ << std::endl;
0892     std::cout << "======================================================" << std::endl;
0893   }
0894 
0895   // checks if pT boundaries are consistent to produce the pT-binned plots
0896   if (check(thePtMax, nFiles_) || check(thePtMin, nFiles_)) {
0897     std::cout << "======================================================" << std::endl;
0898     std::cout << "FitPVResiduals::FitPVResiduals(): the pT binning is different" << std::endl;
0899     std::cout << "won't do the pT analysis..." << std::endl;
0900     std::cout << "======================================================" << std::endl;
0901     minPt_ = -1.;
0902   } else {
0903     if (thePtMin[0] != 0.) {
0904       minPt_ = thePtMin[0];
0905       maxPt_ = thePtMax[0];
0906       mypT_bins = PVValHelper::makeLogBins<float, nPtBins_>(thePtMin[0], thePtMax[0]);
0907       std::cout << "======================================================" << std::endl;
0908       std::cout << "FitPVResiduals::FitPVResiduals(): log bins [" << thePtMin[0] << "," << thePtMax[0] << "]"
0909                 << std::endl;
0910       std::cout << "======================================================" << std::endl;
0911     } else {
0912       std::cout << "======================================================" << std::endl;
0913       std::cout << "FitPVResiduals::FitPVResiduals(): using default bins ";
0914       for (const auto &bin : mypT_bins) {
0915         std::cout << bin << " ";
0916       }
0917       std::cout << std::endl;
0918       std::cout << "======================================================" << std::endl;
0919     }
0920   }
0921 
0922   // check now that all the files have events
0923   bool areAllFilesFull = true;
0924   for (Int_t i = 0; i < nFiles_; i++) {
0925     if (dxyRefit[i]->GetEntries() == 0.) {
0926       areAllFilesFull = false;
0927       break;
0928     }
0929   }
0930 
0931   if (!areAllFilesFull) {
0932     std::cout << "======================================================" << std::endl;
0933     std::cout << "FitPVResiduals::FitPVResiduals(): not all the files have events" << std::endl;
0934     std::cout << "exiting (...to prevent a segmentation fault)" << std::endl;
0935     exit(EXIT_FAILURE);
0936   }
0937 
0938   Double_t highedge = nBins_ - 0.5;
0939   Double_t lowedge = -0.5;
0940 
0941   // DCA absolute
0942 
0943   TH1F *dxyPhiMeanTrend[nFiles_];
0944   TH1F *dxyPhiWidthTrend[nFiles_];
0945 
0946   TH1F *dxPhiMeanTrend[nFiles_];
0947   TH1F *dxPhiWidthTrend[nFiles_];
0948 
0949   TH1F *dyPhiMeanTrend[nFiles_];
0950   TH1F *dyPhiWidthTrend[nFiles_];
0951 
0952   TH1F *dzPhiMeanTrend[nFiles_];
0953   TH1F *dzPhiWidthTrend[nFiles_];
0954 
0955   TH1F *dxyEtaMeanTrend[nFiles_];
0956   TH1F *dxyEtaWidthTrend[nFiles_];
0957 
0958   TH1F *dxEtaMeanTrend[nFiles_];
0959   TH1F *dxEtaWidthTrend[nFiles_];
0960 
0961   TH1F *dyEtaMeanTrend[nFiles_];
0962   TH1F *dyEtaWidthTrend[nFiles_];
0963 
0964   TH1F *dzEtaMeanTrend[nFiles_];
0965   TH1F *dzEtaWidthTrend[nFiles_];
0966 
0967   TH1F *dxyPtMeanTrend[nFiles_];
0968   TH1F *dxyPtWidthTrend[nFiles_];
0969 
0970   TH1F *dzPtMeanTrend[nFiles_];
0971   TH1F *dzPtWidthTrend[nFiles_];
0972 
0973   // vs ladder and module number
0974 
0975   TH1F *dxyLadderMeanTrend[nFiles_];
0976   TH1F *dxyLadderWidthTrend[nFiles_];
0977   TH1F *dzLadderMeanTrend[nFiles_];
0978   TH1F *dzLadderWidthTrend[nFiles_];
0979 
0980   TH1F *dxyModZMeanTrend[nFiles_];
0981   TH1F *dxyModZWidthTrend[nFiles_];
0982   TH1F *dzModZMeanTrend[nFiles_];
0983   TH1F *dzModZWidthTrend[nFiles_];
0984 
0985   // DCA normalized
0986 
0987   TH1F *dxyNormPhiMeanTrend[nFiles_];
0988   TH1F *dxyNormPhiWidthTrend[nFiles_];
0989   TH1F *dzNormPhiMeanTrend[nFiles_];
0990   TH1F *dzNormPhiWidthTrend[nFiles_];
0991 
0992   TH1F *dxyNormEtaMeanTrend[nFiles_];
0993   TH1F *dxyNormEtaWidthTrend[nFiles_];
0994   TH1F *dzNormEtaMeanTrend[nFiles_];
0995   TH1F *dzNormEtaWidthTrend[nFiles_];
0996 
0997   TH1F *dxyNormPtMeanTrend[nFiles_];
0998   TH1F *dxyNormPtWidthTrend[nFiles_];
0999   TH1F *dzNormPtMeanTrend[nFiles_];
1000   TH1F *dzNormPtWidthTrend[nFiles_];
1001 
1002   TH1F *dxyNormLadderMeanTrend[nFiles_];
1003   TH1F *dxyNormLadderWidthTrend[nFiles_];
1004   TH1F *dzNormLadderMeanTrend[nFiles_];
1005   TH1F *dzNormLadderWidthTrend[nFiles_];
1006 
1007   TH1F *dxyNormModZMeanTrend[nFiles_];
1008   TH1F *dxyNormModZWidthTrend[nFiles_];
1009   TH1F *dzNormModZMeanTrend[nFiles_];
1010   TH1F *dzNormModZWidthTrend[nFiles_];
1011 
1012   // 2D maps
1013 
1014   // bias
1015   TH2F *dxyMeanMap[nFiles_];
1016   TH2F *dzMeanMap[nFiles_];
1017   TH2F *dxyNormMeanMap[nFiles_];
1018   TH2F *dzNormMeanMap[nFiles_];
1019 
1020   // width
1021   TH2F *dxyWidthMap[nFiles_];
1022   TH2F *dzWidthMap[nFiles_];
1023   TH2F *dxyNormWidthMap[nFiles_];
1024   TH2F *dzNormWidthMap[nFiles_];
1025 
1026   // trimmed maps
1027 
1028   // bias
1029   TH2F *t_dxyMeanMap[nFiles_];
1030   TH2F *t_dzMeanMap[nFiles_];
1031   TH2F *t_dxyNormMeanMap[nFiles_];
1032   TH2F *t_dzNormMeanMap[nFiles_];
1033 
1034   // width
1035   TH2F *t_dxyWidthMap[nFiles_];
1036   TH2F *t_dzWidthMap[nFiles_];
1037   TH2F *t_dxyNormWidthMap[nFiles_];
1038   TH2F *t_dzNormWidthMap[nFiles_];
1039 
1040   // 2D L1 maps
1041 
1042   // bias
1043   TH2F *dxyMeanL1Map[nFiles_];
1044   TH2F *dzMeanL1Map[nFiles_];
1045   TH2F *dxyNormMeanL1Map[nFiles_];
1046   TH2F *dzNormMeanL1Map[nFiles_];
1047 
1048   // width
1049   TH2F *dxyWidthL1Map[nFiles_];
1050   TH2F *dzWidthL1Map[nFiles_];
1051   TH2F *dxyNormWidthL1Map[nFiles_];
1052   TH2F *dzNormWidthL1Map[nFiles_];
1053 
1054   // trimmed maps
1055 
1056   // bias
1057   TH2F *t_dxyMeanL1Map[nFiles_];
1058   TH2F *t_dzMeanL1Map[nFiles_];
1059   TH2F *t_dxyNormMeanL1Map[nFiles_];
1060   TH2F *t_dzNormMeanL1Map[nFiles_];
1061 
1062   // width
1063   TH2F *t_dxyWidthL1Map[nFiles_];
1064   TH2F *t_dzWidthL1Map[nFiles_];
1065   TH2F *t_dxyNormWidthL1Map[nFiles_];
1066   TH2F *t_dzNormWidthL1Map[nFiles_];
1067 
1068   for (Int_t i = 0; i < nFiles_; i++) {
1069     // DCA trend plots
1070 
1071     dxyPhiMeanTrend[i] = new TH1F(Form("means_dxy_phi_%i", i),
1072                                   "#LT d_{xy} #GT vs #phi sector;track #phi [rad];#LT d_{xy} #GT [#mum]",
1073                                   nBins_,
1074                                   lowedge,
1075                                   highedge);
1076     dxyPhiWidthTrend[i] = new TH1F(Form("widths_dxy_phi_%i", i),
1077                                    "#sigma(d_{xy}) vs #phi sector;track #phi [rad];#sigma(d_{xy}) [#mum]",
1078                                    nBins_,
1079                                    lowedge,
1080                                    highedge);
1081 
1082     dxPhiMeanTrend[i] = new TH1F(Form("means_dx_phi_%i", i),
1083                                  "#LT d_{x} #GT vs #phi sector;track #phi [rad];#LT d_{x} #GT [#mum]",
1084                                  nBins_,
1085                                  lowedge,
1086                                  highedge);
1087     dxPhiWidthTrend[i] = new TH1F(Form("widths_dx_phi_%i", i),
1088                                   "#sigma(d_{x}) vs #phi sector;track #phi [rad];#sigma(d_{x}) [#mum]",
1089                                   nBins_,
1090                                   lowedge,
1091                                   highedge);
1092 
1093     dyPhiMeanTrend[i] = new TH1F(Form("means_dy_phi_%i", i),
1094                                  "#LT d_{y} #GT vs #phi sector;track #phi [rad];#LT d_{y} #GT [#mum]",
1095                                  nBins_,
1096                                  lowedge,
1097                                  highedge);
1098     dyPhiWidthTrend[i] = new TH1F(Form("widths_dy_phi_%i", i),
1099                                   "#sigma(d_{y}) vs #phi sector;track #phi [rad];#sigma(d_{y}) [#mum]",
1100                                   nBins_,
1101                                   lowedge,
1102                                   highedge);
1103 
1104     dzPhiMeanTrend[i] = new TH1F(Form("means_dz_phi_%i", i),
1105                                  "#LT d_{z} #GT vs #phi sector;track #phi [rad];#LT d_{z} #GT [#mum]",
1106                                  nBins_,
1107                                  lowedge,
1108                                  highedge);
1109     dzPhiWidthTrend[i] = new TH1F(Form("widths_dz_phi_%i", i),
1110                                   "#sigma(d_{z}) vs #phi sector;track #phi [rad];#sigma(d_{z}) [#mum]",
1111                                   nBins_,
1112                                   lowedge,
1113                                   highedge);
1114 
1115     dxyEtaMeanTrend[i] = new TH1F(Form("means_dxy_eta_%i", i),
1116                                   "#LT d_{xy} #GT vs #eta sector;track #eta;#LT d_{xy} #GT [#mum]",
1117                                   nBins_,
1118                                   lowedge,
1119                                   highedge);
1120     dxyEtaWidthTrend[i] = new TH1F(Form("widths_dxy_eta_%i", i),
1121                                    "#sigma(d_{xy}) vs #eta sector;track #eta;#sigma(d_{xy}) [#mum]",
1122                                    nBins_,
1123                                    lowedge,
1124                                    highedge);
1125 
1126     dxEtaMeanTrend[i] = new TH1F(Form("means_dx_eta_%i", i),
1127                                  "#LT d_{x} #GT vs #eta sector;track #eta;#LT d_{x} #GT [#mum]",
1128                                  nBins_,
1129                                  lowedge,
1130                                  highedge);
1131     dxEtaWidthTrend[i] = new TH1F(Form("widths_dx_eta_%i", i),
1132                                   "#sigma(d_{x}) vs #eta sector;track #eta;#sigma(d_{x}) [#mum]",
1133                                   nBins_,
1134                                   lowedge,
1135                                   highedge);
1136 
1137     dyEtaMeanTrend[i] = new TH1F(Form("means_dy_eta_%i", i),
1138                                  "#LT d_{y} #GT vs #eta sector;track #eta;#LT d_{y} #GT [#mum]",
1139                                  nBins_,
1140                                  lowedge,
1141                                  highedge);
1142     dyEtaWidthTrend[i] = new TH1F(Form("widths_dy_eta_%i", i),
1143                                   "#sigma(d_{y}) vs #eta sector;track #eta;#sigma(d_{y}) [#mum]",
1144                                   nBins_,
1145                                   lowedge,
1146                                   highedge);
1147 
1148     dzEtaMeanTrend[i] = new TH1F(Form("means_dz_eta_%i", i),
1149                                  "#LT d_{z} #GT vs #eta sector;track #eta;#LT d_{z} #GT [#mum]",
1150                                  nBins_,
1151                                  lowedge,
1152                                  highedge);
1153     dzEtaWidthTrend[i] = new TH1F(Form("widths_dz_eta_%i", i),
1154                                   "#sigma(d_{xy}) vs #eta sector;track #eta;#sigma(d_{z}) [#mum]",
1155                                   nBins_,
1156                                   lowedge,
1157                                   highedge);
1158 
1159     if (minPt_ > 0.) {
1160       dxyPtMeanTrend[i] = new TH1F(Form("means_dxy_pT_%i", i),
1161                                    "#LT d_{xy} #GT vs p_{T} sector;track p_{T} [GeV];#LT d_{xy} #GT [#mum]",
1162                                    mypT_bins.size() - 1,
1163                                    mypT_bins.data());
1164       dxyPtWidthTrend[i] = new TH1F(Form("widths_dxy_pT_%i", i),
1165                                     "#sigma(d_{xy}) vs p_{T} sector;track p_{T} [GeV];#sigma(d_{xy}) [#mum]",
1166                                     mypT_bins.size() - 1,
1167                                     mypT_bins.data());
1168       dzPtMeanTrend[i] = new TH1F(Form("means_dz_pT_%i", i),
1169                                   "#LT d_{z} #GT vs p_{T} sector;track p_{T} [GeV];#LT d_{z} #GT [#mum]",
1170                                   mypT_bins.size() - 1,
1171                                   mypT_bins.data());
1172       dzPtWidthTrend[i] = new TH1F(Form("widths_dz_pT_%i", i),
1173                                    "#sigma(d_{z}) vs p_{T} sector;track p_{T} [GeV];#sigma(d_{z}) [#mum]",
1174                                    mypT_bins.size() - 1,
1175                                    mypT_bins.data());
1176     }
1177 
1178     if (nModZ_ > 0) {
1179       dxyModZMeanTrend[i] = new TH1F(Form("means_dxy_modZ_%i", i),
1180                                      "#LT d_{xy} #GT vs Layer 1 module number;module number;#LT d_{xy} #GT [#mum]",
1181                                      theModZ[i],
1182                                      0.,
1183                                      theModZ[i]);
1184       dxyModZWidthTrend[i] = new TH1F(Form("widths_dxy_modZ_%i", i),
1185                                       "#sigma(d_{xy}) vs Layer 1 module number;module number;#sigma(d_{xy}) [#mum]",
1186                                       theModZ[i],
1187                                       0.,
1188                                       theModZ[i]);
1189       dzModZMeanTrend[i] = new TH1F(Form("means_dz_modZ_%i", i),
1190                                     "#LT d_{z} #GT vs Layer 1 module number;module number;#LT d_{z} #GT [#mum]",
1191                                     theModZ[i],
1192                                     0.,
1193                                     theModZ[i]);
1194       dzModZWidthTrend[i] = new TH1F(Form("widths_dz_modZ_%i", i),
1195                                      "#sigma(d_{z}) vs Layer 1 module number;module number;#sigma(d_{z}) [#mum]",
1196                                      theModZ[i],
1197                                      0.,
1198                                      theModZ[i]);
1199     }
1200 
1201     if (nLadders_ > 0) {
1202       dxyLadderMeanTrend[i] = new TH1F(Form("means_dxy_ladder_%i", i),
1203                                        "#LT d_{xy} #GT vs Layer 1 ladder;ladder number;#LT d_{xy} #GT [#mum]",
1204                                        theLadders[i],
1205                                        0.,
1206                                        theLadders[i]);
1207       dxyLadderWidthTrend[i] = new TH1F(Form("widths_dxy_ladder_%i", i),
1208                                         "#sigma(d_{xy}) vs Layer 1 ladder;ladder number;#sigma(d_{xy}) [#mum]",
1209                                         theLadders[i],
1210                                         0.,
1211                                         theLadders[i]);
1212       dzLadderMeanTrend[i] = new TH1F(Form("means_dz_ladder_%i", i),
1213                                       "#LT d_{z} #GT vs Layer 1 ladder;ladder number;#LT d_{z} #GT [#mum]",
1214                                       theLadders[i],
1215                                       0.,
1216                                       theLadders[i]);
1217       dzLadderWidthTrend[i] = new TH1F(Form("widths_dz_ladder_%i", i),
1218                                        "#sigma(d_{z}) vs Layer 1 ladder;ladder number;#sigma(d_{z}) [#mum]",
1219                                        theLadders[i],
1220                                        0.,
1221                                        theLadders[i]);
1222     }
1223 
1224     // DCA normalized trend plots
1225 
1226     dxyNormPhiMeanTrend[i] =
1227         new TH1F(Form("means_dxyNorm_phi_%i", i),
1228                  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;track #phi [rad];#LT d_{xy}/#sigma_{d_{xy}} #GT",
1229                  nBins_,
1230                  lowedge,
1231                  highedge);
1232     dxyNormPhiWidthTrend[i] =
1233         new TH1F(Form("widths_dxyNorm_phi_%i", i),
1234                  "#sigma(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;track #phi [rad];#sigma(d_{xy}/#sigma_{d_{xy}})",
1235                  nBins_,
1236                  lowedge,
1237                  highedge);
1238     dzNormPhiMeanTrend[i] =
1239         new TH1F(Form("means_dzNorm_phi_%i", i),
1240                  "#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;track #phi [rad];#LT d_{z}/#sigma_{d_{z}} #GT",
1241                  nBins_,
1242                  lowedge,
1243                  highedge);
1244     dzNormPhiWidthTrend[i] =
1245         new TH1F(Form("widths_dzNorm_phi_%i", i),
1246                  "#sigma(d_{z}/#sigma_{d_{z}}) vs #phi sector;track #phi [rad];#sigma(d_{z}/#sigma_{d_{z}})",
1247                  nBins_,
1248                  lowedge,
1249                  highedge);
1250 
1251     dxyNormEtaMeanTrend[i] =
1252         new TH1F(Form("means_dxyNorm_eta_%i", i),
1253                  "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;track #eta;#LT d_{xy}/#sigma_{d_{xy}} #GT",
1254                  nBins_,
1255                  lowedge,
1256                  highedge);
1257     dxyNormEtaWidthTrend[i] =
1258         new TH1F(Form("widths_dxyNorm_eta_%i", i),
1259                  "#sigma(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;track #eta;#sigma(d_{xy}/#sigma_{d_{xy}})",
1260                  nBins_,
1261                  lowedge,
1262                  highedge);
1263     dzNormEtaMeanTrend[i] =
1264         new TH1F(Form("means_dzNorm_eta_%i", i),
1265                  "#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;track #eta;#LT d_{z}/#sigma_{d_{z}} #GT",
1266                  nBins_,
1267                  lowedge,
1268                  highedge);
1269     dzNormEtaWidthTrend[i] =
1270         new TH1F(Form("widths_dzNorm_eta_%i", i),
1271                  "#sigma(d_{z}/#sigma_{d_{z}}) vs #eta sector;track #eta;#sigma(d_{z}/#sigma_{d_{z}})",
1272                  nBins_,
1273                  lowedge,
1274                  highedge);
1275 
1276     if (minPt_ > 0.) {
1277       dxyNormPtMeanTrend[i] =
1278           new TH1F(Form("means_dxyNorm_pT_%i", i),
1279                    "#LT d_{xy}/#sigma_{d_{xy}} #GT vs p_{T} sector;track p_{T} [GeV];#LT d_{xy}/#sigma_{d_{xy}} #GT",
1280                    mypT_bins.size() - 1,
1281                    mypT_bins.data());
1282       dxyNormPtWidthTrend[i] =
1283           new TH1F(Form("widths_dxyNorm_pT_%i", i),
1284                    "#sigma(d_{xy}/#sigma_{d_{xy}}) vs p_{T} sector;track p_{T} [GeV];#sigma(d_{xy}/#sigma_{d_{xy}})",
1285                    mypT_bins.size() - 1,
1286                    mypT_bins.data());
1287       dzNormPtMeanTrend[i] =
1288           new TH1F(Form("means_dzNorm_pT_%i", i),
1289                    "#LT d_{z}/#sigma_{d_{z}} #GT vs p_{T} sector;track p_{T} [GeV];#LT d_{z}/#sigma_{d_{z}} #GT",
1290                    mypT_bins.size() - 1,
1291                    mypT_bins.data());
1292       dzNormPtWidthTrend[i] =
1293           new TH1F(Form("widths_dzNorm_pT_%i", i),
1294                    "#sigma(d_{z}/#sigma_{d_{z}}) vs p_{T} sector;track p_{T} [GeV];#sigma(d_{z}/#sigma_{d_{z}})",
1295                    mypT_bins.size() - 1,
1296                    mypT_bins.data());
1297     }
1298 
1299     if (nLadders_ > 0) {
1300       dxyNormLadderMeanTrend[i] =
1301           new TH1F(Form("means_dxyNorm_ladder_%i", i),
1302                    "#LT d_{xy}/#sigma_{d_{xy}} #GT vs Layer 1 ladder;ladder number;#LT d_{xy}/#sigma_{d_{xy}} #GT",
1303                    theLadders[i],
1304                    0.,
1305                    theLadders[i]);
1306       dxyNormLadderWidthTrend[i] =
1307           new TH1F(Form("widths_dxyNorm_ladder_%i", i),
1308                    "#sigma(d_{xy}/#sigma_{d_{xy}}) vs Layer 1 ladder;ladder number;#sigma(d_{xy}/#sigma_{d_{xy}})",
1309                    theLadders[i],
1310                    0.,
1311                    theLadders[i]);
1312       dzNormLadderMeanTrend[i] =
1313           new TH1F(Form("means_dzNorm_ladder_%i", i),
1314                    "#LT d_{z}/#sigma_{d_{z}} #GT vs Layer 1 ladder;ladder number;#LT d_{z}/#sigma_{d_{z}} #GT",
1315                    theLadders[i],
1316                    0.,
1317                    theLadders[i]);
1318       dzNormLadderWidthTrend[i] =
1319           new TH1F(Form("widths_dzNorm_ladder_%i", i),
1320                    "#sigma(d_{z}/#sigma_{d_{z}}) vs Layer 1 ladder;ladder number;#sigma(d_{z}/#sigma_{d_{z}})",
1321                    theLadders[i],
1322                    0.,
1323                    theLadders[i]);
1324     }
1325 
1326     if (nModZ_ > 0) {
1327       dxyNormModZMeanTrend[i] = new TH1F(
1328           Form("means_dxyNorm_modZ_%i", i),
1329           "#LT d_{xy}/#sigma_{d_{xy}} #GT vs Layer 1 module number;module number;#LT d_{xy}/#sigma_{d_{xy}} #GT",
1330           theModZ[i],
1331           0.,
1332           theModZ[i]);
1333       dxyNormModZWidthTrend[i] = new TH1F(
1334           Form("widths_dxyNorm_modZ_%i", i),
1335           "#sigma(d_{xy}/#sigma_{d_{xy}}) vs Layer 1 module number;module number;#sigma(d_{xy}/#sigma_{d_{xy}})",
1336           theModZ[i],
1337           0.,
1338           theModZ[i]);
1339       dzNormModZMeanTrend[i] =
1340           new TH1F(Form("means_dzNorm_modZ_%i", i),
1341                    "#LT d_{z}/#sigma_{d_{z}} #GT vs Layer 1 module number;module number;#LT d_{z}/#sigma_{d_{z}} #GT",
1342                    theModZ[i],
1343                    0.,
1344                    theModZ[i]);
1345       dzNormModZWidthTrend[i] =
1346           new TH1F(Form("widths_dzNorm_modZ_%i", i),
1347                    "#sigma(d_{z}/#sigma_{d_{z}}) vs Layer 1 module number;module number;#sigma(d_{z}/#sigma_{d_{z}})",
1348                    theModZ[i],
1349                    0.,
1350                    theModZ[i]);
1351     }
1352 
1353     // 2D maps
1354     dxyMeanMap[i] = new TH2F(Form("means_dxy_map_%i", i),
1355                              "#LT d_{xy} #GT map;track #eta;track #phi [rad];#LT d_{xy} #GT [#mum]",
1356                              nBins_,
1357                              lowedge,
1358                              highedge,
1359                              nBins_,
1360                              lowedge,
1361                              highedge);
1362     dzMeanMap[i] = new TH2F(Form("means_dz_map_%i", i),
1363                             "#LT d_{z} #GT map;track #eta;track #phi [rad];#LT d_{z} #GT [#mum]",
1364                             nBins_,
1365                             lowedge,
1366                             highedge,
1367                             nBins_,
1368                             lowedge,
1369                             highedge);
1370     dxyNormMeanMap[i] =
1371         new TH2F(Form("norm_means_dxy_map_%i", i),
1372                  "#LT d_{xy}/#sigma_{d_{xy}} #GT map;track #eta;track #phi [rad];#LT d_{xy}/#sigma_{d_{xy}} #GT",
1373                  nBins_,
1374                  lowedge,
1375                  highedge,
1376                  nBins_,
1377                  lowedge,
1378                  highedge);
1379     dzNormMeanMap[i] =
1380         new TH2F(Form("norm_means_dz_map_%i", i),
1381                  "#LT d_{z}/#sigma_{d_{z}} #GT map;track #eta;track #phi[rad];#LT d_{xy}/#sigma_{d_{z}} #GT",
1382                  nBins_,
1383                  lowedge,
1384                  highedge,
1385                  nBins_,
1386                  lowedge,
1387                  highedge);
1388 
1389     dxyWidthMap[i] = new TH2F(Form("widths_dxy_map_%i", i),
1390                               "#sigma_{d_{xy}} map;track #eta;track #phi [rad];#sigma(d_{xy}) [#mum]",
1391                               nBins_,
1392                               lowedge,
1393                               highedge,
1394                               nBins_,
1395                               lowedge,
1396                               highedge);
1397     dzWidthMap[i] = new TH2F(Form("widths_dz_map_%i", i),
1398                              "#sigma_{d_{z}} map;track #eta;track #phi [rad];#sigma(d_{z}) [#mum]",
1399                              nBins_,
1400                              lowedge,
1401                              highedge,
1402                              nBins_,
1403                              lowedge,
1404                              highedge);
1405     dxyNormWidthMap[i] =
1406         new TH2F(Form("norm_widths_dxy_map_%i", i),
1407                  "width(d_{xy}/#sigma_{d_{xy}}) map;track #eta;track #phi[rad];#sigma(d_{xy}/#sigma_{d_{xy}})",
1408                  nBins_,
1409                  lowedge,
1410                  highedge,
1411                  nBins_,
1412                  lowedge,
1413                  highedge);
1414     dzNormWidthMap[i] =
1415         new TH2F(Form("norm_widths_dz_map_%i", i),
1416                  "width(d_{z}/#sigma_{d_{z}}) map;track #eta;track #phi [rad];#sigma(d_{z}/#sigma_{d_{z}})",
1417                  nBins_,
1418                  lowedge,
1419                  highedge,
1420                  nBins_,
1421                  lowedge,
1422                  highedge);
1423 
1424     // 2D maps L1
1425     dxyMeanL1Map[i] = new TH2F(Form("means_dxy_L1Map_%i", i),
1426                                "#LT d_{xy} #GT map;module number;ladder number;#LT d_{xy} #GT [#mum]",
1427                                nModZ_,
1428                                -0.5,
1429                                nModZ_ - 0.5,
1430                                nLadders_,
1431                                -0.5,
1432                                nLadders_ - 0.5);
1433     dzMeanL1Map[i] = new TH2F(Form("means_dz_L1Map_%i", i),
1434                               "#LT d_{z} #GT map;module number;ladder number;#LT d_{z} #GT [#mum]",
1435                               nModZ_,
1436                               -0.5,
1437                               nModZ_ - 0.5,
1438                               nLadders_,
1439                               -0.5,
1440                               nLadders_ - 0.5);
1441     dxyNormMeanL1Map[i] =
1442         new TH2F(Form("norm_means_dxy_L1Map_%i", i),
1443                  "#LT d_{xy}/#sigma_{d_{xy}} #GT map;module number;ladder number;#LT d_{xy}/#sigma_{d_{xy}} #GT",
1444                  nModZ_,
1445                  -0.5,
1446                  nModZ_ - 0.5,
1447                  nLadders_,
1448                  -0.5,
1449                  nLadders_ - 0.5);
1450     dzNormMeanL1Map[i] =
1451         new TH2F(Form("norm_means_dz_L1Map_%i", i),
1452                  "#LT d_{z}/#sigma_{d_{z}} #GT map;module number;ladder number;#LT d_{xy}/#sigma_{d_{z}} #GT",
1453                  nModZ_,
1454                  -0.5,
1455                  nModZ_ - 0.5,
1456                  nLadders_,
1457                  -0.5,
1458                  nLadders_ - 0.5);
1459 
1460     dxyWidthL1Map[i] = new TH2F(Form("widths_dxy_L1Map_%i", i),
1461                                 "#sigma_{d_{xy}} map;module number;ladder number;#sigma(d_{xy}) [#mum]",
1462                                 nModZ_,
1463                                 -0.5,
1464                                 nModZ_ - 0.5,
1465                                 nLadders_,
1466                                 -0.5,
1467                                 nLadders_ - 0.5);
1468     dzWidthL1Map[i] = new TH2F(Form("widths_dz_L1Map_%i", i),
1469                                "#sigma_{d_{z}} map;module number;ladder number;#sigma(d_{z}) [#mum]",
1470                                nModZ_,
1471                                -0.5,
1472                                nModZ_ - 0.5,
1473                                nLadders_,
1474                                -0.5,
1475                                nLadders_ - 0.5);
1476     dxyNormWidthL1Map[i] =
1477         new TH2F(Form("norm_widths_dxy_L1Map_%i", i),
1478                  "width(d_{xy}/#sigma_{d_{xy}}) map;module number;ladder number;#sigma(d_{xy}/#sigma_{d_{xy}})",
1479                  nModZ_,
1480                  -0.5,
1481                  nModZ_ - 0.5,
1482                  nLadders_,
1483                  -0.5,
1484                  nLadders_ - 0.5);
1485     dzNormWidthL1Map[i] =
1486         new TH2F(Form("norm_widths_dz_L1Map_%i", i),
1487                  "width(d_{z}/#sigma_{d_{z}}) map;module number;ladder number;#sigma(d_{z}/#sigma_{d_{z}})",
1488                  nModZ_,
1489                  -0.5,
1490                  nModZ_ - 0.5,
1491                  nLadders_,
1492                  -0.5,
1493                  nLadders_ - 0.5);
1494 
1495     // DCA absolute
1496 
1497     if (isDebugMode) {
1498       timer.Stop();
1499       std::cout << "check point 3-" << i << " " << timer.CpuTime() << " " << timer.RealTime() << std::endl;
1500       timer.Continue();
1501     }
1502 
1503     FillTrendPlot(dxyPhiMeanTrend[i], dxyPhiResiduals[i], params::MEAN, "phi", nBins_);
1504     FillTrendPlot(dxyPhiWidthTrend[i], dxyPhiResiduals[i], params::WIDTH, "phi", nBins_);
1505 
1506     FillTrendPlot(dxPhiMeanTrend[i], dxPhiResiduals[i], params::MEAN, "phi", nBins_);
1507     FillTrendPlot(dxPhiWidthTrend[i], dxPhiResiduals[i], params::WIDTH, "phi", nBins_);
1508 
1509     FillTrendPlot(dyPhiMeanTrend[i], dyPhiResiduals[i], params::MEAN, "phi", nBins_);
1510     FillTrendPlot(dyPhiWidthTrend[i], dyPhiResiduals[i], params::WIDTH, "phi", nBins_);
1511 
1512     FillTrendPlot(dzPhiMeanTrend[i], dzPhiResiduals[i], params::MEAN, "phi", nBins_);
1513     FillTrendPlot(dzPhiWidthTrend[i], dzPhiResiduals[i], params::WIDTH, "phi", nBins_);
1514 
1515     FillTrendPlot(dxyEtaMeanTrend[i], dxyEtaResiduals[i], params::MEAN, "eta", nBins_);
1516     FillTrendPlot(dxyEtaWidthTrend[i], dxyEtaResiduals[i], params::WIDTH, "eta", nBins_);
1517 
1518     FillTrendPlot(dxEtaMeanTrend[i], dxEtaResiduals[i], params::MEAN, "eta", nBins_);
1519     FillTrendPlot(dxEtaWidthTrend[i], dxEtaResiduals[i], params::WIDTH, "eta", nBins_);
1520 
1521     FillTrendPlot(dyEtaMeanTrend[i], dyEtaResiduals[i], params::MEAN, "eta", nBins_);
1522     FillTrendPlot(dyEtaWidthTrend[i], dyEtaResiduals[i], params::WIDTH, "eta", nBins_);
1523 
1524     FillTrendPlot(dzEtaMeanTrend[i], dzEtaResiduals[i], params::MEAN, "eta", nBins_);
1525     FillTrendPlot(dzEtaWidthTrend[i], dzEtaResiduals[i], params::WIDTH, "eta", nBins_);
1526 
1527     if (minPt_ > 0.) {
1528       FillTrendPlot(dxyPtMeanTrend[i], dxyPtResiduals[i], params::MEAN, "pT", nPtBins_);
1529       FillTrendPlot(dxyPtWidthTrend[i], dxyPtResiduals[i], params::WIDTH, "pT", nPtBins_);
1530       FillTrendPlot(dzPtMeanTrend[i], dzPtResiduals[i], params::MEAN, "pT", nPtBins_);
1531       FillTrendPlot(dzPtWidthTrend[i], dzPtResiduals[i], params::WIDTH, "pT", nPtBins_);
1532     }
1533 
1534     if (nLadders_ > 0) {
1535       FillTrendPlot(dxyLadderMeanTrend[i], dxyLadderResiduals[i], params::MEAN, "else", nLadders_);
1536       FillTrendPlot(dxyLadderWidthTrend[i], dxyLadderResiduals[i], params::WIDTH, "else", nLadders_);
1537       FillTrendPlot(dzLadderMeanTrend[i], dzLadderResiduals[i], params::MEAN, "else", nLadders_);
1538       FillTrendPlot(dzLadderWidthTrend[i], dzLadderResiduals[i], params::WIDTH, "else", nLadders_);
1539     }
1540 
1541     if (nModZ_ > 0) {
1542       FillTrendPlot(dxyModZMeanTrend[i], dxyModZResiduals[i], params::MEAN, "else", nModZ_);
1543       FillTrendPlot(dxyModZWidthTrend[i], dxyModZResiduals[i], params::WIDTH, "else", nModZ_);
1544       FillTrendPlot(dzModZMeanTrend[i], dzModZResiduals[i], params::MEAN, "else", nModZ_);
1545       FillTrendPlot(dzModZWidthTrend[i], dzModZResiduals[i], params::WIDTH, "else", nModZ_);
1546     }
1547 
1548     MakeNiceTrendPlotStyle(dxyPhiMeanTrend[i], colors[i], markers[i]);
1549     MakeNiceTrendPlotStyle(dxyPhiWidthTrend[i], colors[i], markers[i]);
1550     MakeNiceTrendPlotStyle(dxPhiMeanTrend[i], colors[i], markers[i]);
1551     MakeNiceTrendPlotStyle(dxPhiWidthTrend[i], colors[i], markers[i]);
1552     MakeNiceTrendPlotStyle(dyPhiMeanTrend[i], colors[i], markers[i]);
1553     MakeNiceTrendPlotStyle(dyPhiWidthTrend[i], colors[i], markers[i]);
1554     MakeNiceTrendPlotStyle(dzPhiMeanTrend[i], colors[i], markers[i]);
1555     MakeNiceTrendPlotStyle(dzPhiWidthTrend[i], colors[i], markers[i]);
1556 
1557     MakeNiceTrendPlotStyle(dxyEtaMeanTrend[i], colors[i], markers[i]);
1558     MakeNiceTrendPlotStyle(dxyEtaWidthTrend[i], colors[i], markers[i]);
1559     MakeNiceTrendPlotStyle(dxEtaMeanTrend[i], colors[i], markers[i]);
1560     MakeNiceTrendPlotStyle(dxEtaWidthTrend[i], colors[i], markers[i]);
1561     MakeNiceTrendPlotStyle(dyEtaMeanTrend[i], colors[i], markers[i]);
1562     MakeNiceTrendPlotStyle(dyEtaWidthTrend[i], colors[i], markers[i]);
1563     MakeNiceTrendPlotStyle(dzEtaMeanTrend[i], colors[i], markers[i]);
1564     MakeNiceTrendPlotStyle(dzEtaWidthTrend[i], colors[i], markers[i]);
1565 
1566     if (minPt_ > 0.) {
1567       MakeNiceTrendPlotStyle(dxyPtMeanTrend[i], colors[i], markers[i]);
1568       MakeNiceTrendPlotStyle(dxyPtWidthTrend[i], colors[i], markers[i]);
1569       MakeNiceTrendPlotStyle(dzPtMeanTrend[i], colors[i], markers[i]);
1570       MakeNiceTrendPlotStyle(dzPtWidthTrend[i], colors[i], markers[i]);
1571     }
1572 
1573     if (nLadders_ > 0) {
1574       MakeNiceTrendPlotStyle(dxyLadderMeanTrend[i], colors[i], markers[i]);
1575       MakeNiceTrendPlotStyle(dxyLadderWidthTrend[i], colors[i], markers[i]);
1576       MakeNiceTrendPlotStyle(dzLadderMeanTrend[i], colors[i], markers[i]);
1577       MakeNiceTrendPlotStyle(dzLadderWidthTrend[i], colors[i], markers[i]);
1578     }
1579 
1580     if (nModZ_ > 0) {
1581       MakeNiceTrendPlotStyle(dxyModZMeanTrend[i], colors[i], markers[i]);
1582       MakeNiceTrendPlotStyle(dxyModZWidthTrend[i], colors[i], markers[i]);
1583       MakeNiceTrendPlotStyle(dzModZMeanTrend[i], colors[i], markers[i]);
1584       MakeNiceTrendPlotStyle(dzModZWidthTrend[i], colors[i], markers[i]);
1585     }
1586 
1587     // DCA normalized
1588 
1589     FillTrendPlot(dxyNormPhiMeanTrend[i], dxyNormPhiResiduals[i], params::MEAN, "phi", nBins_);
1590     FillTrendPlot(dxyNormPhiWidthTrend[i], dxyNormPhiResiduals[i], params::WIDTH, "phi", nBins_);
1591     FillTrendPlot(dzNormPhiMeanTrend[i], dzNormPhiResiduals[i], params::MEAN, "phi", nBins_);
1592     FillTrendPlot(dzNormPhiWidthTrend[i], dzNormPhiResiduals[i], params::WIDTH, "phi", nBins_);
1593 
1594     FillTrendPlot(dxyNormEtaMeanTrend[i], dxyNormEtaResiduals[i], params::MEAN, "eta", nBins_);
1595     FillTrendPlot(dxyNormEtaWidthTrend[i], dxyNormEtaResiduals[i], params::WIDTH, "eta", nBins_);
1596     FillTrendPlot(dzNormEtaMeanTrend[i], dzNormEtaResiduals[i], params::MEAN, "eta", nBins_);
1597     FillTrendPlot(dzNormEtaWidthTrend[i], dzNormEtaResiduals[i], params::WIDTH, "eta", nBins_);
1598 
1599     if (minPt_ > 0.) {
1600       FillTrendPlot(dxyNormPtMeanTrend[i], dxyNormPtResiduals[i], params::MEAN, "pT", nPtBins_);
1601       FillTrendPlot(dxyNormPtWidthTrend[i], dxyNormPtResiduals[i], params::WIDTH, "pT", nPtBins_);
1602       FillTrendPlot(dzNormPtMeanTrend[i], dzNormPtResiduals[i], params::MEAN, "pT", nPtBins_);
1603       FillTrendPlot(dzNormPtWidthTrend[i], dzNormPtResiduals[i], params::WIDTH, "pT", nPtBins_);
1604     }
1605 
1606     if (nLadders_ > 0) {
1607       FillTrendPlot(dxyNormLadderMeanTrend[i], dxyNormLadderResiduals[i], params::MEAN, "else", nLadders_);
1608       FillTrendPlot(dxyNormLadderWidthTrend[i], dxyNormLadderResiduals[i], params::WIDTH, "else", nLadders_);
1609       FillTrendPlot(dzNormLadderMeanTrend[i], dzNormLadderResiduals[i], params::MEAN, "else", nLadders_);
1610       FillTrendPlot(dzNormLadderWidthTrend[i], dzNormLadderResiduals[i], params::WIDTH, "else", nLadders_);
1611     }
1612 
1613     if (nModZ_ > 0) {
1614       FillTrendPlot(dxyNormModZMeanTrend[i], dxyNormModZResiduals[i], params::MEAN, "else", nModZ_);
1615       FillTrendPlot(dxyNormModZWidthTrend[i], dxyNormModZResiduals[i], params::WIDTH, "else", nModZ_);
1616       FillTrendPlot(dzNormModZMeanTrend[i], dzNormModZResiduals[i], params::MEAN, "else", nModZ_);
1617       FillTrendPlot(dzNormModZWidthTrend[i], dzNormModZResiduals[i], params::WIDTH, "else", nModZ_);
1618     }
1619 
1620     MakeNiceTrendPlotStyle(dxyNormPhiMeanTrend[i], colors[i], markers[i]);
1621     MakeNiceTrendPlotStyle(dxyNormPhiWidthTrend[i], colors[i], markers[i]);
1622     MakeNiceTrendPlotStyle(dzNormPhiMeanTrend[i], colors[i], markers[i]);
1623     MakeNiceTrendPlotStyle(dzNormPhiWidthTrend[i], colors[i], markers[i]);
1624 
1625     MakeNiceTrendPlotStyle(dxyNormEtaMeanTrend[i], colors[i], markers[i]);
1626     MakeNiceTrendPlotStyle(dxyNormEtaWidthTrend[i], colors[i], markers[i]);
1627     MakeNiceTrendPlotStyle(dzNormEtaMeanTrend[i], colors[i], markers[i]);
1628     MakeNiceTrendPlotStyle(dzNormEtaWidthTrend[i], colors[i], markers[i]);
1629 
1630     if (minPt_ > 0.) {
1631       MakeNiceTrendPlotStyle(dxyNormPtMeanTrend[i], colors[i], markers[i]);
1632       MakeNiceTrendPlotStyle(dxyNormPtWidthTrend[i], colors[i], markers[i]);
1633       MakeNiceTrendPlotStyle(dzNormPtMeanTrend[i], colors[i], markers[i]);
1634       MakeNiceTrendPlotStyle(dzNormPtWidthTrend[i], colors[i], markers[i]);
1635     }
1636 
1637     if (nLadders_ > 0) {
1638       MakeNiceTrendPlotStyle(dxyNormLadderMeanTrend[i], colors[i], markers[i]);
1639       MakeNiceTrendPlotStyle(dxyNormLadderWidthTrend[i], colors[i], markers[i]);
1640       MakeNiceTrendPlotStyle(dzNormLadderMeanTrend[i], colors[i], markers[i]);
1641       MakeNiceTrendPlotStyle(dzNormLadderWidthTrend[i], colors[i], markers[i]);
1642     }
1643 
1644     if (nModZ_ > 0) {
1645       MakeNiceTrendPlotStyle(dxyNormModZMeanTrend[i], colors[i], markers[i]);
1646       MakeNiceTrendPlotStyle(dxyNormModZWidthTrend[i], colors[i], markers[i]);
1647       MakeNiceTrendPlotStyle(dzNormModZMeanTrend[i], colors[i], markers[i]);
1648       MakeNiceTrendPlotStyle(dzNormModZWidthTrend[i], colors[i], markers[i]);
1649     }
1650 
1651     // maps
1652 
1653     //TTimeStamp filling1D_done;
1654 
1655     if (do2DMaps) {
1656       if (isDebugMode) {
1657         timer.Stop();
1658         std::cout << "check point 4-" << i << " " << timer.CpuTime() << " " << timer.RealTime() << std::endl;
1659         timer.Continue();
1660       }
1661 
1662       std::vector<std::vector<TH1F *> > v_dxyAbsMap;  //(nBins_, std::vector<TH1F*>(nBins_));
1663       std::vector<std::vector<TH1F *> > v_dzAbsMap;
1664       std::vector<std::vector<TH1F *> > v_dxyNormMap;
1665       std::vector<std::vector<TH1F *> > v_dzNormMap;
1666 
1667       for (Int_t index1 = 0; index1 < nBins_; index1++) {
1668         std::vector<TH1F *> a_temp_vec_xy;
1669         std::vector<TH1F *> n_temp_vec_xy;
1670         std::vector<TH1F *> a_temp_vec_z;
1671         std::vector<TH1F *> n_temp_vec_z;
1672 
1673         for (Int_t index2 = 0; index2 < nBins_; index2++) {
1674           if (isDebugMode)
1675             std::cout << index1 << " " << index2 << " " << (dxyMapResiduals[i][index1][index2])->GetName() << " "
1676                       << (dxyMapResiduals[i][index1][index2])->GetEntries() << std::endl;
1677 
1678           a_temp_vec_xy.push_back(dxyMapResiduals[i][index1][index2]);
1679           n_temp_vec_xy.push_back(dxyNormMapResiduals[i][index1][index2]);
1680           a_temp_vec_z.push_back(dzMapResiduals[i][index1][index2]);
1681           n_temp_vec_z.push_back(dzNormMapResiduals[i][index1][index2]);
1682         }
1683 
1684         v_dxyAbsMap.push_back(a_temp_vec_xy);
1685         v_dzAbsMap.push_back(a_temp_vec_z);
1686         v_dxyNormMap.push_back(n_temp_vec_xy);
1687         v_dzNormMap.push_back(n_temp_vec_z);
1688       }
1689 
1690       FillMap(dxyMeanMap[i], v_dxyAbsMap, params::MEAN);
1691       FillMap(dxyWidthMap[i], v_dxyAbsMap, params::WIDTH);
1692       FillMap(dzMeanMap[i], v_dzAbsMap, params::MEAN);
1693       FillMap(dzWidthMap[i], v_dzAbsMap, params::WIDTH);
1694 
1695       FillMap(dxyNormMeanMap[i], v_dxyNormMap, params::MEAN);
1696       FillMap(dxyNormWidthMap[i], v_dxyNormMap, params::WIDTH);
1697       FillMap(dzNormMeanMap[i], v_dzNormMap, params::MEAN);
1698       FillMap(dzNormWidthMap[i], v_dzNormMap, params::WIDTH);
1699 
1700       if (isDebugMode) {
1701         timer.Stop();
1702         std::cout << "check point 5-" << i << " " << timer.CpuTime() << " " << timer.RealTime() << std::endl;
1703         timer.Continue();
1704       }
1705 
1706       t_dxyMeanMap[i] = trimTheMap(dxyMeanMap[i]).first;
1707       t_dxyWidthMap[i] = trimTheMap(dxyWidthMap[i]).first;
1708       t_dzMeanMap[i] = trimTheMap(dzMeanMap[i]).first;
1709       t_dzWidthMap[i] = trimTheMap(dzWidthMap[i]).first;
1710 
1711       t_dxyNormMeanMap[i] = trimTheMap(dxyNormMeanMap[i]).first;
1712       t_dxyNormWidthMap[i] = trimTheMap(dxyNormWidthMap[i]).first;
1713       t_dzNormMeanMap[i] = trimTheMap(dzNormMeanMap[i]).first;
1714       t_dzNormWidthMap[i] = trimTheMap(dzNormWidthMap[i]).first;
1715 
1716       MakeNiceMapStyle(t_dxyMeanMap[i]);
1717       MakeNiceMapStyle(t_dxyWidthMap[i]);
1718       MakeNiceMapStyle(t_dzMeanMap[i]);
1719       MakeNiceMapStyle(t_dzWidthMap[i]);
1720 
1721       MakeNiceMapStyle(t_dxyNormMeanMap[i]);
1722       MakeNiceMapStyle(t_dxyNormWidthMap[i]);
1723       MakeNiceMapStyle(t_dzNormMeanMap[i]);
1724       MakeNiceMapStyle(t_dzNormWidthMap[i]);
1725 
1726       // clear the vectors of vectors
1727       for (Int_t index1 = 0; index1 < nBins_; index1++) {
1728         v_dxyAbsMap.clear();
1729         v_dzAbsMap.clear();
1730         v_dxyNormMap.clear();
1731         v_dzNormMap.clear();
1732       }
1733 
1734       for (Int_t index1 = 0; index1 < nLadders_; index1++) {
1735         std::vector<TH1F *> a_temp_vec_xy;
1736         std::vector<TH1F *> n_temp_vec_xy;
1737         std::vector<TH1F *> a_temp_vec_z;
1738         std::vector<TH1F *> n_temp_vec_z;
1739 
1740         for (Int_t index2 = 0; index2 < nModZ_; index2++) {
1741           a_temp_vec_xy.push_back(dxyL1MapResiduals[i][index1][index2]);
1742           a_temp_vec_z.push_back(dzL1MapResiduals[i][index1][index2]);
1743           n_temp_vec_xy.push_back(dxyL1NormMapResiduals[i][index1][index2]);
1744           n_temp_vec_z.push_back(dzL1NormMapResiduals[i][index1][index2]);
1745         }
1746 
1747         v_dxyAbsMap.push_back(a_temp_vec_xy);
1748         v_dzAbsMap.push_back(a_temp_vec_z);
1749         v_dxyNormMap.push_back(n_temp_vec_xy);
1750         v_dzNormMap.push_back(n_temp_vec_z);
1751       }
1752 
1753       FillMap(dxyMeanL1Map[i], v_dxyAbsMap, params::MEAN, nModZ_, nLadders_);
1754       FillMap(dxyWidthL1Map[i], v_dxyAbsMap, params::WIDTH, nModZ_, nLadders_);
1755       FillMap(dzMeanL1Map[i], v_dzAbsMap, params::MEAN, nModZ_, nLadders_);
1756       FillMap(dzWidthL1Map[i], v_dzAbsMap, params::WIDTH, nModZ_, nLadders_);
1757 
1758       FillMap(dxyNormMeanL1Map[i], v_dxyNormMap, params::MEAN, nModZ_, nLadders_);
1759       FillMap(dxyNormWidthL1Map[i], v_dxyNormMap, params::WIDTH, nModZ_, nLadders_);
1760       FillMap(dzNormMeanL1Map[i], v_dzNormMap, params::MEAN, nModZ_, nLadders_);
1761       FillMap(dzNormWidthL1Map[i], v_dzNormMap, params::WIDTH, nModZ_, nLadders_);
1762 
1763       if (isDebugMode) {
1764         timer.Stop();
1765         std::cout << "check point 5-" << i << " " << timer.CpuTime() << " " << timer.RealTime() << std::endl;
1766         timer.Continue();
1767       }
1768 
1769       t_dxyMeanL1Map[i] = trimTheMap(dxyMeanL1Map[i]).first;
1770       t_dxyWidthL1Map[i] = trimTheMap(dxyWidthL1Map[i]).first;
1771       t_dzMeanL1Map[i] = trimTheMap(dzMeanL1Map[i]).first;
1772       t_dzWidthL1Map[i] = trimTheMap(dzWidthL1Map[i]).first;
1773 
1774       t_dxyNormMeanL1Map[i] = trimTheMap(dxyNormMeanL1Map[i]).first;
1775       t_dxyNormWidthL1Map[i] = trimTheMap(dxyNormWidthL1Map[i]).first;
1776       t_dzNormMeanL1Map[i] = trimTheMap(dzNormMeanL1Map[i]).first;
1777       t_dzNormWidthL1Map[i] = trimTheMap(dzNormWidthL1Map[i]).first;
1778 
1779       MakeNiceMapStyle(t_dxyMeanL1Map[i]);
1780       MakeNiceMapStyle(t_dxyWidthL1Map[i]);
1781       MakeNiceMapStyle(t_dzMeanL1Map[i]);
1782       MakeNiceMapStyle(t_dzWidthL1Map[i]);
1783 
1784       MakeNiceMapStyle(t_dxyNormMeanL1Map[i]);
1785       MakeNiceMapStyle(t_dxyNormWidthL1Map[i]);
1786       MakeNiceMapStyle(t_dzNormMeanL1Map[i]);
1787       MakeNiceMapStyle(t_dzNormWidthL1Map[i]);
1788 
1789     }  // if do2DMaps
1790 
1791     MakeNiceTrendPlotStyle(dxyRefit[i], colors[i], markers[i]);
1792     MakeNiceTrendPlotStyle(dzRefit[i], colors[i], markers[i]);
1793     MakeNiceTrendPlotStyle(dxySigRefit[i], colors[i], markers[i]);
1794     MakeNiceTrendPlotStyle(dzSigRefit[i], colors[i], markers[i]);
1795   }
1796 
1797   TTimeStamp filling2D_done;
1798 
1799   TString theStrDate = theDate;
1800   TString theStrAlignment = LegLabels[0];
1801 
1802   /*
1803     // in case labels are needed
1804     std::vector<TString> vLabels(LegLabels, LegLabels+10);
1805     vLabels.shrink_to_fit();
1806   */
1807 
1808   for (Int_t j = 1; j < nFiles_; j++) {
1809     theStrAlignment += ("_vs_" + LegLabels[j]);
1810   }
1811 
1812   theStrDate.ReplaceAll(" ", "");
1813   theStrAlignment.ReplaceAll(" ", "_");
1814 
1815   // non-differential
1816   TCanvas *BareResiduals = new TCanvas("BareResiduals", "BareResiduals", 1200, 1200);
1817   arrangeBiasCanvas(BareResiduals, dxyRefit, dxySigRefit, dzRefit, dzSigRefit, nFiles_, LegLabels, theDate, true);
1818 
1819   BareResiduals->SaveAs("ResidualsCanvas_" + theStrDate + theStrAlignment + ".pdf");
1820   BareResiduals->SaveAs("ResidualsCanvas_" + theStrDate + theStrAlignment + ".png");
1821 
1822   // DCA absolute
1823 
1824   TCanvas *dxyPhiTrend = new TCanvas("dxyPhiTrend", "dxyPhiTrend", 1200, 600);
1825   arrangeCanvas(dxyPhiTrend, dxyPhiMeanTrend, dxyPhiWidthTrend, nFiles_, LegLabels, theDate, false, setAutoLimits);
1826 
1827   dxyPhiTrend->SaveAs("dxyPhiTrend_" + theStrDate + theStrAlignment + ".pdf");
1828   dxyPhiTrend->SaveAs("dxyPhiTrend_" + theStrDate + theStrAlignment + ".png");
1829 
1830   TCanvas *dzPhiTrend = new TCanvas("dzPhiTrend", "dzPhiTrend", 1200, 600);
1831   arrangeCanvas(dzPhiTrend, dzPhiMeanTrend, dzPhiWidthTrend, nFiles_, LegLabels, theDate, false, setAutoLimits);
1832 
1833   dzPhiTrend->SaveAs("dzPhiTrend_" + theStrDate + theStrAlignment + ".pdf");
1834   dzPhiTrend->SaveAs("dzPhiTrend_" + theStrDate + theStrAlignment + ".png");
1835 
1836   TCanvas *dxyEtaTrend = new TCanvas("dxyEtaTrend", "dxyEtaTrend", 1200, 600);
1837   arrangeCanvas(dxyEtaTrend, dxyEtaMeanTrend, dxyEtaWidthTrend, nFiles_, LegLabels, theDate, false, setAutoLimits);
1838 
1839   dxyEtaTrend->SaveAs("dxyEtaTrend_" + theStrDate + theStrAlignment + ".pdf");
1840   dxyEtaTrend->SaveAs("dxyEtaTrend_" + theStrDate + theStrAlignment + ".png");
1841 
1842   TCanvas *dzEtaTrend = new TCanvas("dzEtaTrend", "dzEtaTrend", 1200, 600);
1843   arrangeCanvas(dzEtaTrend, dzEtaMeanTrend, dzEtaWidthTrend, nFiles_, LegLabels, theDate, false, setAutoLimits);
1844 
1845   dzEtaTrend->SaveAs("dzEtaTrend_" + theStrDate + theStrAlignment + ".pdf");
1846   dzEtaTrend->SaveAs("dzEtaTrend_" + theStrDate + theStrAlignment + ".png");
1847 
1848   if (nLadders_ > 0) {
1849     TCanvas *dxyLadderTrend = new TCanvas("dxyLadderTrend", "dxyLadderTrend", 600, 600);
1850     arrangeCanvas(
1851         dxyLadderTrend, dxyLadderMeanTrend, dxyLadderWidthTrend, nFiles_, LegLabels, theDate, true, setAutoLimits);
1852 
1853     dxyLadderTrend->SaveAs("dxyLadderTrend_" + theStrDate + theStrAlignment + ".pdf");
1854     dxyLadderTrend->SaveAs("dxyLadderTrend_" + theStrDate + theStrAlignment + ".png");
1855 
1856     delete dxyLadderTrend;
1857   }
1858 
1859   // fit dz vs phi
1860   TCanvas *dzPhiTrendFit = new TCanvas("dzPhiTrendFit", "dzPhiTrendFit", 1200, 600);
1861   arrangeFitCanvas(dzPhiTrendFit, dzPhiMeanTrend, nFiles_, LegLabels, theDate);
1862 
1863   dzPhiTrendFit->SaveAs("dzPhiTrendFit_" + theStrDate + theStrAlignment + ".pdf");
1864   dzPhiTrendFit->SaveAs("dzPhiTrendFit_" + theStrDate + theStrAlignment + ".png");
1865 
1866   if (minPt_ > 0.) {
1867     TCanvas *dxyPtTrend = new TCanvas("dxyPtTrend", "dxyPtTrend", 1200, 600);
1868     arrangeCanvas(dxyPtTrend, dxyPtMeanTrend, dxyPtWidthTrend, nFiles_, LegLabels, theDate, false, setAutoLimits);
1869 
1870     dxyPtTrend->SaveAs("dxyPtTrend_" + theStrDate + theStrAlignment + ".pdf");
1871     dxyPtTrend->SaveAs("dxyPtTrend_" + theStrDate + theStrAlignment + ".png");
1872 
1873     TCanvas *dzPtTrend = new TCanvas("dzPtTrend", "dzPtTrend", 1200, 600);
1874     arrangeCanvas(dzPtTrend, dzPtMeanTrend, dzPtWidthTrend, nFiles_, LegLabels, theDate, false, setAutoLimits);
1875 
1876     dzPtTrend->SaveAs("dzPtTrend_" + theStrDate + theStrAlignment + ".pdf");
1877     dzPtTrend->SaveAs("dzPtTrend_" + theStrDate + theStrAlignment + ".png");
1878 
1879     delete dxyPtTrend;
1880     delete dzPtTrend;
1881   }
1882 
1883   // delete all news
1884 
1885   delete BareResiduals;
1886   delete dxyPhiTrend;
1887   delete dzPhiTrend;
1888   delete dxyEtaTrend;
1889   delete dzEtaTrend;
1890   delete dzPhiTrendFit;
1891 
1892   // DCA normalized
1893 
1894   TCanvas *dxyNormPhiTrend = new TCanvas("dxyNormPhiTrend", "dxyNormPhiTrend", 1200, 600);
1895   arrangeCanvas(
1896       dxyNormPhiTrend, dxyNormPhiMeanTrend, dxyNormPhiWidthTrend, nFiles_, LegLabels, theDate, false, setAutoLimits);
1897 
1898   dxyNormPhiTrend->SaveAs("dxyPhiTrendNorm_" + theStrDate + theStrAlignment + ".pdf");
1899   dxyNormPhiTrend->SaveAs("dxyPhiTrendNorm_" + theStrDate + theStrAlignment + ".png");
1900 
1901   TCanvas *dzNormPhiTrend = new TCanvas("dzNormPhiTrend", "dzNormPhiTrend", 1200, 600);
1902   arrangeCanvas(
1903       dzNormPhiTrend, dzNormPhiMeanTrend, dzNormPhiWidthTrend, nFiles_, LegLabels, theDate, false, setAutoLimits);
1904 
1905   dzNormPhiTrend->SaveAs("dzPhiTrendNorm_" + theStrDate + theStrAlignment + ".pdf");
1906   dzNormPhiTrend->SaveAs("dzPhiTrendNorm_" + theStrDate + theStrAlignment + ".png");
1907 
1908   TCanvas *dxyNormEtaTrend = new TCanvas("dxyNormEtaTrend", "dxyNormEtaTrend", 1200, 600);
1909   arrangeCanvas(
1910       dxyNormEtaTrend, dxyNormEtaMeanTrend, dxyNormEtaWidthTrend, nFiles_, LegLabels, theDate, false, setAutoLimits);
1911 
1912   dxyNormEtaTrend->SaveAs("dxyEtaTrendNorm_" + theStrDate + theStrAlignment + ".pdf");
1913   dxyNormEtaTrend->SaveAs("dxyEtaTrendNorm_" + theStrDate + theStrAlignment + ".png");
1914 
1915   TCanvas *dzNormEtaTrend = new TCanvas("dzNormEtaTrend", "dzNormEtaTrend", 1200, 600);
1916   arrangeCanvas(
1917       dzNormEtaTrend, dzNormEtaMeanTrend, dzNormEtaWidthTrend, nFiles_, LegLabels, theDate, false, setAutoLimits);
1918 
1919   dzNormEtaTrend->SaveAs("dzEtaTrendNorm_" + theStrDate + theStrAlignment + ".pdf");
1920   dzNormEtaTrend->SaveAs("dzEtaTrendNorm_" + theStrDate + theStrAlignment + ".png");
1921 
1922   if (minPt_ > 0.) {
1923     TCanvas *dxyNormPtTrend = new TCanvas("dxyNormPtTrend", "dxyNormPtTrend", 1200, 600);
1924     arrangeCanvas(
1925         dxyNormPtTrend, dxyNormPtMeanTrend, dxyNormPtWidthTrend, nFiles_, LegLabels, theDate, false, setAutoLimits);
1926 
1927     dxyNormPtTrend->SaveAs("dxyPtTrendNorm_" + theStrDate + theStrAlignment + ".pdf");
1928     dxyNormPtTrend->SaveAs("dxyPtTrendNorm_" + theStrDate + theStrAlignment + ".png");
1929 
1930     TCanvas *dzNormPtTrend = new TCanvas("dzNormPtTrend", "dzNormPtTrend", 1200, 600);
1931     arrangeCanvas(
1932         dzNormPtTrend, dzNormPtMeanTrend, dzNormPtWidthTrend, nFiles_, LegLabels, theDate, false, setAutoLimits);
1933 
1934     dzNormPtTrend->SaveAs("dzPtTrendNorm_" + theStrDate + theStrAlignment + ".pdf");
1935     dzNormPtTrend->SaveAs("dzPtTrendNorm_" + theStrDate + theStrAlignment + ".png");
1936 
1937     delete dxyNormPtTrend;
1938     delete dzNormPtTrend;
1939   }
1940 
1941   // delete all news
1942 
1943   delete dxyNormPhiTrend;
1944   delete dzNormPhiTrend;
1945   delete dxyNormEtaTrend;
1946   delete dzNormEtaTrend;
1947 
1948   // Bias plots
1949 
1950   TCanvas *BiasesCanvas = new TCanvas("BiasCanvas", "BiasCanvas", 1200, 1200);
1951   arrangeBiasCanvas(BiasesCanvas,
1952                     dxyPhiMeanTrend,
1953                     dzPhiMeanTrend,
1954                     dxyEtaMeanTrend,
1955                     dzEtaMeanTrend,
1956                     nFiles_,
1957                     LegLabels,
1958                     theDate,
1959                     setAutoLimits);
1960 
1961   BiasesCanvas->SaveAs("BiasesCanvas_" + theStrDate + theStrAlignment + ".pdf");
1962   BiasesCanvas->SaveAs("BiasesCanvas_" + theStrDate + theStrAlignment + ".png");
1963 
1964   // Bias plots (x and y)
1965 
1966   TCanvas *BiasesCanvasXY = new TCanvas("BiasCanvasXY", "BiasCanvasXY", 1200, 1200);
1967   arrangeBiasCanvas(BiasesCanvasXY,
1968                     dxPhiMeanTrend,
1969                     dyPhiMeanTrend,
1970                     dxEtaMeanTrend,
1971                     dyEtaMeanTrend,
1972                     nFiles_,
1973                     LegLabels,
1974                     theDate,
1975                     setAutoLimits);
1976 
1977   BiasesCanvasXY->SaveAs("BiasesCanvasXY_" + theStrDate + theStrAlignment + ".pdf");
1978   BiasesCanvasXY->SaveAs("BiasesCanvasXY_" + theStrDate + theStrAlignment + ".png");
1979 
1980   // Bias plots (ladders and module number)
1981   if (nLadders_ > 0 && nModZ_ > 0) {
1982     TCanvas *BiasesCanvasLayer1 = new TCanvas("BiasCanvasLayer1", "BiasCanvasLayer1", 1200, 1200);
1983     arrangeBiasCanvas(BiasesCanvasLayer1,
1984                       dxyLadderMeanTrend,
1985                       dzLadderMeanTrend,
1986                       dxyModZMeanTrend,
1987                       dzModZMeanTrend,
1988                       nFiles_,
1989                       LegLabels,
1990                       theDate,
1991                       setAutoLimits);
1992 
1993     BiasesCanvasLayer1->SaveAs("BiasesCanvasLayer1_" + theStrDate + theStrAlignment + ".pdf");
1994     BiasesCanvasLayer1->SaveAs("BiasesCanvasLayer1_" + theStrDate + theStrAlignment + ".png");
1995     delete BiasesCanvasLayer1;
1996   }
1997 
1998   TCanvas *dxyPhiBiasCanvas = new TCanvas("dxyPhiBiasCanvas", "dxyPhiBiasCanvas", 600, 600);
1999   TCanvas *dxyEtaBiasCanvas = new TCanvas("dxyEtaBiasCanvas", "dxyEtaBiasCanvas", 600, 600);
2000   TCanvas *dzPhiBiasCanvas = new TCanvas("dzPhiBiasCanvas", "dzPhiBiasCanvas", 600, 600);
2001   TCanvas *dzEtaBiasCanvas = new TCanvas("dzEtaBiasCanvas", "dzEtaBiasCanvas", 600, 600);
2002 
2003   arrangeCanvas(dxyPhiBiasCanvas, dxyPhiMeanTrend, dxyPhiWidthTrend, nFiles_, LegLabels, theDate, true, setAutoLimits);
2004   arrangeCanvas(dzPhiBiasCanvas, dzPhiMeanTrend, dzPhiWidthTrend, nFiles_, LegLabels, theDate, true, setAutoLimits);
2005   arrangeCanvas(dxyEtaBiasCanvas, dxyEtaMeanTrend, dxyEtaWidthTrend, nFiles_, LegLabels, theDate, true, setAutoLimits);
2006   arrangeCanvas(dzEtaBiasCanvas, dzEtaMeanTrend, dzEtaWidthTrend, nFiles_, LegLabels, theDate, true, setAutoLimits);
2007 
2008   dxyPhiBiasCanvas->SaveAs("dxyPhiBiasCanvas_" + theStrDate + theStrAlignment + ".pdf");
2009   dxyEtaBiasCanvas->SaveAs("dxyEtaBiasCanvas_" + theStrDate + theStrAlignment + ".pdf");
2010   dzPhiBiasCanvas->SaveAs("dzPhiBiasCanvas_" + theStrDate + theStrAlignment + ".pdf");
2011   dzEtaBiasCanvas->SaveAs("dzEtaBiasCanvas_" + theStrDate + theStrAlignment + ".pdf");
2012 
2013   dxyPhiBiasCanvas->SaveAs("dxyPhiBiasCanvas_" + theStrDate + theStrAlignment + ".png");
2014   dxyEtaBiasCanvas->SaveAs("dxyEtaBiasCanvas_" + theStrDate + theStrAlignment + ".png");
2015   dzPhiBiasCanvas->SaveAs("dzPhiBiasCanvas_" + theStrDate + theStrAlignment + ".png");
2016   dzEtaBiasCanvas->SaveAs("dzEtaBiasCanvas_" + theStrDate + theStrAlignment + ".png");
2017 
2018   // delete all news
2019 
2020   delete BiasesCanvas;
2021   delete BiasesCanvasXY;
2022   delete dxyPhiBiasCanvas;
2023   delete dxyEtaBiasCanvas;
2024   delete dzPhiBiasCanvas;
2025   delete dzEtaBiasCanvas;
2026 
2027   // Resolution plots
2028   TCanvas *ResolutionsCanvas = new TCanvas("ResolutionsCanvas", "ResolutionsCanvas", 1200, 1200);
2029   arrangeBiasCanvas(ResolutionsCanvas,
2030                     dxyPhiWidthTrend,
2031                     dzPhiWidthTrend,
2032                     dxyEtaWidthTrend,
2033                     dzEtaWidthTrend,
2034                     nFiles_,
2035                     LegLabels,
2036                     theDate,
2037                     setAutoLimits);
2038 
2039   ResolutionsCanvas->SaveAs("ResolutionsCanvas_" + theStrDate + theStrAlignment + ".pdf");
2040   ResolutionsCanvas->SaveAs("ResolutionsCanvas_" + theStrDate + theStrAlignment + ".png");
2041 
2042   TCanvas *ResolutionsCanvasXY = new TCanvas("ResolutionsCanvasXY", "ResolutionsCanvasXY", 1200, 1200);
2043   arrangeBiasCanvas(ResolutionsCanvasXY,
2044                     dxPhiWidthTrend,
2045                     dyPhiWidthTrend,
2046                     dxEtaWidthTrend,
2047                     dyEtaWidthTrend,
2048                     nFiles_,
2049                     LegLabels,
2050                     theDate,
2051                     setAutoLimits);
2052 
2053   ResolutionsCanvasXY->SaveAs("ResolutionsCanvasXY_" + theStrDate + theStrAlignment + ".pdf");
2054   ResolutionsCanvasXY->SaveAs("ResolutionsCanvasXY_" + theStrDate + theStrAlignment + ".png");
2055 
2056   if (nLadders_ > 0 && nModZ_ > 0) {
2057     TCanvas *ResolutionsCanvasLayer1 = new TCanvas("ResolutionsCanvasLayer1", "ResolutionsCanvasLayer1", 1200, 1200);
2058     arrangeBiasCanvas(ResolutionsCanvasLayer1,
2059                       dxyLadderWidthTrend,
2060                       dzLadderWidthTrend,
2061                       dxyModZWidthTrend,
2062                       dzModZWidthTrend,
2063                       nFiles_,
2064                       LegLabels,
2065                       theDate,
2066                       setAutoLimits);
2067 
2068     ResolutionsCanvasLayer1->SaveAs("ResolutionsCanvasLayer1_" + theStrDate + theStrAlignment + ".pdf");
2069     ResolutionsCanvasLayer1->SaveAs("ResolutionsCanvasLayer1_" + theStrDate + theStrAlignment + ".png");
2070     delete ResolutionsCanvasLayer1;
2071   }
2072 
2073   // Pull plots
2074   TCanvas *PullsCanvas = new TCanvas("PullsCanvas", "PullsCanvas", 1200, 1200);
2075   arrangeBiasCanvas(PullsCanvas,
2076                     dxyNormPhiWidthTrend,
2077                     dzNormPhiWidthTrend,
2078                     dxyNormEtaWidthTrend,
2079                     dzNormEtaWidthTrend,
2080                     nFiles_,
2081                     LegLabels,
2082                     theDate,
2083                     setAutoLimits);
2084 
2085   PullsCanvas->SaveAs("PullsCanvas_" + theStrDate + theStrAlignment + ".pdf");
2086   PullsCanvas->SaveAs("PullsCanvas_" + theStrDate + theStrAlignment + ".png");
2087 
2088   if (nLadders_ > 0 && nModZ_ > 0) {
2089     TCanvas *PullsCanvasLayer1 = new TCanvas("PullsCanvasLayer1", "PullsCanvasLayer1", 1200, 1200);
2090     arrangeBiasCanvas(PullsCanvasLayer1,
2091                       dxyNormLadderWidthTrend,
2092                       dzNormLadderWidthTrend,
2093                       dxyNormModZWidthTrend,
2094                       dzNormModZWidthTrend,
2095                       nFiles_,
2096                       LegLabels,
2097                       theDate,
2098                       setAutoLimits);
2099 
2100     PullsCanvasLayer1->SaveAs("PullsCanvasLayer1_" + theStrDate + theStrAlignment + ".pdf");
2101     PullsCanvasLayer1->SaveAs("PullsCanvasLayer1_" + theStrDate + theStrAlignment + ".png");
2102     delete PullsCanvasLayer1;
2103   }
2104 
2105   // delete all news
2106   delete ResolutionsCanvas;
2107   delete ResolutionsCanvasXY;
2108   delete PullsCanvas;
2109 
2110   // 2D Maps
2111 
2112   if (do2DMaps) {
2113     TCanvas *dxyAbsMap = new TCanvas("dxyAbsMap", "dxyAbsMap", 1200, 500 * nFiles_);
2114     arrangeCanvas2D(dxyAbsMap, t_dxyMeanMap, t_dxyWidthMap, nFiles_, LegLabels, theDate);
2115     dxyAbsMap->SaveAs("dxyAbsMap_" + theStrDate + theStrAlignment + ".pdf");
2116     dxyAbsMap->SaveAs("dxyAbsMap_" + theStrDate + theStrAlignment + ".png");
2117 
2118     TCanvas *dzAbsMap = new TCanvas("dzAbsMap", "dzAbsMap", 1200, 500 * nFiles_);
2119     arrangeCanvas2D(dzAbsMap, t_dzMeanMap, t_dzWidthMap, nFiles_, LegLabels, theDate);
2120     dzAbsMap->SaveAs("dzAbsMap_" + theStrDate + theStrAlignment + ".pdf");
2121     dzAbsMap->SaveAs("dzAbsMap_" + theStrDate + theStrAlignment + ".png");
2122 
2123     TCanvas *dxyNormMap = new TCanvas("dxyNormMap", "dxyNormMap", 1200, 500 * nFiles_);
2124     arrangeCanvas2D(dxyNormMap, t_dxyNormMeanMap, t_dxyNormWidthMap, nFiles_, LegLabels, theDate);
2125     dxyNormMap->SaveAs("dxyNormMap_" + theStrDate + theStrAlignment + ".pdf");
2126     dxyNormMap->SaveAs("dxyNormMap_" + theStrDate + theStrAlignment + ".png");
2127 
2128     TCanvas *dzNormMap = new TCanvas("dzNormMap", "dzNormMap", 1200, 500 * nFiles_);
2129     arrangeCanvas2D(dzNormMap, t_dzNormMeanMap, t_dzNormWidthMap, nFiles_, LegLabels, theDate);
2130     dzNormMap->SaveAs("dzNormMap_" + theStrDate + theStrAlignment + ".pdf");
2131     dzNormMap->SaveAs("dzNormMap_" + theStrDate + theStrAlignment + ".png");
2132 
2133     delete dxyAbsMap;
2134     delete dzAbsMap;
2135     delete dxyNormMap;
2136     delete dzNormMap;
2137 
2138     // L1 Map
2139 
2140     TCanvas *dxyAbsL1Map = new TCanvas("dxyAbsL1Map", "dxyAbsL1Map", 1200, 500 * nFiles_);
2141     arrangeCanvas2D(dxyAbsL1Map, t_dxyMeanL1Map, t_dxyWidthL1Map, nFiles_, LegLabels, theDate);
2142     dxyAbsL1Map->SaveAs("dxyAbsL1Map_" + theStrDate + theStrAlignment + ".pdf");
2143     dxyAbsL1Map->SaveAs("dxyAbsL1Map_" + theStrDate + theStrAlignment + ".png");
2144 
2145     TCanvas *dzAbsL1Map = new TCanvas("dzAbsL1Map", "dzAbsL1Map", 1200, 500 * nFiles_);
2146     arrangeCanvas2D(dzAbsL1Map, t_dzMeanL1Map, t_dzWidthL1Map, nFiles_, LegLabels, theDate);
2147     dzAbsL1Map->SaveAs("dzAbsL1Map_" + theStrDate + theStrAlignment + ".pdf");
2148     dzAbsL1Map->SaveAs("dzAbsL1Map_" + theStrDate + theStrAlignment + ".png");
2149 
2150     TCanvas *dxyNormL1Map = new TCanvas("dxyNormL1Map", "dxyNormL1Map", 1200, 500 * nFiles_);
2151     arrangeCanvas2D(dxyNormL1Map, t_dxyNormMeanL1Map, t_dxyNormWidthL1Map, nFiles_, LegLabels, theDate);
2152     dxyNormL1Map->SaveAs("dxyNormL1Map_" + theStrDate + theStrAlignment + ".pdf");
2153     dxyNormL1Map->SaveAs("dxyNormL1Map_" + theStrDate + theStrAlignment + ".png");
2154 
2155     TCanvas *dzNormL1Map = new TCanvas("dzNormL1Map", "dzNormL1Map", 1200, 500 * nFiles_);
2156     arrangeCanvas2D(dzNormL1Map, t_dzNormMeanL1Map, t_dzNormWidthL1Map, nFiles_, LegLabels, theDate);
2157     dzNormL1Map->SaveAs("dzNormL1Map_" + theStrDate + theStrAlignment + ".pdf");
2158     dzNormL1Map->SaveAs("dzNormL1Map_" + theStrDate + theStrAlignment + ".png");
2159 
2160     delete dxyAbsL1Map;
2161     delete dzAbsL1Map;
2162     delete dxyNormL1Map;
2163     delete dzNormL1Map;
2164   }
2165 
2166   delete thePlotLimits;
2167 
2168   // delete everything in the source list
2169   for (std::vector<PVValidationVariables *>::iterator it = sourceList.begin(); it != sourceList.end(); ++it) {
2170     delete (*it);
2171   }
2172 
2173   TTimeStamp plotting_done;
2174 
2175   std::cout << " ======   TIMING REPORT ====== " << std::endl;
2176   std::cout << "time tp initialize = " << initialization_done.AsDouble() - start_time.AsDouble() << "s" << std::endl;
2177   std::cout << "time to cache      = " << caching_done.AsDouble() - initialization_done.AsDouble() << "s" << std::endl;
2178   //  std::cout<<"time to fill 1D    = "<<filling1D_done.AsDouble()-caching_done.AsDouble()<<"s"<<std::endl;
2179   std::cout << "time to fit        = " << filling2D_done.AsDouble() - caching_done.AsDouble() << "s" << std::endl;
2180   std::cout << "time to plot       = " << plotting_done.AsDouble() - filling2D_done.AsDouble() << "s" << std::endl;
2181 
2182   timer.Stop();
2183   timer.Print();
2184 }
2185 
2186 //*************************************************************
2187 void arrangeBiasCanvas(TCanvas *canv,
2188                        TH1F *dxyPhiMeanTrend[100],
2189                        TH1F *dzPhiMeanTrend[100],
2190                        TH1F *dxyEtaMeanTrend[100],
2191                        TH1F *dzEtaMeanTrend[100],
2192                        Int_t nFiles,
2193                        TString LegLabels[10],
2194                        TString theDate,
2195                        bool setAutoLimits) {
2196   //*************************************************************
2197 
2198   TLegend *lego = new TLegend(0.22, 0.80, 0.79, 0.91);
2199   // might be useful if many objects are compared
2200   if (nFiles > 3) {
2201     lego->SetNColumns(2);
2202   }
2203 
2204   lego->SetFillColor(10);
2205   if (nFiles > 3) {
2206     lego->SetTextSize(0.032);
2207   } else {
2208     lego->SetTextSize(0.042);
2209   }
2210   lego->SetTextFont(42);
2211   lego->SetFillColor(10);
2212   lego->SetLineColor(10);
2213   lego->SetShadowColor(10);
2214 
2215   TPaveText *ptDate = new TPaveText(0.20, 0.95, 0.50, 0.99, "blNDC");
2216   //ptDate->SetFillColor(kYellow);
2217   ptDate->SetFillColor(10);
2218   ptDate->SetBorderSize(1);
2219   ptDate->SetLineColor(kBlue);
2220   ptDate->SetLineWidth(1);
2221   ptDate->SetTextFont(32);
2222   TText *textDate = ptDate->AddText(theDate);
2223   textDate->SetTextSize(0.04);
2224   textDate->SetTextColor(kBlue);
2225 
2226   canv->SetFillColor(10);
2227   canv->Divide(2, 2);
2228 
2229   canv->cd(1)->SetBottomMargin(0.14);
2230   canv->cd(1)->SetLeftMargin(0.18);
2231   canv->cd(1)->SetRightMargin(0.01);
2232   canv->cd(1)->SetTopMargin(0.06);
2233 
2234   canv->cd(2)->SetBottomMargin(0.14);
2235   canv->cd(2)->SetLeftMargin(0.18);
2236   canv->cd(2)->SetRightMargin(0.01);
2237   canv->cd(2)->SetTopMargin(0.06);
2238 
2239   canv->cd(3)->SetBottomMargin(0.14);
2240   canv->cd(3)->SetLeftMargin(0.18);
2241   canv->cd(3)->SetRightMargin(0.01);
2242   canv->cd(3)->SetTopMargin(0.06);
2243 
2244   canv->cd(4)->SetBottomMargin(0.14);
2245   canv->cd(4)->SetLeftMargin(0.18);
2246   canv->cd(4)->SetRightMargin(0.01);
2247   canv->cd(4)->SetTopMargin(0.06);
2248 
2249   TH1F *dBiasTrend[4][nFiles];
2250 
2251   for (Int_t i = 0; i < nFiles; i++) {
2252     dBiasTrend[0][i] = dxyPhiMeanTrend[i];
2253     dBiasTrend[1][i] = dzPhiMeanTrend[i];
2254     dBiasTrend[2][i] = dxyEtaMeanTrend[i];
2255     dBiasTrend[3][i] = dzEtaMeanTrend[i];
2256   }
2257 
2258   Double_t absmin[4] = {999., 999., 999., 999.};
2259   Double_t absmax[4] = {-999., -999. - 999., -999.};
2260 
2261   for (Int_t k = 0; k < 4; k++) {
2262     canv->cd(k + 1);
2263 
2264     for (Int_t i = 0; i < nFiles; i++) {
2265       if (TString(canv->GetName()).Contains("BareResiduals")) {
2266         dBiasTrend[k][i]->Scale(1. / dBiasTrend[k][i]->GetSumOfWeights());
2267       }
2268 
2269       if (dBiasTrend[k][i]->GetMaximum() > absmax[k])
2270         absmax[k] = dBiasTrend[k][i]->GetMaximum();
2271       if (dBiasTrend[k][i]->GetMinimum() < absmin[k])
2272         absmin[k] = dBiasTrend[k][i]->GetMinimum();
2273     }
2274 
2275     Double_t safeDelta = (absmax[k] - absmin[k]) / 8.;
2276 
2277     // if(safeDelta<0.1*absmax[k]) safeDelta*=16;
2278 
2279     Double_t theExtreme = std::max(absmax[k], TMath::Abs(absmin[k]));
2280 
2281     for (Int_t i = 0; i < nFiles; i++) {
2282       if (i == 0) {
2283         TString theTitle = dBiasTrend[k][i]->GetName();
2284 
2285         //std::cout << theTitle<< " --->" << safeDelta << std::endl;
2286 
2287         // if the autoLimits are not set
2288         if (!setAutoLimits) {
2289           params::measurement range = getTheRangeUser(dBiasTrend[k][i], thePlotLimits);
2290           dBiasTrend[k][i]->GetYaxis()->SetRangeUser(range.first, range.second);
2291 
2292         } else {
2293           if (theTitle.Contains("width")) {
2294             if (theTitle.Contains("Norm"))
2295               safeDelta = (theTitle.Contains("ladder") == true || theTitle.Contains("modZ") == true) ? 1. : safeDelta;
2296             else
2297               safeDelta = (theTitle.Contains("ladder") == true || theTitle.Contains("modZ") == true) ? safeDelta * 10.
2298                                                                                                      : safeDelta;
2299 
2300             dBiasTrend[k][i]->GetYaxis()->SetRangeUser(0., theExtreme + (safeDelta / 2.));
2301           } else {
2302             if (theTitle.Contains("Norm")) {
2303               dBiasTrend[k][i]->GetYaxis()->SetRangeUser(std::min(-0.48, absmin[k] - (safeDelta / 2.)),
2304                                                          std::max(0.48, absmax[k] + (safeDelta / 2.)));
2305             } else if (theTitle.Contains("h_probe")) {
2306               TGaxis::SetMaxDigits(4);
2307               dBiasTrend[k][i]->GetYaxis()->SetRangeUser(0., theExtreme + (safeDelta * 2.));
2308             } else {
2309               safeDelta = (theTitle.Contains("ladder") == true || theTitle.Contains("modZ") == true) ? safeDelta * 10.
2310                                                                                                      : safeDelta;
2311 
2312               dBiasTrend[k][i]->GetYaxis()->SetRangeUser(-theExtreme - (safeDelta / 2.), theExtreme + (safeDelta / 2.));
2313             }
2314           }
2315         }
2316 
2317         if (TString(canv->GetName()).Contains("BareResiduals") && (k == 0 || k == 2)) {
2318           dBiasTrend[k][i]->GetXaxis()->SetRangeUser(-0.11, 0.11);
2319         }
2320 
2321         dBiasTrend[k][i]->Draw("e1");
2322         makeNewXAxis(dBiasTrend[k][i]);
2323         Int_t nbins = dBiasTrend[k][i]->GetNbinsX();
2324         Double_t lowedge = dBiasTrend[k][i]->GetBinLowEdge(1);
2325         Double_t highedge = dBiasTrend[k][i]->GetBinLowEdge(nbins + 1);
2326 
2327         /*
2328       TH1F* zeros = DrawZero(dBiasTrend[k][i],nbins,lowedge,highedge,1);
2329       zeros->Draw("PLsame"); 
2330     */
2331 
2332         Double_t theC = -1.;
2333 
2334         if (theTitle.Contains("width")) {
2335           if (theTitle.Contains("Norm")) {
2336             theC = 1.;
2337           } else {
2338             theC = -1.;
2339           }
2340         } else {
2341           theC = 0.;
2342         }
2343 
2344         TH1F *theConst = DrawConstant(dBiasTrend[k][i], nbins, lowedge, highedge, 1, theC);
2345         theConst->Draw("PLsame");
2346 
2347       } else {
2348         if (TString(canv->GetName()).Contains("BareResiduals") && (k == 0 || k == 2)) {
2349           dBiasTrend[k][i]->GetXaxis()->SetRangeUser(-0.11, 0.11);
2350         }
2351 
2352         dBiasTrend[k][i]->Draw("e1sames");
2353       }
2354 
2355       if (k == 0) {
2356         lego->AddEntry(dBiasTrend[k][i], LegLabels[i]);
2357       }
2358     }
2359 
2360     lego->Draw();
2361 
2362     TPad *current_pad = static_cast<TPad *>(canv->GetPad(k + 1));
2363     CMS_lumi(current_pad, 4, 33);
2364     if (theDate != "")
2365       ptDate->Draw("same");
2366   }
2367 }
2368 
2369 //*************************************************************
2370 void arrangeCanvas(TCanvas *canv,
2371                    TH1F *meanplots[100],
2372                    TH1F *widthplots[100],
2373                    Int_t nFiles,
2374                    TString LegLabels[10],
2375                    TString theDate,
2376                    bool onlyBias,
2377                    bool setAutoLimits) {
2378   //*************************************************************
2379 
2380   TPaveText *ali = new TPaveText(0.18, 0.85, 0.50, 0.93, "NDC");
2381   ali->SetFillColor(10);
2382   ali->SetTextColor(1);
2383   ali->SetTextFont(42);
2384   ali->SetMargin(0.);
2385   ali->SetLineColor(10);
2386   ali->SetShadowColor(10);
2387   TText *alitext = ali->AddText("Alignment: PCL");
2388   alitext->SetTextSize(0.04);
2389 
2390   TLegend *lego = new TLegend(0.22, 0.80, 0.78, 0.91);
2391   // in case many objects are compared
2392   if (nFiles > 3) {
2393     lego->SetNColumns(2);
2394   }
2395   // TLegend *lego = new TLegend(0.18,0.77,0.50,0.86);
2396   lego->SetFillColor(10);
2397   if (nFiles > 3) {
2398     lego->SetTextSize(0.03);
2399   } else {
2400     lego->SetTextSize(0.04);
2401   }
2402   lego->SetTextFont(42);
2403   lego->SetFillColor(10);
2404   lego->SetLineColor(10);
2405   lego->SetShadowColor(10);
2406 
2407   TPaveText *ptDate = NULL;
2408 
2409   canv->SetFillColor(10);
2410 
2411   if (!onlyBias) {
2412     ptDate = new TPaveText(0.20, 0.95, 0.50, 0.99, "blNDC");
2413   } else {
2414     ptDate = new TPaveText(0.20, 0.95, 0.50, 0.99, "blNDC");
2415   }
2416 
2417   //ptDate->SetFillColor(kYellow);
2418   ptDate->SetFillColor(10);
2419   ptDate->SetBorderSize(1);
2420   ptDate->SetLineColor(kBlue);
2421   ptDate->SetLineWidth(1);
2422   ptDate->SetTextFont(42);
2423   TText *textDate = ptDate->AddText(theDate);
2424   textDate->SetTextSize(0.04);
2425   textDate->SetTextColor(kBlue);
2426 
2427   if (!onlyBias) {
2428     canv->Divide(2, 1);
2429 
2430     canv->cd(1)->SetBottomMargin(0.14);
2431     canv->cd(1)->SetLeftMargin(0.17);
2432     canv->cd(1)->SetRightMargin(0.02);
2433     canv->cd(1)->SetTopMargin(0.06);
2434 
2435     canv->cd(2)->SetBottomMargin(0.14);
2436     canv->cd(2)->SetLeftMargin(0.17);
2437     canv->cd(2)->SetRightMargin(0.02);
2438     canv->cd(2)->SetTopMargin(0.06);
2439     canv->cd(1);
2440 
2441   } else {
2442     canv->cd()->SetBottomMargin(0.14);
2443     canv->cd()->SetLeftMargin(0.17);
2444     canv->cd()->SetRightMargin(0.02);
2445     canv->cd()->SetTopMargin(0.06);
2446     canv->cd();
2447   }
2448 
2449   Double_t absmin(999.);
2450   Double_t absmax(-999.);
2451 
2452   for (Int_t i = 0; i < nFiles; i++) {
2453     if (meanplots[i]->GetMaximum() > absmax)
2454       absmax = meanplots[i]->GetMaximum();
2455     if (meanplots[i]->GetMinimum() < absmin)
2456       absmin = meanplots[i]->GetMinimum();
2457   }
2458 
2459   Double_t safeDelta = (absmax - absmin) / 2.;
2460   Double_t theExtreme = std::max(absmax, TMath::Abs(absmin));
2461 
2462   for (Int_t i = 0; i < nFiles; i++) {
2463     if (i == 0) {
2464       // if the autoLimits are not set
2465       if (!setAutoLimits) {
2466         params::measurement range = getTheRangeUser(meanplots[i], thePlotLimits);
2467         meanplots[i]->GetYaxis()->SetRangeUser(range.first, range.second);
2468 
2469       } else {
2470         TString theTitle = meanplots[i]->GetName();
2471         if (theTitle.Contains("Norm")) {
2472           meanplots[i]->GetYaxis()->SetRangeUser(std::min(-0.48, absmin - safeDelta),
2473                                                  std::max(0.48, absmax + safeDelta));
2474         } else {
2475           if (!onlyBias) {
2476             meanplots[i]->GetYaxis()->SetRangeUser(absmin - safeDelta, absmax + safeDelta);
2477           } else {
2478             meanplots[i]->GetYaxis()->SetRangeUser(-theExtreme - (TMath::Abs(absmin) / 10.),
2479                                                    theExtreme + (TMath::Abs(absmax / 10.)));
2480           }
2481         }
2482       }
2483 
2484       meanplots[i]->Draw("e1");
2485       if (TString(meanplots[i]->GetName()).Contains("pT")) {
2486         //meanplots[i]->Draw("HIST][same");
2487         gPad->SetLogx();
2488         gPad->SetGridx();
2489         gPad->SetGridy();
2490       } else {
2491         makeNewXAxis(meanplots[i]);
2492       }
2493 
2494       if (onlyBias) {
2495         canv->cd();
2496         Int_t nbins = meanplots[i]->GetNbinsX();
2497         Double_t lowedge = meanplots[i]->GetBinLowEdge(1);
2498         Double_t highedge = meanplots[i]->GetBinLowEdge(nbins + 1);
2499 
2500         TH1F *hzero = DrawZero(meanplots[i], nbins, lowedge, highedge, 2);
2501         hzero->Draw("PLsame");
2502       }
2503     } else
2504       meanplots[i]->Draw("e1sames");
2505     //if(TString(meanplots[i]->GetName()).Contains("pT")){
2506     //  meanplots[i]->Draw("HIST][same");
2507     // }
2508 
2509     lego->AddEntry(meanplots[i], LegLabels[i]);
2510   }
2511 
2512   lego->Draw();
2513 
2514   //ali->Draw("same");
2515   //ptDate->Draw("same");
2516 
2517   TPad *current_pad;
2518   if (!onlyBias) {
2519     current_pad = static_cast<TPad *>(canv->GetPad(1));
2520   } else {
2521     current_pad = static_cast<TPad *>(canv->GetPad(0));
2522   }
2523 
2524   CMS_lumi(current_pad, 4, 33);
2525   if (theDate != "")
2526     ptDate->Draw("same");
2527 
2528   if (!onlyBias) {
2529     canv->cd(2);
2530     Double_t absmax2(-999.);
2531 
2532     for (Int_t i = 0; i < nFiles; i++) {
2533       if (widthplots[i]->GetMaximum() > absmax2)
2534         absmax2 = widthplots[i]->GetMaximum();
2535     }
2536 
2537     Double_t safeDelta2 = absmax2 / 3.;
2538 
2539     for (Int_t i = 0; i < nFiles; i++) {
2540       widthplots[i]->GetXaxis()->SetLabelOffset(999);
2541       widthplots[i]->GetXaxis()->SetTickLength(0);
2542 
2543       if (i == 0) {
2544         if (!setAutoLimits) {
2545           params::measurement range = getTheRangeUser(widthplots[i], thePlotLimits);
2546           widthplots[i]->GetYaxis()->SetRangeUser(range.first, range.second);
2547         } else {
2548           widthplots[i]->SetMinimum(0.5);
2549           widthplots[i]->SetMaximum(absmax2 + safeDelta2);
2550         }
2551 
2552         widthplots[i]->Draw("e1");
2553         if (TString(widthplots[i]->GetName()).Contains("pT")) {
2554           //widthplots[i]->Draw("HIST][same");
2555           gPad->SetGridx();
2556           gPad->SetGridy();
2557         }
2558         makeNewXAxis(widthplots[i]);
2559       } else {
2560         widthplots[i]->Draw("e1sames");
2561         if (TString(widthplots[i]->GetName()).Contains("pT")) {
2562           //widthplots[i]->Draw("HIST][same");
2563         }
2564       }
2565     }
2566 
2567     lego->Draw();
2568 
2569     TPad *current_pad2 = static_cast<TPad *>(canv->GetPad(2));
2570     CMS_lumi(current_pad2, 4, 33);
2571     if (theDate != "")
2572       ptDate->Draw("same");
2573   }
2574 }
2575 
2576 //*************************************************************
2577 void arrangeCanvas2D(
2578     TCanvas *canv, TH2F *meanmaps[100], TH2F *widthmaps[100], Int_t nFiles, TString LegLabels[10], TString theDate)
2579 //*************************************************************
2580 {
2581   TLegend *lego = new TLegend(0.18, 0.75, 0.58, 0.92);
2582   lego->SetFillColor(10);
2583   lego->SetTextSize(0.05);
2584   lego->SetTextFont(42);
2585   lego->SetFillColor(10);
2586   lego->SetLineColor(10);
2587   lego->SetShadowColor(10);
2588 
2589   TPaveText *pt[nFiles];
2590   TPaveText *pt2[nFiles];
2591   TPaveText *pt3[nFiles];
2592 
2593   for (Int_t i = 0; i < nFiles; i++) {
2594     pt[i] = new TPaveText(0.13, 0.95, 0.191, 0.975, "NDC");
2595     //pt[i] = new TPaveText(gPad->GetUxmin(),gPad->GetUymax()+0.3,gPad->GetUxmin()+0.6,gPad->GetUymax()+0.3,"NDC");
2596     //std::cout<<"gPad->GetUymax():"<<gPad->GetUymax()<<std::endl;
2597     //pt[i] = new TPaveText(gPad->GetLeftMargin(),0.95,gPad->GetLeftMargin()+0.3,0.98,"NDC");
2598     pt[i]->SetFillColor(10);
2599     pt[i]->SetTextColor(1);
2600     pt[i]->SetTextFont(61);
2601     pt[i]->SetTextAlign(22);
2602     TText *text1 = pt[i]->AddText("CMS");  // preliminary 2015 p-p data, #sqrt{s}=8 TeV "+LegLabels[i]);
2603     text1->SetTextSize(0.05);
2604     //delete text1;
2605 
2606     //float extraOverCmsTextSize  = 0.76;
2607 
2608     pt2[i] = new TPaveText(0.21, 0.95, 0.25, 0.975, "NDC");
2609     pt2[i]->SetFillColor(10);
2610     pt2[i]->SetTextColor(1);
2611     //pt[i]->SetTextSize(0.05);
2612     pt2[i]->SetTextFont(52);
2613     pt2[i]->SetTextAlign(12);
2614     // TText *text2 = pt2->AddText("run: "+theDate);
2615     TText *text2 = pt2[i]->AddText("INTERNAL");
2616     text2->SetTextSize(0.06 * extraOverCmsTextSize);
2617 
2618     pt3[i] = new TPaveText(0.55, 0.955, 0.95, 0.98, "NDC");
2619     pt3[i]->SetFillColor(10);
2620     pt3[i]->SetTextColor(kBlue);
2621     pt3[i]->SetTextFont(61);
2622     pt3[i]->SetTextAlign(22);
2623     // TText *text2 = pt2->AddText("run: "+theDate);
2624     TText *text3 = pt3[i]->AddText(LegLabels[i]);
2625     text3->SetTextSize(0.05);
2626   }
2627 
2628   canv->SetFillColor(10);
2629   canv->Divide(2, nFiles);
2630 
2631   Double_t absmin(999.);
2632   Double_t absmax(-999.);
2633   Double_t maxwidth(-999.);
2634 
2635   for (Int_t i = 0; i < nFiles; i++) {
2636     if (widthmaps[i]->GetMaximum() > maxwidth)
2637       maxwidth = widthmaps[i]->GetMaximum();
2638     if (meanmaps[i]->GetMaximum() > absmax)
2639       absmax = meanmaps[i]->GetMaximum();
2640     if (meanmaps[i]->GetMinimum() < absmin)
2641       absmin = meanmaps[i]->GetMinimum();
2642   }
2643 
2644   /*
2645   const Int_t nLevels = 255;
2646   Double_t levels[nLevels];
2647 
2648   for(int i = 0; i < nLevels; i++) {
2649     levels[i] = absmin + (absmax - absmin) / (nLevels - 1) * (i);
2650   }
2651   */
2652 
2653   //Double_t theExtreme = std::min(std::abs(absmin),std::abs(absmax));
2654 
2655   for (Int_t i = 0; i < nFiles; i++) {
2656     canv->cd(2 * i + 1)->SetBottomMargin(0.13);
2657     canv->cd(2 * i + 1)->SetLeftMargin(0.12);
2658     canv->cd(2 * i + 1)->SetRightMargin(0.19);
2659     canv->cd(2 * i + 1)->SetTopMargin(0.08);
2660 
2661     //meanmaps[i]->SetContour((sizeof(levels)/sizeof(Double_t)), levels);
2662     meanmaps[i]->GetZaxis()->SetRangeUser(absmin, absmax);
2663     //meanmaps[i]->GetZaxis()->SetRangeUser(-theExtreme,theExtreme);
2664     meanmaps[i]->Draw("colz1");
2665 
2666     //TH2F* cloned = (TH2F*)meanmaps[i]->DrawClone("colz");// draw "axes", "contents", "statistics box
2667     //makeNewPairOfAxes(cloned);
2668     //meanmaps[i]->GetZaxis()->SetRangeUser(absmin, absmax); // ... set the range ...
2669     //meanmaps[i]->Draw("colzsame"); // draw the "color palette"
2670 
2671     makeNewPairOfAxes(meanmaps[i]);
2672 
2673     pt[i]->Draw("same");
2674     pt2[i]->Draw("same");
2675     pt3[i]->Draw("same");
2676 
2677     canv->cd(2 * (i + 1))->SetBottomMargin(0.13);
2678     canv->cd(2 * (i + 1))->SetLeftMargin(0.12);
2679     canv->cd(2 * (i + 1))->SetRightMargin(0.19);
2680     canv->cd(2 * (i + 1))->SetTopMargin(0.08);
2681 
2682     widthmaps[i]->Draw("colz1");
2683     makeNewPairOfAxes(widthmaps[i]);
2684 
2685     widthmaps[i]->GetZaxis()->SetRangeUser(0., maxwidth);
2686 
2687     pt[i]->Draw("same");
2688     pt2[i]->Draw("same");
2689     pt3[i]->Draw("same");
2690   }
2691 }
2692 
2693 //*************************************************************
2694 void arrangeFitCanvas(TCanvas *canv, TH1F *meanplots[100], Int_t nFiles, TString LegLabels[10], TString theDate)
2695 //*************************************************************
2696 {
2697   canv->SetBottomMargin(0.14);
2698   canv->SetLeftMargin(0.1);
2699   canv->SetRightMargin(0.02);
2700   canv->SetTopMargin(0.08);
2701 
2702   TLegend *lego = new TLegend(0.12, 0.80, 0.82, 0.89);
2703   lego->SetFillColor(10);
2704   lego->SetTextSize(0.035);
2705   lego->SetTextFont(42);
2706   lego->SetFillColor(10);
2707   lego->SetLineColor(10);
2708   if (nFiles > 3) {
2709     lego->SetNColumns(2);
2710   }
2711   lego->SetShadowColor(10);
2712 
2713   TPaveText *ptDate = new TPaveText(0.12, 0.95, 0.50, 0.99, "blNDC");
2714   //ptDate->SetFillColor(kYellow);
2715   ptDate->SetFillColor(10);
2716   ptDate->SetBorderSize(1);
2717   ptDate->SetLineColor(kBlue);
2718   ptDate->SetLineWidth(1);
2719   ptDate->SetTextFont(32);
2720   TText *textDate = ptDate->AddText(theDate);
2721   textDate->SetTextSize(0.04);
2722   textDate->SetTextColor(kBlue);
2723 
2724   TF1 *fleft[nFiles];
2725   TF1 *fright[nFiles];
2726   TF1 *fall[nFiles];
2727 
2728   TF1 *FitDzUp[nFiles];
2729   TF1 *FitDzDown[nFiles];
2730 
2731   for (Int_t j = 0; j < nFiles; j++) {
2732     Double_t deltaZ(0);
2733     Double_t sigmadeltaZ(-1);
2734 
2735     TCanvas *theNewCanvas2 = new TCanvas("NewCanvas2", "Fitting Canvas 2", 800, 600);
2736     theNewCanvas2->Divide(2, 1);
2737 
2738     TH1F *hnewUp = (TH1F *)meanplots[j]->Clone("hnewUp_dz_phi");
2739     TH1F *hnewDown = (TH1F *)meanplots[j]->Clone("hnewDown_dz_phi");
2740 
2741     fleft[j] = new TF1(Form("fleft_%i", j), fULine, _boundMin, _boundSx, 1);
2742     fright[j] = new TF1(Form("fright_%i", j), fULine, _boundDx, _boundMax, 1);
2743     fall[j] = new TF1(Form("fall_%i", j), fDLine, _boundSx, _boundDx, 1);
2744 
2745     FitULine(hnewUp);
2746     FitDzUp[j] = (TF1 *)hnewUp->GetListOfFunctions()->FindObject("lineUp");
2747     if (FitDzUp[j]) {
2748       fleft[j]->SetParameters(FitDzUp[j]->GetParameters());
2749       fleft[j]->SetParErrors(FitDzUp[j]->GetParErrors());
2750       hnewUp->GetListOfFunctions()->Add(fleft[j]);
2751       fright[j]->SetParameters(FitDzUp[j]->GetParameters());
2752       fright[j]->SetParErrors(FitDzUp[j]->GetParErrors());
2753       hnewUp->GetListOfFunctions()->Add(fright[j]);
2754       FitDzUp[j]->Delete();
2755 
2756       theNewCanvas2->cd(1);
2757       MakeNiceTF1Style(fright[j], meanplots[j]->GetLineColor());
2758       MakeNiceTF1Style(fleft[j], meanplots[j]->GetLineColor());
2759       fright[j]->Draw("same");
2760       fleft[j]->Draw("same");
2761     }
2762 
2763     FitDLine(hnewDown);
2764     FitDzDown[j] = (TF1 *)hnewDown->GetListOfFunctions()->FindObject("lineDown");
2765 
2766     if (FitDzDown[j]) {
2767       fall[j]->SetParameters(FitDzDown[j]->GetParameters());
2768       fall[j]->SetParErrors(FitDzDown[j]->GetParErrors());
2769       hnewDown->GetListOfFunctions()->Add(fall[j]);
2770       FitDzDown[j]->Delete();
2771       theNewCanvas2->cd(2);
2772       MakeNiceTF1Style(fall[j], meanplots[j]->GetLineColor());
2773       fall[j]->Draw("same");
2774       canv->cd();
2775       hnewUp->GetYaxis()->SetTitleOffset(0.7);
2776       if (j == 0) {
2777         hnewUp->Draw();
2778         makeNewXAxis(hnewUp);
2779       } else {
2780         hnewUp->Draw("same");
2781         makeNewXAxis(hnewUp);
2782       }
2783       fright[j]->Draw("same");
2784       fleft[j]->Draw("same");
2785       fall[j]->Draw("same");
2786     }
2787 
2788     if (j == nFiles - 1) {
2789       theNewCanvas2->Close();
2790     }
2791 
2792     deltaZ = (fright[j]->GetParameter(0) - fall[j]->GetParameter(0)) / 2;
2793     sigmadeltaZ = 0.5 * TMath::Sqrt(fright[j]->GetParError(0) * fright[j]->GetParError(0) +
2794                                     fall[j]->GetParError(0) * fall[j]->GetParError(0));
2795     TString MYOUT = Form(" : #Delta z = %.f #pm %.f #mum", deltaZ, sigmadeltaZ);
2796 
2797     lego->AddEntry(meanplots[j], LegLabels[j] + MYOUT);
2798 
2799     if (j == nFiles - 1) {
2800       outfile << deltaZ << "|" << sigmadeltaZ << std::endl;
2801     }
2802 
2803     delete theNewCanvas2;
2804   }
2805 
2806   //TkAlStyle::drawStandardTitle(Coll0T15);
2807   lego->Draw("same");
2808   CMS_lumi(canv, 4, 33);
2809   if (theDate != "")
2810     ptDate->Draw("same");
2811   //pt->Draw("same");
2812 }
2813 
2814 //*************************************************************
2815 std::pair<params::measurement, params::measurement> fitStudentTResiduals(TH1 *hist)
2816 //*************************************************************
2817 {
2818   hist->SetMarkerStyle(21);
2819   hist->SetMarkerSize(0.8);
2820   hist->SetStats(1);
2821 
2822   double dx = hist->GetBinWidth(1);
2823   double nmax = hist->GetBinContent(hist->GetMaximumBin());
2824   double xmax = hist->GetBinCenter(hist->GetMaximumBin());
2825   double nn = 7 * nmax;
2826 
2827   int nb = hist->GetNbinsX();
2828   double n1 = hist->GetBinContent(1);
2829   double n9 = hist->GetBinContent(nb);
2830   double bg = 0.5 * (n1 + n9);
2831 
2832   double x1 = hist->GetBinCenter(1);
2833   double x9 = hist->GetBinCenter(nb);
2834 
2835   // create a TF1 with the range from x1 to x9 and 5 parameters
2836 
2837   TF1 *tp0Fcn = new TF1("tmp", tp0Fit, x1, x9, 5);
2838 
2839   tp0Fcn->SetParName(0, "mean");
2840   tp0Fcn->SetParName(1, "sigma");
2841   tp0Fcn->SetParName(2, "nu");
2842   tp0Fcn->SetParName(3, "area");
2843   tp0Fcn->SetParName(4, "BG");
2844 
2845   tp0Fcn->SetNpx(500);
2846   tp0Fcn->SetLineWidth(2);
2847   //tp0Fcn->SetLineColor(kMagenta);
2848   //tp0Fcn->SetLineColor(kGreen);
2849   tp0Fcn->SetLineColor(kRed);
2850 
2851   // set start values for some parameters:
2852 
2853   tp0Fcn->SetParameter(0, xmax);    // peak position
2854   tp0Fcn->SetParameter(1, 4 * dx);  // width
2855   tp0Fcn->SetParameter(2, 2.2);     // nu
2856   tp0Fcn->SetParameter(3, nn);      // N
2857   tp0Fcn->SetParameter(4, bg);
2858 
2859   hist->Fit("tmp", "R", "ep");
2860   // h->Fit("tmp","V+","ep");
2861 
2862   hist->Draw("histepsame");  // data again on top
2863 
2864   float res_mean = tp0Fcn->GetParameter(0);
2865   float res_width = tp0Fcn->GetParameter(1);
2866 
2867   float res_mean_err = tp0Fcn->GetParError(0);
2868   float res_width_err = tp0Fcn->GetParError(1);
2869 
2870   params::measurement resultM;
2871   params::measurement resultW;
2872 
2873   resultM = std::make_pair(res_mean, res_mean_err);
2874   resultW = std::make_pair(res_width, res_width_err);
2875 
2876   std::pair<params::measurement, params::measurement> result;
2877 
2878   result = std::make_pair(resultM, resultW);
2879   return result;
2880 }
2881 
2882 //*************************************************************
2883 Double_t tp0Fit(Double_t *x, Double_t *par5)
2884 //*************************************************************
2885 {
2886   static int nn = 0;
2887   nn++;
2888   static double dx = 0.1;
2889   static double b1 = 0;
2890   if (nn == 1)
2891     b1 = x[0];
2892   if (nn == 2)
2893     dx = x[0] - b1;
2894   //
2895   //--  Mean and width:
2896   //
2897   double xm = par5[0];
2898   double t = (x[0] - xm) / par5[1];
2899   double tt = t * t;
2900   //
2901   //--  exponent:
2902   //
2903   double rn = par5[2];
2904   double xn = 0.5 * (rn + 1.0);
2905   //
2906   //--  Normalization needs Gamma function:
2907   //
2908   double pk = 0.0;
2909 
2910   if (rn > 0.0) {
2911     double pi = 3.14159265358979323846;
2912     double aa = dx / par5[1] / sqrt(rn * pi) * TMath::Gamma(xn) / TMath::Gamma(0.5 * rn);
2913 
2914     pk = par5[3] * aa * exp(-xn * log(1.0 + tt / rn));
2915   }
2916 
2917   return pk + par5[4];
2918 }
2919 
2920 //*************************************************************
2921 params::measurement getMedian(TH1F *histo)
2922 //*************************************************************
2923 {
2924   Double_t median = 999;
2925   int nbins = histo->GetNbinsX();
2926 
2927   //extract median from histogram
2928   double *x = new double[nbins];
2929   double *y = new double[nbins];
2930   for (int j = 0; j < nbins; j++) {
2931     x[j] = histo->GetBinCenter(j + 1);
2932     y[j] = histo->GetBinContent(j + 1);
2933   }
2934   median = TMath::Median(nbins, x, y);
2935 
2936   delete[] x;
2937   x = 0;
2938   delete[] y;
2939   y = 0;
2940 
2941   params::measurement result;
2942   result = std::make_pair(median, median / TMath::Sqrt(histo->GetEntries()));
2943 
2944   return result;
2945 }
2946 
2947 //*************************************************************
2948 params::measurement getMAD(TH1F *histo)
2949 //*************************************************************
2950 {
2951   int nbins = histo->GetNbinsX();
2952   Double_t median = getMedian(histo).first;
2953   Double_t x_lastBin = histo->GetBinLowEdge(nbins + 1);
2954   const char *HistoName = histo->GetName();
2955   TString Finalname = Form("resMed%s", HistoName);
2956   TH1F *newHisto = new TH1F(Finalname, Finalname, nbins, 0., x_lastBin);
2957   Double_t *residuals = new Double_t[nbins];
2958   Double_t *weights = new Double_t[nbins];
2959 
2960   for (int j = 0; j < nbins; j++) {
2961     residuals[j] = TMath::Abs(median - histo->GetBinCenter(j + 1));
2962     weights[j] = histo->GetBinContent(j + 1);
2963     newHisto->Fill(residuals[j], weights[j]);
2964   }
2965 
2966   Double_t theMAD = (getMedian(newHisto).first) * 1.4826;
2967   newHisto->Delete("");
2968 
2969   params::measurement result;
2970   result = std::make_pair(theMAD, theMAD / histo->GetEntries());
2971 
2972   return result;
2973 }
2974 
2975 //*************************************************************
2976 std::pair<params::measurement, params::measurement> fitResiduals(TH1 *hist, bool singleTime)
2977 //*************************************************************
2978 {
2979   assert(hist != nullptr);
2980 
2981   if (hist->GetEntries() < 10) {
2982     // std::cout<<"hist name: "<<hist->GetName() << std::endl;
2983     return std::make_pair(std::make_pair(0., 0.), std::make_pair(0., 0.));
2984   }
2985 
2986   float maxHist = hist->GetXaxis()->GetXmax();
2987   float minHist = hist->GetXaxis()->GetXmin();
2988   float mean = hist->GetMean();
2989   float sigma = hist->GetRMS();
2990 
2991   if (TMath::IsNaN(mean) || TMath::IsNaN(sigma)) {
2992     mean = 0;
2993     //sigma= - hist->GetXaxis()->GetBinLowEdge(1) + hist->GetXaxis()->GetBinLowEdge(hist->GetNbinsX()+1);
2994     sigma = -minHist + maxHist;
2995     std::cout << "FitPVResiduals::fitResiduals(): histogram" << hist->GetName() << " mean or sigma are NaN!!"
2996               << std::endl;
2997   }
2998 
2999   TF1 func("tmp", "gaus", mean - 2. * sigma, mean + 2. * sigma);
3000   if (0 == hist->Fit(&func, "QNR")) {  // N: do not blow up file by storing fit!
3001     mean = func.GetParameter(1);
3002     sigma = func.GetParameter(2);
3003 
3004     if (!singleTime) {
3005       // second fit: three sigma of first fit around mean of first fit
3006       func.SetRange(std::max(mean - 2 * sigma, minHist), std::min(mean + 2 * sigma, maxHist));
3007       // I: integral gives more correct results if binning is too wide
3008       // L: Likelihood can treat empty bins correctly (if hist not weighted...)
3009       if (0 == hist->Fit(&func, "Q0LR")) {
3010         if (hist->GetFunction(func.GetName())) {  // Take care that it is later on drawn:
3011           hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
3012         }
3013       }
3014     }
3015   }
3016 
3017   /*
3018   float res_mean  = func.GetParameter(1);
3019   float res_width = func.GetParameter(2);
3020   
3021   float res_mean_err  = func.GetParError(1);
3022   float res_width_err = func.GetParError(2);
3023 
3024   params::measurement resultM;
3025   params::measurement resultW;
3026 
3027   resultM = std::make_pair(res_mean,res_mean_err);
3028   resultW = std::make_pair(res_width,res_width_err);
3029 
3030   std::pair<params::measurement, params::measurement  > result;
3031   
3032   result = std::make_pair(resultM,resultW);
3033   */
3034   return std::make_pair(std::make_pair(func.GetParameter(1), func.GetParError(1)),
3035                         std::make_pair(func.GetParameter(2), func.GetParError(2)));
3036 }
3037 
3038 //*************************************************************
3039 Double_t DoubleSidedCB(double *x, double *par) {
3040   //*************************************************************
3041 
3042   double m = x[0];
3043   double m0 = par[0];
3044   double sigma = par[1];
3045   double alphaL = par[2];
3046   double alphaR = par[3];
3047   double nL = par[4];
3048   double nR = par[5];
3049   double N = par[6];
3050 
3051   Double_t arg = m - m0;
3052 
3053   if (arg < 0.0) {
3054     Double_t t = (m - m0) / sigma;               //t < 0
3055     Double_t absAlpha = fabs((Double_t)alphaL);  //slightly redundant since alpha > 0 anyway, but never mind
3056     if (t >= -absAlpha) {                        //-absAlpha <= t < 0
3057       return N * exp(-0.5 * t * t);
3058     } else {
3059       Double_t a = TMath::Power(nL / absAlpha, nL) * exp(-0.5 * absAlpha * absAlpha);
3060       Double_t b = nL / absAlpha - absAlpha;
3061       return N * (a / TMath::Power(b - t, nL));  //b - t
3062     }
3063   } else {
3064     Double_t t = (m - m0) / sigma;  //t > 0
3065     Double_t absAlpha = fabs((Double_t)alphaR);
3066     if (t <= absAlpha) {  //0 <= t <= absAlpha
3067       return N * exp(-0.5 * t * t);
3068     } else {
3069       Double_t a = TMath::Power(nR / absAlpha, nR) * exp(-0.5 * absAlpha * absAlpha);
3070       Double_t b = nR / absAlpha - absAlpha;
3071       return N * (a / TMath::Power(b + t, nR));  //b + t
3072     }
3073   }
3074 }
3075 
3076 //*************************************************************
3077 std::pair<params::measurement, params::measurement> fitResidualsCB(TH1 *hist)
3078 //*************************************************************
3079 {
3080   //hist->Rebin(2);
3081 
3082   float mean = hist->GetMean();
3083   float sigma = hist->GetRMS();
3084   //int   nbinsX   = hist->GetNbinsX();
3085   float nentries = hist->GetEntries();
3086   float meanerr = sigma / TMath::Sqrt(nentries);
3087   float sigmaerr = TMath::Sqrt(sigma * sigma * TMath::Sqrt(2 / nentries));
3088 
3089   float lowBound = hist->GetXaxis()->GetBinLowEdge(1);
3090   float highBound = hist->GetXaxis()->GetBinLowEdge(hist->GetNbinsX() + 1);
3091 
3092   if (TMath::IsNaN(mean) || TMath::IsNaN(sigma)) {
3093     mean = 0;
3094     sigma = -lowBound + highBound;
3095   }
3096 
3097   TF1 func("tmp", "gaus", mean - 1. * sigma, mean + 1. * sigma);
3098   if (0 == hist->Fit(&func, "QNR")) {  // N: do not blow up file by storing fit!
3099     mean = func.GetParameter(1);
3100     sigma = func.GetParameter(2);
3101   }
3102 
3103   // first round
3104   TF1 *doubleCB = new TF1("myDoubleCB", DoubleSidedCB, lowBound, highBound, 7);
3105   doubleCB->SetParameters(mean, sigma, 1.5, 1.5, 2.5, 2.5, 100);
3106   doubleCB->SetParLimits(0, mean - meanerr, mean + meanerr);
3107   doubleCB->SetParLimits(1, 0., sigma + 2 * sigmaerr);
3108   doubleCB->SetParLimits(2, 0., 30.);
3109   doubleCB->SetParLimits(3, 0., 30.);
3110   doubleCB->SetParLimits(4, 0., 50.);
3111   doubleCB->SetParLimits(5, 0., 50.);
3112   doubleCB->SetParLimits(6, 0., 100 * nentries);
3113 
3114   doubleCB->SetParNames("#mu", "#sigma", "#alpha_{L}", "#alpha_{R}", "n_{L}", "n_{R}", "N");
3115   doubleCB->SetLineColor(kRed);
3116   doubleCB->SetNpx(1000);
3117   // doubleCB->SetRange(0.8*lowBound,0.8*highBound);
3118 
3119   hist->Fit(doubleCB, "QM");
3120 
3121   // second round
3122 
3123   float p0 = doubleCB->GetParameter(0);
3124   float p1 = doubleCB->GetParameter(1);
3125   float p2 = doubleCB->GetParameter(2);
3126   float p3 = doubleCB->GetParameter(3);
3127   float p4 = doubleCB->GetParameter(4);
3128   float p5 = doubleCB->GetParameter(5);
3129   float p6 = doubleCB->GetParameter(6);
3130 
3131   float p0err = doubleCB->GetParError(0);
3132   float p1err = doubleCB->GetParError(1);
3133   float p2err = doubleCB->GetParError(2);
3134   float p3err = doubleCB->GetParError(3);
3135   float p4err = doubleCB->GetParError(4);
3136   float p5err = doubleCB->GetParError(5);
3137   float p6err = doubleCB->GetParError(6);
3138 
3139   if ((doubleCB->GetChisquare() / doubleCB->GetNDF()) > 5) {
3140     std::cout << "------------------------" << std::endl;
3141     std::cout << "chi2 1st:" << doubleCB->GetChisquare() << std::endl;
3142 
3143     //std::cout<<"p0: "<<p0<<"+/-"<<p0err<<std::endl;
3144     //std::cout<<"p1: "<<p1<<"+/-"<<p1err<<std::endl;
3145     //std::cout<<"p2: "<<p2<<"+/-"<<p2err<<std::endl;
3146     //std::cout<<"p3: "<<p3<<"+/-"<<p3err<<std::endl;
3147     //std::cout<<"p4: "<<p4<<"+/-"<<p4err<<std::endl;
3148     //std::cout<<"p5: "<<p5<<"+/-"<<p5err<<std::endl;
3149     //std::cout<<"p6: "<<p6<<"+/-"<<p6err<<std::endl;
3150 
3151     doubleCB->SetParameters(p0, p1, 3, 3, 6, 6, p6);
3152     doubleCB->SetParLimits(0, p0 - 2 * p0err, p0 + 2 * p0err);
3153     doubleCB->SetParLimits(1, p1 - 2 * p1err, p0 + 2 * p1err);
3154     doubleCB->SetParLimits(2, p2 - 2 * p2err, p0 + 2 * p2err);
3155     doubleCB->SetParLimits(3, p3 - 2 * p3err, p0 + 2 * p3err);
3156     doubleCB->SetParLimits(4, p4 - 2 * p4err, p0 + 2 * p4err);
3157     doubleCB->SetParLimits(5, p5 - 2 * p5err, p0 + 2 * p5err);
3158     doubleCB->SetParLimits(6, p6 - 2 * p6err, p0 + 2 * p6err);
3159 
3160     hist->Fit(doubleCB, "MQ");
3161 
3162     //gMinuit->Command("SCAn 1");
3163     //TGraph *gr = (TGraph*)gMinuit->GetPlot();
3164     //gr->SetMarkerStyle(21);
3165     //gr->Draw("alp");
3166 
3167     std::cout << "chi2 2nd:" << doubleCB->GetChisquare() << std::endl;
3168   }
3169 
3170   float res_mean = doubleCB->GetParameter(0);
3171   float res_width = doubleCB->GetParameter(1);
3172 
3173   float res_mean_err = doubleCB->GetParError(0);
3174   float res_width_err = doubleCB->GetParError(1);
3175 
3176   params::measurement resultM;
3177   params::measurement resultW;
3178 
3179   resultM = std::make_pair(res_mean, res_mean_err);
3180   resultW = std::make_pair(res_width, res_width_err);
3181 
3182   std::pair<params::measurement, params::measurement> result;
3183 
3184   result = std::make_pair(resultM, resultW);
3185   return result;
3186 }
3187 
3188 //*************************************************************
3189 void FillTrendPlot(TH1F *trendPlot, TH1F *residualsPlot[100], params::estimator fitPar_, TString var_, Int_t myBins_)
3190 //*************************************************************
3191 {
3192   //std::cout<<"trendPlot name: "<<trendPlot->GetName()<<std::endl;
3193 
3194   // float phiInterval = (360.)/myBins_;
3195   float phiInterval = (2 * TMath::Pi() / myBins_);
3196   float etaInterval = 5. / myBins_;
3197 
3198   for (int i = 0; i < myBins_; ++i) {
3199     //int binn = i+1;
3200 
3201     char phipositionString[129];
3202     // float phiposition = (-180+i*phiInterval)+(phiInterval/2);
3203     float phiposition = (-TMath::Pi() + i * phiInterval) + (phiInterval / 2);
3204     sprintf(phipositionString, "%.1f", phiposition);
3205 
3206     char etapositionString[129];
3207     float etaposition = (-etaRange + i * etaInterval) + (etaInterval / 2);
3208     sprintf(etapositionString, "%.1f", etaposition);
3209 
3210     std::pair<params::measurement, params::measurement> myFit =
3211         std::make_pair(std::make_pair(0., 0.), std::make_pair(0., 0.));
3212 
3213     if (((TString)trendPlot->GetName()).Contains("Norm")) {
3214       myFit = fitResiduals(residualsPlot[i]);
3215     } else {
3216       // myFit = fitStudentTResiduals(residualsPlot[i]);
3217       myFit = fitResiduals(residualsPlot[i]);
3218 
3219       /*
3220     if(TString(residualsPlot[i]->GetName()).Contains("dx") ||
3221     TString(residualsPlot[i]->GetName()).Contains("dy")) {
3222     std::cout<<residualsPlot[i]->GetName() << " " << myFit.first.first << "+/- " << myFit.first.second  << std::endl; 
3223     }
3224       */
3225     }
3226 
3227     switch (fitPar_) {
3228       case params::MEAN: {
3229         float mean_ = myFit.first.first;
3230         float meanErr_ = myFit.first.second;
3231         trendPlot->SetBinContent(i + 1, mean_);
3232         trendPlot->SetBinError(i + 1, meanErr_);
3233         break;
3234       }
3235       case params::WIDTH: {
3236         float width_ = myFit.second.first;
3237         float widthErr_ = myFit.second.second;
3238         trendPlot->SetBinContent(i + 1, width_);
3239         trendPlot->SetBinError(i + 1, widthErr_);
3240         break;
3241       }
3242       case params::MEDIAN: {
3243         float median_ = getMedian(residualsPlot[i]).first;
3244         float medianErr_ = getMedian(residualsPlot[i]).second;
3245         trendPlot->SetBinContent(i + 1, median_);
3246         trendPlot->SetBinError(i + 1, medianErr_);
3247         break;
3248       }
3249       case params::MAD: {
3250         float mad_ = getMAD(residualsPlot[i]).first;
3251         float madErr_ = getMAD(residualsPlot[i]).second;
3252         trendPlot->SetBinContent(i + 1, mad_);
3253         trendPlot->SetBinError(i + 1, madErr_);
3254         break;
3255       }
3256       default:
3257         std::cout << "PrimaryVertexValidation::FillTrendPlot() " << fitPar_ << " unknown estimator!" << std::endl;
3258         break;
3259     }
3260   }
3261 
3262   //trendPlot->GetXaxis()->LabelsOption("h");
3263 
3264   if (fitPar_ == params::MEAN || fitPar_ == params::MEDIAN) {
3265     TString res;
3266     if (TString(residualsPlot[0]->GetName()).Contains("dxy"))
3267       res = "dxy";
3268     else if (TString(residualsPlot[0]->GetName()).Contains("dx"))
3269       res = "dx";
3270     else if (TString(residualsPlot[0]->GetName()).Contains("dy"))
3271       res = "dy";
3272     else if (TString(residualsPlot[0]->GetName()).Contains("dz"))
3273       res = "dz";
3274     else if (TString(residualsPlot[0]->GetName()).Contains("IP2D"))
3275       res = "IP2D";
3276     else if (TString(residualsPlot[0]->GetName()).Contains("resz"))
3277       res = "resz";
3278 
3279     TCanvas *fitOutput = new TCanvas(Form("fitOutput_%s_%s_%s", res.Data(), var_.Data(), trendPlot->GetName()),
3280                                      Form("fitOutput_%s_%s", res.Data(), var_.Data()),
3281                                      1200,
3282                                      1200);
3283     fitOutput->Divide(5, 5);
3284 
3285     TCanvas *fitPulls = new TCanvas(Form("fitPulls_%s_%s_%s", res.Data(), var_.Data(), trendPlot->GetName()),
3286                                     Form("fitPulls_%s_%s", res.Data(), var_.Data()),
3287                                     1200,
3288                                     1200);
3289     fitPulls->Divide(5, 5);
3290 
3291     TH1F *residualsPull[myBins_];
3292 
3293     for (Int_t i = 0; i < myBins_; i++) {
3294       TF1 *tmp1 = (TF1 *)residualsPlot[i]->GetListOfFunctions()->FindObject("tmp");
3295       if (tmp1 && residualsPlot[i]->GetEntries() > 0. && residualsPlot[i]->GetMinimum() > 0.) {
3296         fitOutput->cd(i + 1)->SetLogy();
3297       }
3298       fitOutput->cd(i + 1)->SetBottomMargin(0.16);
3299       //fitOutput->cd(i+1)->SetTopMargin(0.05);
3300       //residualsPlot[i]->Sumw2();
3301       MakeNicePlotStyle(residualsPlot[i]);
3302       residualsPlot[i]->SetMarkerStyle(20);
3303       residualsPlot[i]->SetMarkerSize(1.);
3304       residualsPlot[i]->SetStats(0);
3305       //residualsPlot[i]->GetXaxis()->SetRangeUser(-3*(tmp1->GetParameter(1)),3*(tmp1->GetParameter(1)));
3306       residualsPlot[i]->Draw("e1");
3307       residualsPlot[i]->GetYaxis()->UnZoom();
3308 
3309       //std::cout<<"*********************"<<std::endl;
3310       //std::cout<<"fitOutput->cd("<<i+1<<")"<<std::endl;
3311       //std::cout<<"residualsPlot["<<i<<"]->GetTitle() = "<<residualsPlot[i]->GetTitle()<<std::endl;
3312 
3313       // -- for chi2 ----
3314       TPaveText *pt = new TPaveText(0.13, 0.78, 0.33, 0.88, "NDC");
3315       pt->SetFillColor(10);
3316       pt->SetTextColor(1);
3317       pt->SetTextSize(0.07);
3318       pt->SetTextFont(42);
3319       pt->SetTextAlign(22);
3320 
3321       //TF1 *tmp1 = (TF1*)residualsPlot[i]->GetListOfFunctions()->FindObject("tmp");
3322       TString MYOUT;
3323       if (tmp1) {
3324         MYOUT = Form("#chi^{2}/ndf=%.1f", tmp1->GetChisquare() / tmp1->GetNDF());
3325       } else {
3326         MYOUT = "!! no plot !!";
3327       }
3328 
3329       TText *text1 = pt->AddText(MYOUT);
3330       text1->SetTextFont(72);
3331       text1->SetTextColor(kBlue);
3332       pt->Draw("same");
3333 
3334       // -- for bins --
3335 
3336       TPaveText *title = new TPaveText(0.1, 0.93, 0.8, 0.95, "NDC");
3337       title->SetFillColor(10);
3338       title->SetTextColor(1);
3339       title->SetTextSize(0.07);
3340       title->SetTextFont(42);
3341       title->SetTextAlign(22);
3342 
3343       //TText *text2 = title->AddText(residualsPlot[i]->GetTitle());
3344       //text2->SetTextFont(72);
3345       //text2->SetTextColor(kBlue);
3346 
3347       title->Draw("same");
3348 
3349       fitPulls->cd(i + 1);
3350       fitPulls->cd(i + 1)->SetBottomMargin(0.15);
3351       fitPulls->cd(i + 1)->SetLeftMargin(0.15);
3352       fitPulls->cd(i + 1)->SetRightMargin(0.05);
3353 
3354       residualsPull[i] = (TH1F *)residualsPlot[i]->Clone(Form("pull_%s", residualsPlot[i]->GetName()));
3355       for (Int_t nbin = 1; nbin <= residualsPull[i]->GetNbinsX(); nbin++) {
3356         if (residualsPlot[i]->GetBinContent(nbin) != 0 && tmp1) {
3357           residualsPull[i]->SetBinContent(
3358               nbin,
3359               (residualsPlot[i]->GetBinContent(nbin) - tmp1->Eval(residualsPlot[i]->GetBinCenter(nbin))) /
3360                   residualsPlot[i]->GetBinContent(nbin));
3361           residualsPull[i]->SetBinError(nbin, 0.1);
3362         }
3363       }
3364 
3365       TF1 *toDel = (TF1 *)residualsPull[i]->FindObject("tmp");
3366       if (toDel)
3367         residualsPull[i]->GetListOfFunctions()->Remove(toDel);
3368       residualsPull[i]->SetMarkerStyle(20);
3369       residualsPull[i]->SetMarkerSize(1.);
3370       residualsPull[i]->SetStats(0);
3371 
3372       residualsPull[i]->GetYaxis()->SetTitle("(res-fit)/res");
3373       // residualsPull[i]->SetOptTitle(1);
3374       residualsPull[i]->GetXaxis()->SetLabelFont(42);
3375       residualsPull[i]->GetYaxis()->SetLabelFont(42);
3376       residualsPull[i]->GetYaxis()->SetLabelSize(.07);
3377       residualsPull[i]->GetXaxis()->SetLabelSize(.07);
3378       residualsPull[i]->GetYaxis()->SetTitleSize(.07);
3379       residualsPull[i]->GetXaxis()->SetTitleSize(.07);
3380       residualsPull[i]->GetXaxis()->SetTitleOffset(0.9);
3381       residualsPull[i]->GetYaxis()->SetTitleOffset(1.2);
3382       residualsPull[i]->GetXaxis()->SetTitleFont(42);
3383       residualsPull[i]->GetYaxis()->SetTitleFont(42);
3384 
3385       residualsPull[i]->Draw("e1");
3386       residualsPull[i]->GetYaxis()->UnZoom();
3387     }
3388 
3389     TString tpName = trendPlot->GetName();
3390 
3391     TString FitNameToSame = Form("fitOutput_%s", (tpName.ReplaceAll("means_", "").Data()));
3392     //fitOutput->SaveAs(FitNameToSame+".pdf");
3393     //fitOutput->SaveAs(FitNameToSame+".png");
3394     TString PullNameToSave = Form("fitPulls_%s", (tpName.ReplaceAll("means_", "").Data()));
3395     //fitPulls->SaveAs(PullNameToSave+".pdf");
3396     //fitPulls->SaveAs(PullNameToSave+".png");
3397 
3398     if (isDebugMode) {
3399       fitOutput->SaveAs(Form("fitOutput_%s_%s_%s.pdf", res.Data(), var_.Data(), trendPlot->GetName()));
3400       fitOutput->SaveAs(Form("fitOutput_%s.pdf", (((TString)trendPlot->GetName()).ReplaceAll("means_", "")).Data()));
3401       fitPulls->SaveAs(Form("fitPulls_%s.pdf", (((TString)trendPlot->GetName()).ReplaceAll("means_", "")).Data()));
3402       fitOutput->SaveAs(Form("fitOutput_%s.png", (((TString)trendPlot->GetName()).ReplaceAll("means_", "")).Data()));
3403     }
3404 
3405     delete fitOutput;
3406     delete fitPulls;
3407   }
3408 }
3409 
3410 //*************************************************************
3411 void FillMap_old(TH2F *trendMap, TH1F *residualsMapPlot[48][48], params::estimator fitPar_)
3412 //*************************************************************
3413 {
3414   float phiInterval = (360.) / nBins_;
3415   float etaInterval = (etaRange * 2.0) / nBins_;
3416 
3417   switch (fitPar_) {
3418     case params::MEAN: {
3419       for (int i = 0; i < nBins_; ++i) {
3420         char phipositionString[129];
3421         float phiposition = (-180 + i * phiInterval) + (phiInterval / 2);
3422         sprintf(phipositionString, "%.f", phiposition);
3423         trendMap->GetYaxis()->SetBinLabel(i + 1, phipositionString);
3424 
3425         for (int j = 0; j < nBins_; ++j) {
3426           char etapositionString[129];
3427           float etaposition = (-etaRange + j * etaInterval) + (etaInterval / 2);
3428           sprintf(etapositionString, "%.1f", etaposition);
3429 
3430           if (i == 0) {
3431             trendMap->GetXaxis()->SetBinLabel(j + 1, etapositionString);
3432           }
3433 
3434           std::pair<params::measurement, params::measurement> myFit =
3435               std::make_pair(std::make_pair(0., 0.), std::make_pair(0., 0.));
3436 
3437           myFit = fitResiduals(residualsMapPlot[i][j], true);
3438 
3439           float mean_ = myFit.first.first;
3440           float meanErr_ = myFit.first.second;
3441 
3442           trendMap->SetBinContent(j + 1, i + 1, mean_);
3443           trendMap->SetBinError(j + 1, i + 1, meanErr_);
3444         }
3445       }
3446 
3447       break;
3448     }
3449 
3450     case params::WIDTH: {
3451       for (int i = 0; i < nBins_; ++i) {
3452         char phipositionString[129];
3453         float phiposition = (-180 + i * phiInterval) + (phiInterval / 2);
3454         sprintf(phipositionString, "%.f", phiposition);
3455         trendMap->GetYaxis()->SetBinLabel(i + 1, phipositionString);
3456 
3457         for (int j = 0; j < nBins_; ++j) {
3458           char etapositionString[129];
3459           float etaposition = (-etaRange + j * etaInterval) + (etaInterval / 2);
3460           sprintf(etapositionString, "%.1f", etaposition);
3461 
3462           if (i == 0) {
3463             trendMap->GetXaxis()->SetBinLabel(j + 1, etapositionString);
3464           }
3465 
3466           std::pair<params::measurement, params::measurement> myFit =
3467               std::make_pair(std::make_pair(0., 0.), std::make_pair(0., 0.));
3468           myFit = fitResiduals(residualsMapPlot[i][j], true);
3469 
3470           float width_ = myFit.second.first;
3471           float widthErr_ = myFit.second.second;
3472           trendMap->SetBinContent(j + 1, i + 1, width_);
3473           trendMap->SetBinError(j + 1, i + 1, widthErr_);
3474         }
3475       }
3476       break;
3477     }
3478     case params::MEDIAN: {
3479       for (int i = 0; i < nBins_; ++i) {
3480         char phipositionString[129];
3481         float phiposition = (-180 + i * phiInterval) + (phiInterval / 2);
3482         sprintf(phipositionString, "%.f", phiposition);
3483         trendMap->GetYaxis()->SetBinLabel(i + 1, phipositionString);
3484 
3485         for (int j = 0; j < nBins_; ++j) {
3486           char etapositionString[129];
3487           float etaposition = (-etaRange + j * etaInterval) + (etaInterval / 2);
3488           sprintf(etapositionString, "%.1f", etaposition);
3489 
3490           if (i == 0) {
3491             trendMap->GetXaxis()->SetBinLabel(j + 1, etapositionString);
3492           }
3493 
3494           float median_ = getMedian(residualsMapPlot[i][j]).first;
3495           float medianErr_ = getMedian(residualsMapPlot[i][j]).second;
3496           trendMap->SetBinContent(j + 1, i + 1, median_);
3497           trendMap->SetBinError(j + 1, i + 1, medianErr_);
3498         }
3499       }
3500       break;
3501     }
3502     case params::MAD: {
3503       for (int i = 0; i < nBins_; ++i) {
3504         char phipositionString[129];
3505         float phiposition = (-180 + i * phiInterval) + (phiInterval / 2);
3506         sprintf(phipositionString, "%.f", phiposition);
3507         trendMap->GetYaxis()->SetBinLabel(i + 1, phipositionString);
3508 
3509         for (int j = 0; j < nBins_; ++j) {
3510           char etapositionString[129];
3511           float etaposition = (-etaRange + j * etaInterval) + (etaInterval / 2);
3512           sprintf(etapositionString, "%.1f", etaposition);
3513 
3514           if (i == 0) {
3515             trendMap->GetXaxis()->SetBinLabel(j + 1, etapositionString);
3516           }
3517 
3518           float mad_ = getMAD(residualsMapPlot[i][j]).first;
3519           float madErr_ = getMAD(residualsMapPlot[i][j]).second;
3520           trendMap->SetBinContent(j + 1, i + 1, mad_);
3521           trendMap->SetBinError(j + 1, i + 1, madErr_);
3522         }
3523       }
3524       break;
3525     }
3526     default:
3527       std::cout << "FitPVResiduals::FillMap() " << fitPar_ << " unknown estimator!" << std::endl;
3528       break;
3529   }
3530 }
3531 
3532 //*************************************************************
3533 void FillMap(TH2F *trendMap,
3534              std::vector<std::vector<TH1F *> > residualsMapPlot,
3535              params::estimator fitPar_,
3536              const int nBinsX,
3537              const int nBinsY)
3538 //*************************************************************
3539 {
3540   float phiInterval = (360.) / nBinsY;
3541   float etaInterval = 5. / nBinsX;
3542 
3543   for (int i = 0; i < nBinsY; ++i) {
3544     char phipositionString[129];
3545     float phiposition = (-180 + i * phiInterval) + (phiInterval / 2);
3546     sprintf(phipositionString, "%.f", phiposition);
3547 
3548     trendMap->GetYaxis()->SetBinLabel(i + 1, phipositionString);
3549 
3550     for (int j = 0; j < nBinsX; ++j) {
3551       //std::cout<<"(i,j)="<<i<<","<<j<<std::endl;
3552 
3553       char etapositionString[129];
3554       float etaposition = (-etaRange + j * etaInterval) + (etaInterval / 2);
3555       sprintf(etapositionString, "%.1f", etaposition);
3556 
3557       if (i == 0) {
3558         trendMap->GetXaxis()->SetBinLabel(j + 1, etapositionString);
3559       }
3560 
3561       std::pair<params::measurement, params::measurement> myFit =
3562           std::make_pair(std::make_pair(0., 0.), std::make_pair(0., 0.));
3563 
3564       myFit = fitResiduals(residualsMapPlot[i][j], true);
3565 
3566       // check if plot is normalized
3567       bool isNormalized = false;
3568       if (((TString)trendMap->GetName()).Contains("Norm"))
3569         isNormalized = true;
3570 
3571       switch (fitPar_) {
3572         case params::MEAN: {
3573           Double_t mean_ = myFit.first.first;
3574 
3575           // do not allow crazy values
3576           if (!isNormalized)
3577             mean_ = (mean_ > 0.) ? std::min(mean_, 100.) : std::max(mean_, -100.);
3578           else
3579             mean_ = (mean_ > 0.) ? std::min(mean_, 2.) : std::max(mean_, -2.);
3580 
3581           float meanErr_ = myFit.first.second;
3582           //std::cout<<"bin i: "<<i<<" bin j: "<<j<<" mean: "<<mean_<<"+/-"<<meanErr_<<endl;
3583           trendMap->SetBinContent(j + 1, i + 1, mean_);
3584           trendMap->SetBinError(j + 1, i + 1, meanErr_);
3585           break;
3586         }
3587         case params::WIDTH: {
3588           Double_t width_ = myFit.second.first;
3589 
3590           // do not allow crazy values
3591           if (!isNormalized)
3592             width_ = std::min(width_, 1500.);
3593           else
3594             width_ = std::min(width_, 3.);
3595 
3596           float widthErr_ = myFit.second.second;
3597           trendMap->SetBinContent(j + 1, i + 1, width_);
3598           trendMap->SetBinError(j + 1, i + 1, widthErr_);
3599           break;
3600           //std::cout<<"bin i: "<<i<<" bin j: "<<j<<" width: "<<width_<<"+/-"<<widthErr_<<endl;
3601         }
3602         case params::MEDIAN: {
3603           float median_ = getMedian(residualsMapPlot[i][j]).first;
3604           float medianErr_ = getMedian(residualsMapPlot[i][j]).second;
3605           trendMap->SetBinContent(j + 1, i + 1, median_);
3606           trendMap->SetBinError(j + 1, i + 1, medianErr_);
3607           break;
3608         }
3609         case params::MAD: {
3610           float mad_ = getMAD(residualsMapPlot[i][j]).first;
3611           float madErr_ = getMAD(residualsMapPlot[i][j]).second;
3612           trendMap->SetBinContent(j + 1, i + 1, mad_);
3613           trendMap->SetBinError(j + 1, i + 1, madErr_);
3614           break;
3615         }
3616         default:
3617           std::cout << "FitPVResiduals::FillMap() " << fitPar_ << " unknown estimator!" << std::endl;
3618           break;
3619       }  // closes the switch statement
3620     }    // closes loop on eta bins
3621   }      // cloeses loop on phi bins
3622 }
3623 
3624 /*--------------------------------------------------------------------*/
3625 void MakeNiceTrendPlotStyle(TH1 *hist, Int_t color, Int_t style)
3626 /*--------------------------------------------------------------------*/
3627 {
3628   hist->SetStats(kFALSE);
3629   hist->SetLineWidth(2);
3630   hist->GetXaxis()->CenterTitle(true);
3631   hist->GetYaxis()->CenterTitle(true);
3632   hist->GetXaxis()->SetTitleFont(42);
3633   hist->GetYaxis()->SetTitleFont(42);
3634   hist->GetXaxis()->SetTitleSize(0.065);
3635   hist->GetYaxis()->SetTitleSize(0.065);
3636   hist->GetXaxis()->SetTitleOffset(1.0);
3637   hist->GetYaxis()->SetTitleOffset(1.2);
3638   hist->GetXaxis()->SetLabelFont(42);
3639   hist->GetYaxis()->SetLabelFont(42);
3640   hist->GetYaxis()->SetLabelSize(.05);
3641   hist->GetXaxis()->SetLabelSize(.07);
3642   //hist->GetXaxis()->SetNdivisions(505);
3643   if (color != 8) {
3644     hist->SetMarkerSize(1.0);
3645   } else {
3646     hist->SetLineWidth(3);
3647     hist->SetMarkerSize(0.0);
3648   }
3649   hist->SetMarkerStyle(style);
3650   hist->SetLineColor(color);
3651   hist->SetMarkerColor(color);
3652 }
3653 
3654 /*--------------------------------------------------------------------*/
3655 void MakeNicePlotStyle(TH1 *hist)
3656 /*--------------------------------------------------------------------*/
3657 {
3658   hist->SetStats(kFALSE);
3659   hist->SetLineWidth(2);
3660   hist->GetXaxis()->SetNdivisions(505);
3661   hist->GetXaxis()->CenterTitle(true);
3662   hist->GetYaxis()->CenterTitle(true);
3663   hist->GetXaxis()->SetTitleFont(42);
3664   hist->GetYaxis()->SetTitleFont(42);
3665   hist->GetXaxis()->SetTitleSize(0.07);
3666   hist->GetYaxis()->SetTitleSize(0.07);
3667   hist->GetXaxis()->SetTitleOffset(0.9);
3668   hist->GetYaxis()->SetTitleOffset(1.3);
3669   hist->GetXaxis()->SetLabelFont(42);
3670   hist->GetYaxis()->SetLabelFont(42);
3671   hist->GetYaxis()->SetLabelSize(.07);
3672   hist->GetXaxis()->SetLabelSize(.07);
3673 }
3674 
3675 /*--------------------------------------------------------------------*/
3676 void MakeNiceMapStyle(TH2 *hist)
3677 /*--------------------------------------------------------------------*/
3678 {
3679   hist->SetStats(kFALSE);
3680   hist->GetXaxis()->CenterTitle(true);
3681   hist->GetYaxis()->CenterTitle(true);
3682   hist->GetZaxis()->CenterTitle(true);
3683   hist->GetXaxis()->SetTitleFont(42);
3684   hist->GetYaxis()->SetTitleFont(42);
3685   hist->GetXaxis()->LabelsOption("v");
3686   hist->GetZaxis()->SetTitleFont(42);
3687   hist->GetXaxis()->SetTitleSize(0.06);
3688   hist->GetYaxis()->SetTitleSize(0.06);
3689   hist->GetZaxis()->SetTitleSize(0.06);
3690   hist->GetXaxis()->SetTitleOffset(1.1);
3691   hist->GetZaxis()->SetTitleOffset(1.1);
3692   hist->GetYaxis()->SetTitleOffset(1.0);
3693   hist->GetXaxis()->SetLabelFont(42);
3694   hist->GetYaxis()->SetLabelFont(42);
3695   hist->GetZaxis()->SetLabelFont(42);
3696   hist->GetYaxis()->SetLabelSize(.05);
3697   hist->GetXaxis()->SetLabelSize(.05);
3698   hist->GetZaxis()->SetLabelSize(.05);
3699 }
3700 
3701 /*--------------------------------------------------------------------*/
3702 std::pair<TH2F *, TH2F *> trimTheMap(TH2 *hist) {
3703   /*--------------------------------------------------------------------*/
3704 
3705   Int_t nXCells = hist->GetNbinsX();
3706   Int_t nYCells = hist->GetNbinsY();
3707   Int_t nCells = nXCells * nYCells;
3708 
3709   Double_t min = 9999.;
3710   Double_t max = -9999.;
3711 
3712   for (Int_t nX = 1; nX <= nXCells; nX++) {
3713     for (Int_t nY = 1; nY <= nYCells; nY++) {
3714       Double_t binContent = hist->GetBinContent(nX, nY);
3715       if (binContent > max)
3716         max = binContent;
3717       if (binContent < min)
3718         min = binContent;
3719     }
3720   }
3721 
3722   TH1F *histContentByCell =
3723       new TH1F(Form("histContentByCell_%s", hist->GetName()), "histContentByCell", nCells, min, max);
3724 
3725   for (Int_t nX = 1; nX <= nXCells; nX++) {
3726     for (Int_t nY = 1; nY <= nYCells; nY++) {
3727       histContentByCell->Fill(hist->GetBinContent(nX, nY));
3728     }
3729   }
3730 
3731   Double_t theMeanOfCells = histContentByCell->GetMean();
3732   Double_t theRMSOfCells = histContentByCell->GetRMS();
3733   params::measurement theMAD = getMAD(histContentByCell);
3734 
3735   if (isDebugMode) {
3736     std::cout << std::setw(24) << std::left << hist->GetName() << "| mean: " << std::setw(10) << std::setprecision(4)
3737               << theMeanOfCells << "| min: " << std::setw(10) << std::setprecision(4) << min
3738               << "| max: " << std::setw(10) << std::setprecision(4) << max << "| rms: " << std::setw(10)
3739               << std::setprecision(4) << theRMSOfCells << "| mad: " << std::setw(10) << std::setprecision(4)
3740               << theMAD.first << std::endl;
3741   }
3742 
3743   TCanvas *cCheck = new TCanvas(Form("cCheck_%s", hist->GetName()), Form("cCheck_%s", hist->GetName()), 1200, 1000);
3744 
3745   cCheck->Divide(2, 2);
3746   for (Int_t i = 1; i <= 4; i++) {
3747     cCheck->cd(i)->SetBottomMargin(0.13);
3748     cCheck->cd(i)->SetLeftMargin(0.12);
3749     if (i % 2 == 1)
3750       cCheck->cd(i)->SetRightMargin(0.19);
3751     else
3752       cCheck->cd(i)->SetRightMargin(0.07);
3753     cCheck->cd(i)->SetTopMargin(0.08);
3754   }
3755 
3756   cCheck->cd(1);
3757   hist->SetStats(kFALSE);
3758   hist->Draw("colz");
3759   //makeNewPairOfAxes(hist);
3760 
3761   cCheck->cd(2)->SetLogy();
3762   MakeNicePlotStyle(histContentByCell);
3763   histContentByCell->SetStats(kTRUE);
3764   histContentByCell->GetYaxis()->SetTitleOffset(0.9);
3765   histContentByCell->Draw();
3766 
3767   //Double_t theNewMin = theMeanOfCells-theRMSOfCells;
3768   //Double_t theNewMax = theMeanOfCells+theRMSOfCells;
3769 
3770   Double_t theNewMin = theMeanOfCells - theMAD.first * 3;
3771   Double_t theNewMax = theMeanOfCells + theMAD.first * 3;
3772 
3773   TArrow *l0 =
3774       new TArrow(theMeanOfCells, cCheck->GetUymin(), theMeanOfCells, histContentByCell->GetMaximum(), 0.3, "|>");
3775   l0->SetAngle(60);
3776   l0->SetLineColor(kRed);
3777   l0->SetLineWidth(4);
3778   l0->Draw("same");
3779 
3780   TArrow *l1 = new TArrow(theNewMin, cCheck->GetUymin(), theNewMin, histContentByCell->GetMaximum(), 0.3, "|>");
3781   l1->SetAngle(60);
3782   l1->SetLineColor(kBlue);
3783   l1->SetLineWidth(4);
3784   l1->Draw("same");
3785 
3786   TArrow *l2 = new TArrow(theNewMax, cCheck->GetUymin(), theNewMax, histContentByCell->GetMaximum(), 0.3, "|>");
3787   l2->SetAngle(60);
3788   l2->SetLineColor(kBlue);
3789   l2->SetLineWidth(4);
3790   l2->Draw("same");
3791 
3792   TH2F *histoTrimmed = new TH2F(Form("%s_trimmed", hist->GetName()),
3793                                 Form("Trimmed %s;%s;%s;%s",
3794                                      hist->GetTitle(),
3795                                      hist->GetXaxis()->GetTitle(),
3796                                      hist->GetYaxis()->GetTitle(),
3797                                      hist->GetZaxis()->GetTitle()),
3798                                 hist->GetNbinsX(),
3799                                 hist->GetXaxis()->GetXmin(),
3800                                 hist->GetXaxis()->GetXmax(),
3801                                 hist->GetNbinsY(),
3802                                 hist->GetYaxis()->GetXmin(),
3803                                 hist->GetYaxis()->GetXmax());
3804 
3805   TH2F *histoMissed = new TH2F(Form("%s_Missed", hist->GetName()),
3806                                Form("Missed %s", hist->GetTitle()),
3807                                hist->GetNbinsX(),
3808                                hist->GetXaxis()->GetXmin(),
3809                                hist->GetXaxis()->GetXmax(),
3810                                hist->GetNbinsY(),
3811                                hist->GetYaxis()->GetXmin(),
3812                                hist->GetYaxis()->GetXmax());
3813 
3814   for (Int_t nX = 1; nX <= nXCells; nX++) {
3815     for (Int_t nY = 1; nY <= nYCells; nY++) {
3816       Double_t binContent = hist->GetBinContent(nX, nY);
3817       Double_t binError = hist->GetBinError(nX, nY);
3818 
3819       if (binContent == 0. && binError == 0.) {
3820         histoMissed->SetBinContent(nX, nY, 1);
3821       } else if (binContent <= theNewMin) {
3822         histoTrimmed->SetBinContent(nX, nY, theNewMin);
3823       } else if (binContent >= theNewMax) {
3824         histoTrimmed->SetBinContent(nX, nY, theNewMax);
3825       } else {
3826         histoTrimmed->SetBinContent(nX, nY, binContent);
3827       }
3828     }
3829   }
3830 
3831   cCheck->cd(3);
3832   histoTrimmed->SetStats(kFALSE);
3833   histoTrimmed->Draw("COLZ1");
3834   histoMissed->SetFillColor(kRed);
3835   gStyle->SetPaintTextFormat("0.1f");
3836   //histoMissed->SetMarkerSize(1.8);
3837   histoMissed->SetFillColor(kMagenta);
3838   histoMissed->SetMarkerColor(kMagenta);
3839   histoMissed->Draw("boxsame");
3840   //makeNewPairOfAxes(histoTrimmed);
3841 
3842   cCheck->cd(4);
3843   histoMissed->SetStats(kFALSE);
3844   histoMissed->Draw("box");
3845   makeNewPairOfAxes(histoMissed);
3846 
3847   if (isDebugMode) {
3848     cCheck->SaveAs(Form("cCheck_%s.png", hist->GetName()));
3849     cCheck->SaveAs(Form("cCheck_%s.pdf", hist->GetName()));
3850 
3851     std::cout << "histo:" << std::setw(25) << hist->GetName() << " old min: " << std::setw(10) << hist->GetMinimum()
3852               << " old max: " << std::setw(10) << hist->GetMaximum();
3853     std::cout << " | new min: " << std::setw(15) << hist->GetMinimum() << " new max: " << std::setw(10)
3854               << hist->GetMaximum() << std::endl;
3855   }
3856 
3857   delete histContentByCell;
3858   // hist->GetZaxis()->SetRangeUser(1.0001*theNewMin,0.999*theNewMax);
3859   // hist->SetMinimum(1.001*theNewMin);
3860   // hist->SetMaximum(0.999*theNewMax);
3861   delete cCheck;
3862 
3863   return std::make_pair(histoTrimmed, histoMissed);
3864 }
3865 
3866 /*--------------------------------------------------------------------*/
3867 void setStyle() {
3868   /*--------------------------------------------------------------------*/
3869 
3870   writeExtraText = true;  // if extra text
3871   lumi_13TeV = "p-p collisions";
3872   extraText = "Internal";
3873 
3874   TH1::StatOverflows(kTRUE);
3875   gStyle->SetOptTitle(0);
3876   gStyle->SetOptStat("e");
3877   //gStyle->SetPadTopMargin(0.05);
3878   //gStyle->SetPadBottomMargin(0.15);
3879   //gStyle->SetPadLeftMargin(0.17);
3880   //gStyle->SetPadRightMargin(0.02);
3881   gStyle->SetPadBorderMode(0);
3882   gStyle->SetTitleFillColor(10);
3883   gStyle->SetTitleFont(42);
3884   gStyle->SetTitleColor(1);
3885   gStyle->SetTitleTextColor(1);
3886   gStyle->SetTitleFontSize(0.06);
3887   gStyle->SetTitleBorderSize(0);
3888   gStyle->SetStatColor(kWhite);
3889   gStyle->SetStatFont(42);
3890   gStyle->SetStatFontSize(0.05);  ///---> gStyle->SetStatFontSize(0.025);
3891   gStyle->SetStatTextColor(1);
3892   gStyle->SetStatFormat("6.4g");
3893   gStyle->SetStatBorderSize(1);
3894   gStyle->SetPadTickX(1);  // To get tick marks on the opposite side of the frame
3895   gStyle->SetPadTickY(1);
3896   gStyle->SetPadBorderMode(0);
3897   gStyle->SetOptFit(1);
3898   gStyle->SetNdivisions(510);
3899 
3900   // this is the standard palette
3901   const Int_t NRGBs = 5;
3902   const Int_t NCont = 255;
3903 
3904   Double_t stops[NRGBs] = {0.00, 0.34, 0.61, 0.84, 1.00};
3905   Double_t red[NRGBs] = {0.00, 0.00, 0.87, 1.00, 0.51};
3906   Double_t green[NRGBs] = {0.00, 0.81, 1.00, 0.20, 0.00};
3907   Double_t blue[NRGBs] = {0.51, 1.00, 0.12, 0.00, 0.00};
3908   TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
3909   gStyle->SetNumberContours(NCont);
3910 
3911   /*
3912   // try an alternative palette
3913   const Int_t NRGBs = 6;
3914   const Int_t NCont = 999;
3915 
3916   Double_t stops[NRGBs] = { 0.00, 0.1, 0.34, 0.61, 0.84, 1.00 };
3917   Double_t red[NRGBs]   = { 0.99, 0.0, 0.00, 0.87, 1.00, 0.51 };
3918   Double_t green[NRGBs] = { 0.00, 0.0, 0.81, 1.00, 0.20, 0.00 };
3919   Double_t blue[NRGBs]  = { 0.99, 0.0, 1.00, 0.12, 0.00, 0.00 };
3920 
3921   TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
3922   gStyle->SetNumberContours(NCont);
3923 
3924   */
3925 
3926   /*
3927   const Int_t NRGBs = 9;
3928   const Int_t NCont = 255;
3929  
3930   Double_t stops[NRGBs] = { 0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000};
3931 
3932   // dark body radiator
3933   // Double_t red[NRGBs]   = { 0./255., 45./255., 99./255., 156./255., 212./255., 230./255., 237./255., 234./255., 242./255.};
3934   // Double_t green[NRGBs] = { 0./255.,  0./255.,  0./255.,  45./255., 101./255., 168./255., 238./255., 238./255., 243./255.};
3935   // Double_t blue[NRGBs]  = { 0./255.,  1./255.,  1./255.,   3./255.,   9./255.,   8./255.,  11./255.,  95./255., 230./255.};
3936   
3937   // printable on grey
3938   //Double_t red[9]   = { 0./255.,   0./255.,   0./255.,  70./255., 148./255., 231./255., 235./255., 237./255., 244./255.};
3939   //Double_t green[9] = { 0./255.,   0./255.,   0./255.,   0./255.,   0./255.,  69./255.,  67./255., 216./255., 244./255.};
3940   //Double_t blue[9]  = { 0./255., 102./255., 228./255., 231./255., 177./255., 124./255., 137./255.,  20./255., 244./255.};
3941 
3942   // thermometer
3943   //Double_t red[9]   = {  34./255.,  70./255., 129./255., 187./255., 225./255., 226./255., 216./255., 193./255., 179./255.};
3944   //Double_t green[9] = {  48./255.,  91./255., 147./255., 194./255., 226./255., 229./255., 196./255., 110./255.,  12./255.};
3945   //Double_t blue[9]  = { 234./255., 212./255., 216./255., 224./255., 206./255., 110./255.,  53./255.,  40./255.,  29./255.};
3946 
3947   // visible spectrum
3948   Double_t red[9]   = { 18./255.,  72./255.,   5./255.,  23./255.,  29./255., 201./255., 200./255., 98./255., 29./255.};
3949   Double_t green[9] = {  0./255.,   0./255.,  43./255., 167./255., 211./255., 117./255.,   0./255.,  0./255.,  0./255.};
3950   Double_t blue[9]  = { 51./255., 203./255., 177./255.,  26./255.,  10./255.,   9./255.,   8./255.,  3./255.,  0./255.};
3951 
3952   TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
3953   gStyle->SetNumberContours(NCont);
3954   */
3955 }
3956 
3957 /*--------------------------------------------------------------------*/
3958 TH1F *DrawZero(TH1F *hist, Int_t nbins, Double_t lowedge, Double_t highedge, Int_t iter)
3959 /*--------------------------------------------------------------------*/
3960 {
3961   TH1F *hzero = new TH1F(
3962       Form("hzero_%s_%i", hist->GetName(), iter), Form("hzero_%s_%i", hist->GetName(), iter), nbins, lowedge, highedge);
3963   for (Int_t i = 0; i < hzero->GetNbinsX(); i++) {
3964     hzero->SetBinContent(i, 0.);
3965     hzero->SetBinError(i, 0.);
3966   }
3967   hzero->SetLineWidth(2);
3968   hzero->SetLineStyle(9);
3969   hzero->SetLineColor(kMagenta);
3970 
3971   return hzero;
3972 }
3973 
3974 /*--------------------------------------------------------------------*/
3975 TH1F *DrawConstant(TH1F *hist, Int_t nbins, Double_t lowedge, Double_t highedge, Int_t iter, Double_t theConst)
3976 /*--------------------------------------------------------------------*/
3977 {
3978   TH1F *hzero = new TH1F(Form("hconst_%s_%i", hist->GetName(), iter),
3979                          Form("hconst_%s_%i", hist->GetName(), iter),
3980                          nbins,
3981                          lowedge,
3982                          highedge);
3983   for (Int_t i = 0; i <= hzero->GetNbinsX(); i++) {
3984     hzero->SetBinContent(i, theConst);
3985     hzero->SetBinError(i, 0.);
3986   }
3987   hzero->SetLineWidth(2);
3988   hzero->SetLineStyle(9);
3989   hzero->SetLineColor(kMagenta);
3990 
3991   return hzero;
3992 }
3993 
3994 /*--------------------------------------------------------------------*/
3995 void makeNewXAxis(TH1F *h)
3996 /*--------------------------------------------------------------------*/
3997 {
3998   TString myTitle = h->GetName();
3999   float axmin = -999;
4000   float axmax = 999.;
4001   int ndiv = 510;
4002   if (myTitle.Contains("eta")) {
4003     axmin = -etaRange;
4004     axmax = etaRange;
4005     ndiv = 505;
4006   } else if (myTitle.Contains("phi")) {
4007     axmin = -TMath::Pi();
4008     axmax = TMath::Pi();
4009     ndiv = 510;
4010   } else if (myTitle.Contains("pT")) {
4011     axmin = minPt_;
4012     axmax = maxPt_;
4013     ndiv = 510;
4014   } else if (myTitle.Contains("ladder")) {
4015     axmin = 0.5;
4016     axmax = nLadders_ + 0.5;
4017   } else if (myTitle.Contains("modZ")) {
4018     axmin = 0.5;
4019     axmax = nModZ_ + 0.5;
4020   } else if (myTitle.Contains("h_probe")) {
4021     ndiv = 505;
4022     axmin = h->GetXaxis()->GetBinCenter(h->GetXaxis()->GetFirst());
4023     axmax = h->GetXaxis()->GetBinCenter(h->GetXaxis()->GetLast());
4024   } else {
4025     std::cout << "unrecognized variable for histogram title: " << myTitle << std::endl;
4026   }
4027 
4028   // Remove the current axis
4029   h->GetXaxis()->SetLabelOffset(999);
4030   h->GetXaxis()->SetTickLength(0);
4031 
4032   // Redraw the new axis
4033   gPad->Update();
4034 
4035   TGaxis *newaxis =
4036       new TGaxis(gPad->GetUxmin(), gPad->GetUymin(), gPad->GetUxmax(), gPad->GetUymin(), axmin, axmax, ndiv, "SDH");
4037 
4038   TGaxis *newaxisup =
4039       new TGaxis(gPad->GetUxmin(), gPad->GetUymax(), gPad->GetUxmax(), gPad->GetUymax(), axmin, axmax, ndiv, "-SDH");
4040 
4041   newaxis->SetLabelOffset(0.02);
4042   newaxis->SetLabelFont(42);
4043   newaxis->SetLabelSize(0.05);
4044 
4045   newaxisup->SetLabelOffset(-0.02);
4046   newaxisup->SetLabelFont(42);
4047   newaxisup->SetLabelSize(0);
4048 
4049   newaxis->Draw();
4050   newaxisup->Draw();
4051 }
4052 
4053 /*--------------------------------------------------------------------*/
4054 void makeNewPairOfAxes(TH2F *h)
4055 /*--------------------------------------------------------------------*/
4056 {
4057   TString myTitle = h->GetName();
4058   // fake defaults
4059   float axmin = -999;
4060   float axmax = 999.;
4061   float aymin = -999;
4062   float aymax = 999.;
4063   int ndivx = h->GetXaxis()->GetNdivisions();
4064   int ndivy = h->GetYaxis()->GetNdivisions();
4065 
4066   if (!myTitle.Contains("L1Map")) {
4067     ndivx = 505;
4068     ndivy = 510;
4069     axmin = -etaRange;
4070     axmax = etaRange;
4071     aymin = -TMath::Pi();
4072     aymax = TMath::Pi();
4073   } else {
4074     // this is a L1 map
4075     axmin = 0.5;
4076     axmax = nModZ_ + 0.5;
4077     aymin = 0.5;
4078     aymax = nLadders_ + 0.5;
4079   }
4080 
4081   // Remove the current axis
4082   h->GetXaxis()->SetLabelOffset(999);
4083   h->GetXaxis()->SetTickLength(0);
4084 
4085   h->GetYaxis()->SetLabelOffset(999);
4086   h->GetYaxis()->SetTickLength(0);
4087 
4088   // Redraw the new axis
4089   gPad->Update();
4090 
4091   TGaxis *newXaxis =
4092       new TGaxis(gPad->GetUxmin(), gPad->GetUymin(), gPad->GetUxmax(), gPad->GetUymin(), axmin, axmax, ndivx, "SDH");
4093 
4094   TGaxis *newXaxisup =
4095       new TGaxis(gPad->GetUxmin(), gPad->GetUymax(), gPad->GetUxmax(), gPad->GetUymax(), axmin, axmax, ndivx, "-SDH");
4096 
4097   TGaxis *newYaxisR =
4098       new TGaxis(gPad->GetUxmin(), gPad->GetUymin(), gPad->GetUxmin(), gPad->GetUymax(), aymin, aymax, ndivy, "SDH");
4099 
4100   TGaxis *newYaxisL =
4101       new TGaxis(gPad->GetUxmax(), gPad->GetUymin(), gPad->GetUxmax(), gPad->GetUymax(), aymin, aymax, ndivy, "-SDH");
4102 
4103   newXaxis->SetLabelOffset(0.02);
4104   newXaxis->SetLabelFont(42);
4105   newXaxis->SetLabelSize(0.055);
4106 
4107   newXaxisup->SetLabelOffset(-0.02);
4108   newXaxisup->SetLabelFont(42);
4109   newXaxisup->SetLabelSize(0);
4110 
4111   newXaxis->Draw();
4112   newXaxisup->Draw();
4113 
4114   newYaxisR->SetLabelOffset(0.02);
4115   newYaxisR->SetLabelFont(42);
4116   newYaxisR->SetLabelSize(0.055);
4117 
4118   newYaxisL->SetLabelOffset(-0.02);
4119   newYaxisL->SetLabelFont(42);
4120   newYaxisL->SetLabelSize(0);
4121 
4122   newYaxisR->Draw();
4123   newYaxisL->Draw();
4124 }
4125 
4126 /*--------------------------------------------------------------------*/
4127 Double_t fDLine(Double_t *x, Double_t *par)
4128 /*--------------------------------------------------------------------*/
4129 {
4130   if (x[0] < _boundSx && x[0] > _boundDx) {
4131     TF1::RejectPoint();
4132     return 0;
4133   }
4134   return par[0];
4135 }
4136 
4137 /*--------------------------------------------------------------------*/
4138 Double_t fULine(Double_t *x, Double_t *par)
4139 /*--------------------------------------------------------------------*/
4140 {
4141   if (x[0] >= _boundSx && x[0] <= _boundDx) {
4142     TF1::RejectPoint();
4143     return 0;
4144   }
4145   return par[0];
4146 }
4147 
4148 /*--------------------------------------------------------------------*/
4149 void FitULine(TH1 *hist)
4150 /*--------------------------------------------------------------------*/
4151 {
4152   // define fitting function
4153   TF1 func1("lineUp", fULine, _boundMin, _boundMax, 1);
4154   //TF1 func1("lineUp","pol0",-0.5,11.5);
4155 
4156   if (0 == hist->Fit(&func1, "QR")) {
4157     if (hist->GetFunction(func1.GetName())) {  // Take care that it is later on drawn:
4158       hist->GetFunction(func1.GetName())->ResetBit(TF1::kNotDraw);
4159     }
4160     //std::cout<<"FitPVResiduals() fit Up done!"<<std::endl;
4161   }
4162 }
4163 
4164 /*--------------------------------------------------------------------*/
4165 void FitDLine(TH1 *hist)
4166 /*--------------------------------------------------------------------*/
4167 {
4168   // define fitting function
4169   // TF1 func1("lineDown",fDLine,-0.5,11.5,1);
4170 
4171   TF1 func2("lineDown", "pol0", _boundSx, _boundDx);
4172   func2.SetRange(_boundSx, _boundDx);
4173 
4174   if (0 == hist->Fit(&func2, "QR")) {
4175     if (hist->GetFunction(func2.GetName())) {  // Take care that it is later on drawn:
4176       hist->GetFunction(func2.GetName())->ResetBit(TF1::kNotDraw);
4177     }
4178     // std::cout<<"FitPVResiduals() fit Down done!"<<std::endl;
4179   }
4180 }
4181 
4182 /*--------------------------------------------------------------------*/
4183 void MakeNiceTF1Style(TF1 *f1, Int_t color)
4184 /*--------------------------------------------------------------------*/
4185 {
4186   f1->SetLineColor(color);
4187   f1->SetLineWidth(3);
4188   f1->SetLineStyle(2);
4189 }
4190 
4191 /*--------------------------------------------------------------------*/
4192 params::measurement getTheRangeUser(TH1F *thePlot, Limits *lims, bool tag)
4193 /*--------------------------------------------------------------------*/
4194 {
4195   TString theTitle = thePlot->GetName();
4196   theTitle.ToLower();
4197 
4198   /*
4199     Double_t m_dxyPhiMax     = 40;
4200     Double_t m_dzPhiMax      = 40;
4201     Double_t m_dxyEtaMax     = 40;
4202     Double_t m_dzEtaMax      = 40;
4203     
4204     Double_t m_dxyPhiNormMax = 0.5;
4205     Double_t m_dzPhiNormMax  = 0.5;
4206     Double_t m_dxyEtaNormMax = 0.5;
4207     Double_t m_dzEtaNormMax  = 0.5;
4208     
4209     Double_t w_dxyPhiMax     = 150;
4210     Double_t w_dzPhiMax      = 150;
4211     Double_t w_dxyEtaMax     = 150;
4212     Double_t w_dzEtaMax      = 1000;
4213     
4214     Double_t w_dxyPhiNormMax = 1.8;
4215     Double_t w_dzPhiNormMax  = 1.8;
4216     Double_t w_dxyEtaNormMax = 1.8;
4217     Double_t w_dzEtaNormMax  = 1.8;   
4218   */
4219 
4220   params::measurement result;
4221 
4222   if (theTitle.Contains("norm")) {
4223     if (theTitle.Contains("means")) {
4224       if (theTitle.Contains("dxy") || theTitle.Contains("dx") || theTitle.Contains("dy")) {
4225         if (theTitle.Contains("phi") || theTitle.Contains("pT") || theTitle.Contains("ladder")) {
4226           result = std::make_pair(-lims->get_dxyPhiNormMax().first, lims->get_dxyPhiNormMax().first);
4227         } else if (theTitle.Contains("eta") || theTitle.Contains("mod")) {
4228           result = std::make_pair(-lims->get_dxyEtaNormMax().first, lims->get_dxyEtaNormMax().first);
4229         } else {
4230           result = std::make_pair(-0.8, 0.8);
4231         }
4232       } else if (theTitle.Contains("dz")) {
4233         if (theTitle.Contains("phi") || theTitle.Contains("pT") || theTitle.Contains("ladder")) {
4234           result = std::make_pair(-lims->get_dzPhiNormMax().first, lims->get_dzPhiNormMax().first);
4235         } else if (theTitle.Contains("eta") || theTitle.Contains("mod")) {
4236           result = std::make_pair(-lims->get_dzEtaNormMax().first, lims->get_dzEtaNormMax().first);
4237         } else {
4238           result = std::make_pair(-0.8, 0.8);
4239         }
4240       }
4241     } else if (theTitle.Contains("widths")) {
4242       if (theTitle.Contains("dxy") || theTitle.Contains("dx") || theTitle.Contains("dy")) {
4243         if (theTitle.Contains("phi") || theTitle.Contains("ladder")) {
4244           result = std::make_pair(0., lims->get_dxyPhiNormMax().second);
4245         } else if (theTitle.Contains("eta") || theTitle.Contains("mod")) {
4246           result = std::make_pair(0., lims->get_dxyEtaNormMax().second);
4247         } else {
4248           result = std::make_pair(0., 2.);
4249         }
4250       } else if (theTitle.Contains("dz")) {
4251         if (theTitle.Contains("phi") || theTitle.Contains("ladder")) {
4252           result = std::make_pair(0., lims->get_dzPhiNormMax().second);
4253         } else if (theTitle.Contains("eta") || theTitle.Contains("mod")) {
4254           result = std::make_pair(0., lims->get_dzEtaNormMax().second);
4255         } else {
4256           result = std::make_pair(0., 2.);
4257         }
4258       }
4259     }
4260   } else {
4261     if (theTitle.Contains("means")) {
4262       if (theTitle.Contains("dxy") || theTitle.Contains("dx") || theTitle.Contains("dy")) {
4263         if (theTitle.Contains("phi") || theTitle.Contains("pT") || theTitle.Contains("ladder")) {
4264           result = std::make_pair(-lims->get_dxyPhiMax().first, lims->get_dxyPhiMax().first);
4265         } else if (theTitle.Contains("eta") || theTitle.Contains("mod")) {
4266           result = std::make_pair(-lims->get_dxyEtaMax().first, lims->get_dxyEtaMax().first);
4267         } else {
4268           result = std::make_pair(-40., 40.);
4269         }
4270       } else if (theTitle.Contains("dz")) {
4271         if (theTitle.Contains("phi") || theTitle.Contains("pT") || theTitle.Contains("ladder")) {
4272           result = std::make_pair(-lims->get_dzPhiMax().first, lims->get_dzPhiMax().first);
4273         } else if (theTitle.Contains("eta") || theTitle.Contains("mod")) {
4274           result = std::make_pair(-lims->get_dzEtaMax().first, lims->get_dzEtaMax().first);
4275         } else {
4276           result = std::make_pair(-80., 80.);
4277         }
4278       }
4279     } else if (theTitle.Contains("widths")) {
4280       if (theTitle.Contains("dxy") || theTitle.Contains("dx") || theTitle.Contains("dy")) {
4281         if (theTitle.Contains("phi") || theTitle.Contains("ladder")) {
4282           result = std::make_pair(0., lims->get_dxyPhiMax().second);
4283         } else if (theTitle.Contains("eta") || theTitle.Contains("mod")) {
4284           result = std::make_pair(0., lims->get_dxyEtaMax().second);
4285         } else {
4286           result = std::make_pair(0., 150.);
4287         }
4288       } else if (theTitle.Contains("dz")) {
4289         if (theTitle.Contains("phi") || theTitle.Contains("ladder")) {
4290           result = std::make_pair(0., lims->get_dzPhiMax().second);
4291         } else if (theTitle.Contains("eta") || theTitle.Contains("mod")) {
4292           result = std::make_pair(0., lims->get_dzEtaMax().second);
4293         } else {
4294           result = std::make_pair(0., 300.);
4295         }
4296       }
4297     }
4298   }
4299 
4300   if (tag)
4301     std::cout << theTitle << " " << result.first << " " << result.second << std::endl;
4302   return result;
4303 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!