Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!