Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-15 23:06:23 

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