OfflineValidation/macros/FitPVResiduals.C

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