OfflineValidation/src/GeometryComparisonPlotter.cc

0001 //#include "GeometryComparisonPlotter.h"
0002 #include "Alignment/OfflineValidation/interface/GeometryComparisonPlotter.h"
0003
0004 /***********************************************************************************/
0005 /*                            GEOMETRY COMPARISON PLOTTER                          */
0006 /* See the talk of 15 January 2015 for short documentation and the example script. */
0007 /* This code is highly commented if need be to upgrade it.                         */
0008 /* Any further question is to be asked to Patrick Connor (patrick.connor@desy.de). */
0009 /*                                             Thanks a million <3                 */
0010 /***********************************************************************************/
0011
0012 // NOTE: look for "TO DO" as a keyword to now what should be upgraded in later versions....
0013
0014 // modes
0015 //#define TALKATIVE  // get some comments while processing
0016 //#define DEBUG     // get a lot of comments while processing + canvases -> resource-consuming!
0017
0018 // MACROS
0019 #define INSIDE_VECTOR(vector)                          \
0020   std::cout << #vector << "={";                        \
0021   for (unsigned int i = 0; i < vector.size() - 1; i++) \
0022     std::cout << vector[i] << ",";                     \
0023   std::cout << vector.back() << "}";
0024 #define CHECK_MAP_CONTENT(m, type)                                            \
0025   for (std::map<TString, type>::iterator it = m.begin(); it != m.end(); it++) \
0026     std::cout << __FILE__ << ":" << __LINE__ << ":Info: " << #m << "[" << it->first << "]=" << it->second << std::endl;
0027 /*
0028     TString _sublevel_names[NB_SUBLEVELS],
0029             _output_directory,
0030             _output_filename,
0031             _print_option,
0032             _module_plot_option,
0033             _alignment_name,
0034             _reference_name;
0035     bool _print,
0036          _legend,
0037          _write,
0038          _print_only_global,
0039          _make_profile_plots,
0040          _batchMode,
0041          _1dModule,
0042          _2dModule;
0043     int _levelCut,
0044         _grid_x,
0045         _grid_y,
0046         _window_width,
0047         _window_height;
0048  */
0049
0050 // CONSTRUCTOR AND DESTRUCTOR
0051 GeometryComparisonPlotter::GeometryComparisonPlotter(TString tree_file_name,
0052                                                      TString output_directory,
0053                                                      TString modulesToPlot,
0054                                                      TString alignmentName,
0055                                                      TString referenceName,
0056                                                      bool printOnlyGlobal,
0057                                                      bool makeProfilePlots,
0058                                                      int canvas_idx)
0059     : _output_directory(output_directory + TString(output_directory.EndsWith("/") ? "" : "/")),
0060       _output_filename("comparison.root"),
0061       _print_option("pdf"),
0062       _module_plot_option(modulesToPlot),
0063       _alignment_name(alignmentName),
0064       _reference_name(referenceName),
0065       _print(true),   // print the graphs in a file (e.g. pdf)
0066       _legend(true),  // print the graphs in a file (e.g. pdf)
0067       _write(true),   // write the graphs in a root file
0068       _print_only_global(printOnlyGlobal),
0069       _make_profile_plots(makeProfilePlots),
0070       _batchMode(
0071 #ifdef DEBUG
0072           false  // false = display canvases (very time- and resource-consuming)
0073 #else
0074           true  // true = no canvases
0075 #endif
0076           ),
0077       _1dModule(true),           // cut on 1d modules
0078       _2dModule(true),           // cut on 2d modules
0079       _levelCut(DEFAULT_LEVEL),  // module level (see branch of same name)
0080       _grid_x(0),                // by default no display the grid in the canvases
0081       _grid_y(0),                // by default no display the grid in the canvases
0082       _window_width(DEFAULT_WINDOW_WIDTH),
0083       _window_height(DEFAULT_WINDOW_HEIGHT),
0084       _canvas_index(canvas_idx) {
0085 #ifdef TALKATIVE
0086   std::cout << ">>> TALKATIVE MODE ACTIVATED <<<" << std::endl;
0087 #endif
0088 #ifdef DEBUG
0089   std::cout << ">>> DEBUG MODE ACTIVATED <<<" << std::endl;
0090   std::cout << __FILE__ << ":" << __LINE__ << ":Info: inside constructor of GeometryComparisonPlotter utility"
0091             << std::endl;
0092 #endif
0093   // _canvas_index
0094   //_canvas_index = 0;
0095
0096   //_sublevel_names = {"PXB", "PXF", "TIB", "TID", "TOB", "TEC"}; // C++11
0097   _sublevel_names[0] = TString("PXB");
0098   _sublevel_names[1] = TString("PXF");
0099   _sublevel_names[2] = TString("TIB");
0100   _sublevel_names[3] = TString("TID");
0101   _sublevel_names[4] = TString("TOB");
0102   _sublevel_names[5] = TString("TEC");
0103   // TO DO: handle other structures
0104
0105   // read tree
0106   tree_file = new TFile(tree_file_name, "UPDATE");
0107   data = (TTree *)tree_file->Get("alignTree");
0108   // int branches
0109   data->SetBranchAddress("id", &branch_i["id"]);
0110   data->SetBranchAddress("inModuleList", &branch_i["inModuleList"]);
0111   data->SetBranchAddress("badModuleQuality", &branch_i["badModuleQuality"]);
0112   data->SetBranchAddress("mid", &branch_i["mid"]);
0113   data->SetBranchAddress("level", &branch_i["level"]);
0114   data->SetBranchAddress("mlevel", &branch_i["mlevel"]);
0115   data->SetBranchAddress("sublevel", &branch_i["sublevel"]);
0116   data->SetBranchAddress("useDetId", &branch_i["useDetId"]);
0117   data->SetBranchAddress("detDim", &branch_i["detDim"]);
0118   // float branches
0119   data->SetBranchAddress("x", &branch_f["x"]);
0120   data->SetBranchAddress("y", &branch_f["y"]);
0121   data->SetBranchAddress("z", &branch_f["z"]);
0122   data->SetBranchAddress("alpha", &branch_f["alpha"]);
0123   data->SetBranchAddress("beta", &branch_f["beta"]);
0124   data->SetBranchAddress("gamma", &branch_f["gamma"]);
0125   data->SetBranchAddress("phi", &branch_f["phi"]);
0126   data->SetBranchAddress("eta", &branch_f["eta"]);
0127   data->SetBranchAddress("r", &branch_f["r"]);
0128   data->SetBranchAddress("dx", &branch_f["dx"]);
0129   data->SetBranchAddress("dy", &branch_f["dy"]);
0130   data->SetBranchAddress("dz", &branch_f["dz"]);
0131   data->SetBranchAddress("dphi", &branch_f["dphi"]);
0132   data->SetBranchAddress("dr", &branch_f["dr"]);
0133   data->SetBranchAddress("dalpha", &branch_f["dalpha"]);
0134   data->SetBranchAddress("dbeta", &branch_f["dbeta"]);
0135   data->SetBranchAddress("dgamma", &branch_f["dgamma"]);
0136   if (data->GetBranch("rdphi") ==
0137       nullptr)  // in the case of rdphi branch not existing, it is created from r and dphi branches
0138   {
0139 #ifdef TALKATIVE
0140     std::cout << __FILE__ << ":" << __LINE__
0141               << ":Info: computing the rdphi branch from r and dphi branches (assuming they exist...)" << std::endl;
0142 #endif
0143     TBranch *br_rdphi = data->Branch("rdphi", &branch_f["rdphi"], "rdphi/F");
0144     for (unsigned int ientry = 0; ientry < data->GetEntries(); ientry++) {
0145       data->GetEntry(ientry);
0146       branch_f["rdphi"] = branch_f["r"] * branch_f["dphi"];
0147       br_rdphi->Fill();
0148     }
0149   } else
0150     data->SetBranchAddress("rdphi", &branch_f["rdphi"]);
0151
0152 #ifdef DEBUG
0153   std::cout << __FILE__ << ":" << __LINE__ << ":Info: branch addresses set" << std::endl;
0154 #endif
0155
0156   // style
0157   gROOT->Reset();
0158
0159   data->SetMarkerSize(0.5);
0160   data->SetMarkerStyle(6);
0161
0162   gStyle->SetOptStat("emr");
0163   gStyle->SetTitleAlign(22);
0164   gStyle->SetTitleX(0.5);
0165   gStyle->SetTitleY(0.97);
0166   gStyle->SetTitleFont(62);
0167   //gStyle->SetOptTitle(0);
0168
0169   gStyle->SetTextFont(132);
0170   gStyle->SetTextSize(0.08);
0171   gStyle->SetLabelFont(132, "x");
0172   gStyle->SetLabelFont(132, "y");
0173   gStyle->SetLabelFont(132, "z");
0174   gStyle->SetTitleSize(0.08, "x");
0175   gStyle->SetTitleSize(0.08, "y");
0176   gStyle->SetTitleSize(0.08, "z");
0177   gStyle->SetLabelSize(0.08, "x");
0178   gStyle->SetLabelSize(0.08, "y");
0179   gStyle->SetLabelSize(0.08, "z");
0180
0181   gStyle->SetMarkerStyle(8);
0182   gStyle->SetHistLineWidth(2);
0183   gStyle->SetLineStyleString(2, "[12 12]");  // postscript dashes
0184
0185   gStyle->SetFrameBorderMode(0);
0186   gStyle->SetCanvasBorderMode(0);
0187   gStyle->SetPadBorderMode(0);
0188   gStyle->SetPadColor(0);
0189   gStyle->SetCanvasColor(0);
0190   gStyle->SetTitleColor(1);
0191   gStyle->SetStatColor(0);
0192   gStyle->SetStatBorderSize(1);
0193   gStyle->SetFrameFillColor(0);
0194
0195   gStyle->SetPadTickX(1);
0196   gStyle->SetPadTickY(1);
0197
0198   gStyle->SetPadTopMargin(0.1);
0199   gStyle->SetPadRightMargin(0.05);
0200   gStyle->SetPadBottomMargin(0.16);
0201   gStyle->SetPadLeftMargin(0.18);
0202
0203 #ifdef DEBUG
0204   std::cout << __FILE__ << ":" << __LINE__ << ":Info: end of constructor" << std::endl;
0205 #endif
0206 }
0207
0208 GeometryComparisonPlotter::~GeometryComparisonPlotter() {
0209 #ifdef DEBUG
0210   std::cout << __FILE__ << ":" << __LINE__ << ":Info: in destructor of the GeometryComparisonPlotter utility"
0211             << std::endl;
0212 #endif
0213   tree_file->Close();
0214 #ifdef DEBUG
0215   std::cout << __FILE__ << ":" << __LINE__ << ":Info: ending." << std::endl;
0216 #endif
0217 }
0218
0219 // MAIN METHOD
0220 void GeometryComparisonPlotter::MakePlots(
0221     std::vector<TString> x,  // axes to combine to plot
0222     std::vector<TString> y,  // every combination (except the ones such that x=y) will be perfomed
0223     pt::ptree CFG)           // property tree storing the min and max values under <var>_min and <var>_max
0224 {
0225   /// -1) check that only existing branches are called
0226   // (we use a macro to avoid copy/paste)
0227 #define CHECK_BRANCHES(branchname_vector)                                                       \
0228   for (unsigned int i = 0; i < branchname_vector.size(); i++) {                                 \
0229     if (branch_f.find(branchname_vector[i]) == branch_f.end()) {                                \
0230       std::cout << __FILE__ << ":" << __LINE__ << ":Error: The branch " << branchname_vector[i] \
0231                 << " is not recognised." << std::endl;                                          \
0232       return;                                                                                   \
0233     }                                                                                           \
0234   }
0235   CHECK_BRANCHES(x);
0236   CHECK_BRANCHES(y);
0237
0238   const unsigned int nentries = data->GetEntries();
0239
0240 #ifdef TALKATIVE
0241   std::cout << __FILE__ << ":" << __LINE__ << ":Info: ";
0242   INSIDE_VECTOR(x);
0243   std::cout << std::endl << __FILE__ << ":" << __LINE__ << ":Info: ";
0244   INSIDE_VECTOR(y);
0245   std::cout << std::endl;
0246 #endif
0247
0248   /// 0) min and max values
0249   // the max and min of the graphs are computed from the tree if they have not been manually input yet
0250   // (we use a macro to avoid copy/paste)
0251 #define LIMITS(axes_vector)                                                          \
0252   for (unsigned int i = 0; i < axes_vector.size(); i++) {                            \
0253     if (_SF.find(axes_vector[i]) == _SF.end())                                       \
0254       _SF[axes_vector[i]] = 1.;                                                      \
0255     if (_min.find(axes_vector[i]) == _min.end())                                     \
0256       _min[axes_vector[i]] = _SF[axes_vector[i]] * data->GetMinimum(axes_vector[i]); \
0257     if (_max.find(axes_vector[i]) == _max.end())                                     \
0258       _max[axes_vector[i]] = _SF[axes_vector[i]] * data->GetMaximum(axes_vector[i]); \
0259   }
0260   LIMITS(x);
0261   LIMITS(y);
0262
0263 #ifdef TALKATIVE
0264   CHECK_MAP_CONTENT(_min, float);
0265   CHECK_MAP_CONTENT(_max, float);
0266   CHECK_MAP_CONTENT(_SF, float);
0267 #endif
0268
0269   /// 1) declare TGraphs and Histograms for profile plots if these are to be plotted
0270   // the idea is to produce at the end a table of 8 TMultiGraphs and histograms:
0271   // - 0=Tracker, with color code for the different sublevels
0272   // - 1..6=different sublevels, with color code for z < or > 0
0273   // - 7=only pixel with color code for BPIX and FPIX
0274
0275   // (convention: the six first (resp. last) correspond to z>0 (resp. z<0))
0276   // Modules with bad quality and in a list of modules that is given
0277   // by the user (e.g. list of bad/untouched modules, default: empty list)
0278   // are stored in seperate graphs and might be plotted (depends on the module
0279   // plot option, default: all modules plotted)
0280   // This means that 3*2*6 TGraphs will be filled during the loop on the TTree,
0281   // and will be arranged differently with different color codes in the TMultiGraphs
0282
0283   // For the profile plots
0284   // Either all modules, only good modules or good modules + those in a given list will be plotted
0285   // This means that 2*6 TH2F will be filled during the loop on the TTree,
0286   // and will be arranged differently with different color codes in the Histograms
0287 #ifndef NB_SUBLEVELS
0288 #define NB_SUBLEVELS 6
0289 #endif
0290 #define NB_Z_SLICES 2
0291 #define NB_MODULE_QUALITY 3
0292 #define COLOR_CODE(icolor) int(icolor / 4) + icolor + 1
0293
0294   TGraph *graphs[x.size()][y.size()][NB_SUBLEVELS * NB_Z_SLICES * NB_MODULE_QUALITY];
0295   long int ipoint[x.size()][y.size()][NB_SUBLEVELS * NB_Z_SLICES * NB_MODULE_QUALITY];
0296
0297   TMultiGraph
0298       *mgraphs[x.size()][y.size()]
0299               [2 + NB_SUBLEVELS];  // the 0th is for global plots, the 1..6th for sublevel plots, 7th for pixel only
0300   TCanvas *c[x.size()][y.size()][2 + NB_SUBLEVELS], *c_global[2 + NB_SUBLEVELS];
0301   _canvas_index++;  // this static index is a safety used in case the MakePlots method is used several times to avoid overloading
0302
0303   // histograms for profile plots,
0304   // 2D-hists to store the data
0305   // 1D-hists to calculate mean and sigma of y-values for each x-bin of the 2D-hists and for the final profile hist
0306   TH2F *histos2D[x.size()][y.size()][NB_SUBLEVELS * NB_Z_SLICES];
0307   TH1F *histos[x.size()][y.size()][NB_SUBLEVELS * NB_Z_SLICES];
0308   TH1F *histosYValues[x.size()][y.size()]
0309                      [NB_SUBLEVELS * NB_Z_SLICES];  // Used to calculate the mean and RMS for each x-bin of the 2D-hist
0310   TH1F *histosTracker
0311       [x.size()][y.size()]
0312       [NB_SUBLEVELS *
0313        NB_Z_SLICES];  // for the tracker plots all histos are copied to avoid using the same hists in different canvas
0314
0315   TCanvas *c_hist[x.size()][y.size()][2 + NB_SUBLEVELS], *c_global_hist[2 + NB_SUBLEVELS];
0316
0317   unsigned int nXBins;  // Sensible number of x-bins differs depending on the variable
0318
0319   for (unsigned int ic = 0; ic <= NB_SUBLEVELS + 1; ic++) {
0320     c_global[ic] = new TCanvas(
0321         TString::Format(
0322             "global_%s_%d", ic == 0 ? "tracker" : (ic == 7 ? "pixel" : _sublevel_names[ic - 1].Data()), _canvas_index),
0323         TString::Format("Global overview of the %s variables",
0324                         ic == 0 ? "tracker" : (ic == 7 ? "pixel" : _sublevel_names[ic - 1].Data())),
0325         _window_width,
0326         _window_height);
0327     c_global[ic]->Divide(x.size(), y.size());
0328
0329     if (_make_profile_plots) {
0330       c_global_hist[ic] =
0331           new TCanvas(TString::Format("global_profile_plots_%s_%d",
0332                                       ic == 0 ? "tracker" : (ic == 7 ? "pixel" : _sublevel_names[ic - 1].Data()),
0333                                       _canvas_index),
0334                       TString::Format("Global overview profile plots of the %s variables",
0335                                       ic == 0 ? "tracker" : (ic == 7 ? "pixel" : _sublevel_names[ic - 1].Data())),
0336                       _window_width,
0337                       _window_height);
0338       c_global_hist[ic]->Divide(x.size(), y.size());
0339     }
0340   }
0341
0342   for (unsigned int ix = 0; ix < x.size(); ix++) {
0343     for (unsigned int iy = 0; iy < y.size(); iy++) {
0344       //if (x[ix] == y[iy]) continue;       // do not plot graphs like (r,r) or (phi,phi)
0345       for (unsigned int igraph = 0; igraph < NB_SUBLEVELS * NB_Z_SLICES * NB_MODULE_QUALITY; igraph++) {
0346         // declaring
0347         ipoint[ix][iy][igraph] =
0348             0;  // the purpose of an index for every graph is to avoid thousands of points at the origin of each
0349         graphs[ix][iy][igraph] = new TGraph();
0350
0351         graphs[ix][iy][igraph]->SetMarkerColor(COLOR_CODE(igraph));
0352         graphs[ix][iy][igraph]->SetMarkerStyle(6);
0353         // pimping
0354         graphs[ix][iy][igraph]->SetName(
0355             x[ix] + y[iy] + _sublevel_names[igraph % NB_SUBLEVELS] +
0356             TString(igraph % (NB_SUBLEVELS * NB_Z_SLICES) >= NB_SUBLEVELS ? "n"
0357                                                                           : "p")  // graphs for negative/positive  z
0358             + TString(igraph >= NB_SUBLEVELS * NB_Z_SLICES ? (igraph >= 2 * NB_SUBLEVELS * NB_Z_SLICES ? "bad" : "list")
0359                                                            : "good"));  // graphs for good, bad modules and from a list
0360         graphs[ix][iy][igraph]->SetTitle(
0361             _sublevel_names[igraph % NB_SUBLEVELS] +
0362             TString(igraph % (NB_SUBLEVELS * NB_Z_SLICES) >= NB_SUBLEVELS ? " at z<0" : " at z>=0") +
0363             TString(igraph >= NB_SUBLEVELS * NB_Z_SLICES
0364                         ? (igraph >= 2 * NB_SUBLEVELS * NB_Z_SLICES ? " bad modules" : " in list")
0365                         : " good modules") +
0366             TString(";") + LateXstyle(x[ix]) + " /" + _units[x[ix]] + TString(";") + LateXstyle(y[iy]) + " /" +
0367             _units[y[iy]]);
0368         graphs[ix][iy][igraph]->SetMarkerStyle(
0369             igraph >= NB_SUBLEVELS * NB_Z_SLICES
0370                 ? (igraph >= 2 * NB_SUBLEVELS * NB_Z_SLICES ? 4 : 5)
0371                 : 6);  // empty circle for bad modules, X for those in list, dot for good ones
0372       }
0373     }
0374   }
0375
0376   // Use seperate loop for the profile histograms since we do not produce histograms for the different module qualities
0377   if (_make_profile_plots) {
0378     for (unsigned int ix = 0; ix < x.size(); ix++) {
0379       if (x[ix] == "phi")
0380         nXBins = 10;
0381       else
0382         nXBins = 40;
0383
0384       for (unsigned int iy = 0; iy < y.size(); iy++) {
0385         for (unsigned int igraph = 0; igraph < NB_SUBLEVELS * NB_Z_SLICES; igraph++) {
0386           // declaring
0387           histos2D[ix][iy][igraph] =
0388               new TH2F("2Dhist" + x[ix] + y[iy] + _sublevel_names[igraph % NB_SUBLEVELS] +
0389                            TString(igraph % (NB_SUBLEVELS * NB_Z_SLICES) >= NB_SUBLEVELS ? "n" : "p") +
0390                            std::to_string(_canvas_index),
0391                        "",
0392                        nXBins,
0393                        _min[x[ix]],
0394                        _max[x[ix]],
0395                        1000,
0396                        _min[y[iy]],
0397                        _max[y[iy]] + 1.);
0398         }
0399       }
0400     }
0401   }
0402
0403 #ifdef DEBUG
0404   std::cout << __FILE__ << ":" << __LINE__ << ":Info: Creation of the TGraph[" << x.size() << "][" << y.size() << "]["
0405             << NB_SUBLEVELS * NB_Z_SLICES * NB_MODULE_QUALITY << "] ended." << std::endl;
0406 #endif
0407
0408   /// 2) loop on the TTree data
0409 #ifdef DEBUG
0410   std::cout << __FILE__ << ":" << __LINE__ << ":Info: Looping on the TTree" << std::endl;
0411 #endif
0412 #ifdef TALKATIVE
0413   unsigned int progress = 0;
0414   std::cout << __FILE__ << ":" << __LINE__ << ":Info: 0%" << std::endl;
0415 #endif
0416   for (unsigned int ientry = 0; ientry < nentries; ientry++) {
0417 #ifdef TALKATIVE
0418     if (10 * ientry / nentries != progress) {
0419       progress = 10 * ientry / nentries;
0420       std::cout << __FILE__ << ":" << __LINE__ << ":Info: " << 10 * progress << "%" << std::endl;
0421     }
0422 #endif
0423     // load current tree entry
0424     data->GetEntry(ientry);
0425
0426     // CUTS on entry
0427     if (branch_i["level"] != _levelCut)
0428       continue;
0429     if (!_1dModule && branch_i["detDim"] == 1)
0430       continue;
0431     if (!_2dModule && branch_i["detDim"] == 2)
0432       continue;
0433
0434     // loop on the different couples of variables to plot in a graph
0435     for (unsigned int ix = 0; ix < x.size(); ix++) {
0436       // CUTS on x[ix]
0437       if (_SF[x[ix]] * branch_f[x[ix]] > _max[x[ix]] || _SF[x[ix]] * branch_f[x[ix]] < _min[x[ix]]) {
0438         //#ifdef DEBUG
0439         //                std::cout << "branch_f[x[ix]]=" << branch_f[x[ix]] << std::endl;
0440         //#endif
0441         continue;
0442       }
0443
0444       for (unsigned int iy = 0; iy < y.size(); iy++) {
0445         // CUTS on y[iy]
0446         //if (x[ix] == y[iy])                                                   continue; // TO DO: handle display when such a case occurs
0447         if (branch_i["sublevel"] < 1 || branch_i["sublevel"] > NB_SUBLEVELS)
0448           continue;
0449
0450         // FILLING histograms take even those outside the plotted range into account
0451         if (_make_profile_plots) {
0452           if (_module_plot_option == "all") {
0453             const short int igraph = (branch_i["sublevel"] - 1) + (branch_f["z"] >= 0 ? 0 : NB_SUBLEVELS);
0454             histos2D[ix][iy][igraph]->Fill(_SF[x[ix]] * branch_f[x[ix]], _SF[y[iy]] * branch_f[y[iy]]);
0455           } else if (_module_plot_option == "good" && branch_i["badModuleQuality"] == 0) {
0456             const short int igraph = (branch_i["sublevel"] - 1) + (branch_f["z"] >= 0 ? 0 : NB_SUBLEVELS);
0457             histos2D[ix][iy][igraph]->Fill(_SF[x[ix]] * branch_f[x[ix]], _SF[y[iy]] * branch_f[y[iy]]);
0458           } else if (_module_plot_option == "list" &&
0459                      (branch_i["inModuleList"] == 1 || branch_i["badModuleQuality"] == 0)) {
0460             const short int igraph = (branch_i["sublevel"] - 1) + (branch_f["z"] >= 0 ? 0 : NB_SUBLEVELS);
0461             histos2D[ix][iy][igraph]->Fill(_SF[x[ix]] * branch_f[x[ix]], _SF[y[iy]] * branch_f[y[iy]]);
0462           }
0463         }
0464
0465         // restrict scatter plots to chosen range
0466         if (_SF[y[iy]] * branch_f[y[iy]] > _max[y[iy]] || _SF[y[iy]] * branch_f[y[iy]] < _min[y[iy]]) {
0467           //#ifdef DEBUG
0468           //                    std::cout << "branch_f[y[iy]]=" << branch_f[y[iy]] << std::endl;
0469           //#endif
0470           continue;
0471         }
0472
0473         // FILLING GRAPH
0474         if (y.size() >= x.size()) {
0475           if (branch_i["inModuleList"] == 0 && branch_i["badModuleQuality"] == 0) {
0476             const short int igraph = (branch_i["sublevel"] - 1) + (branch_f["z"] >= 0 ? 0 : NB_SUBLEVELS);
0477             graphs[ix][iy][igraph]->SetPoint(
0478                 ipoint[ix][iy][igraph], _SF[x[ix]] * branch_f[x[ix]], _SF[y[iy]] * branch_f[y[iy]]);
0479             ipoint[ix][iy][igraph]++;
0480           }
0481           if (branch_i["inModuleList"] > 0) {
0482             const short int igraph =
0483                 (branch_i["sublevel"] - 1) + (branch_f["z"] >= 0 ? 0 : NB_SUBLEVELS) + NB_SUBLEVELS * NB_Z_SLICES;
0484             graphs[ix][iy][igraph]->SetPoint(
0485                 ipoint[ix][iy][igraph], _SF[x[ix]] * branch_f[x[ix]], _SF[y[iy]] * branch_f[y[iy]]);
0486             ipoint[ix][iy][igraph]++;
0487           }
0488           if (branch_i["badModuleQuality"] > 0) {
0489             const short int igraph =
0490                 (branch_i["sublevel"] - 1) + (branch_f["z"] >= 0 ? 0 : NB_SUBLEVELS) + 2 * NB_SUBLEVELS * NB_Z_SLICES;
0491             graphs[ix][iy][igraph]->SetPoint(
0492                 ipoint[ix][iy][igraph], _SF[x[ix]] * branch_f[x[ix]], _SF[y[iy]] * branch_f[y[iy]]);
0493             ipoint[ix][iy][igraph]++;
0494           }
0495         } else {
0496           if (branch_i["inModuleList"] == 0 && branch_i["badModuleQuality"] == 0) {
0497             const short int igraph = (branch_i["sublevel"] - 1) + (branch_f["z"] >= 0 ? 0 : NB_SUBLEVELS);
0498             graphs[iy][ix][igraph]->SetPoint(
0499                 ipoint[iy][ix][igraph], _SF[x[ix]] * branch_f[x[ix]], _SF[y[iy]] * branch_f[y[iy]]);
0500             ipoint[iy][ix][igraph]++;
0501           }
0502           if (branch_i["inModuleList"] > 0) {
0503             const short int igraph =
0504                 (branch_i["sublevel"] - 1) + (branch_f["z"] >= 0 ? 0 : NB_SUBLEVELS) + NB_SUBLEVELS * NB_Z_SLICES;
0505             graphs[iy][ix][igraph]->SetPoint(
0506                 ipoint[iy][ix][igraph], _SF[x[ix]] * branch_f[x[ix]], _SF[y[iy]] * branch_f[y[iy]]);
0507             ipoint[iy][ix][igraph]++;
0508           }
0509           if (branch_i["badModuleQuality"] > 0) {
0510             const short int igraph =
0511                 (branch_i["sublevel"] - 1) + (branch_f["z"] >= 0 ? 0 : NB_SUBLEVELS) + 2 * NB_SUBLEVELS * NB_Z_SLICES;
0512             graphs[iy][ix][igraph]->SetPoint(
0513                 ipoint[ix][iy][igraph], _SF[x[ix]] * branch_f[x[ix]], _SF[y[iy]] * branch_f[y[iy]]);
0514             ipoint[iy][ix][igraph]++;
0515           }
0516         }
0517       }
0518     }
0519   }
0520 #ifdef TALKATIVE
0521   std::cout << __FILE__ << ":" << __LINE__ << ":Info: 100%\tLoop ended" << std::endl;
0522 #endif
0523
0524   /// 3) merge TGraph objects into TMultiGraph objects, then draw, print and write (according to the options _batchMode, _print and _write respectively)
0525   gROOT->SetBatch(_batchMode);  // if true, then equivalent to "root -b", i.e. no canvas
0526   if (_write) {                 // opening the file to write the graphs
0527     output = new TFile(_output_directory + TString(_output_filename),
0528                        "UPDATE");  // possibly existing file will be updated, otherwise created
0529     if (output->IsZombie()) {
0530       std::cout << __FILE__ << ":" << __LINE__ << ":Error: Opening of " << _output_directory + TString(_output_filename)
0531                 << " failed" << std::endl;
0532       exit(-1);
0533     }
0534 #ifdef TALKATIVE
0535     std::cout << __FILE__ << ":" << __LINE__ << ":Info: output file is "
0536               << _output_directory + TString(_output_filename) << std::endl;
0537 #endif
0538   }
0539   // declaring TMultiGraphs and TCanvas
0540   // Usually more y variables than x variables
0541   // creating TLegend
0542   TLegend *legend = MakeLegend(.1, .92, .9, 1., NB_SUBLEVELS);
0543   if (_write)
0544     legend->Write();
0545
0546   // check which modules are supposed to be plotted
0547   unsigned int n_module_types = 1;
0548   if (_module_plot_option == "all") {
0549     n_module_types = 3;  //plot all modules (good, list and bad )
0550   } else if (_module_plot_option == "list") {
0551     n_module_types = 2;  // plot good modules and those in the list
0552   } else if (_module_plot_option == "good") {
0553     n_module_types = 1;  // only plot the modules that are neither bad or in the list
0554   }
0555
0556 #define INDEX_IN_GLOBAL_CANVAS(i1, i2) 1 + i1 + i2 *x.size()
0557   // running on the TGraphs to produce the TMultiGraph and draw/print them
0558   for (unsigned int ix = 0; ix < x.size(); ix++) {
0559 #ifdef DEBUG
0560     std::cout << __FILE__ << ":" << __LINE__ << ":Info: x[" << ix << "]=" << x[ix] << std::endl;
0561 #endif
0562
0563     // looping on Y axes
0564     for (unsigned int iy = 0; iy < y.size(); iy++) {
0565       TString min_branch = TString::Format("%s_min", y[iy].Data());
0566       TString max_branch = TString::Format("%s_max", y[iy].Data());
0567
0568 #ifdef DEBUG
0569       std::cout << __FILE__ << ":" << __LINE__ << ":Info: x[" << ix << "]=" << x[ix] << " and y[" << iy << "]=" << y[iy]
0570                 << "\t-> creating TMultiGraph" << std::endl;
0571 #endif
0572       mgraphs[ix][iy][0] = new TMultiGraph(
0573           TString::Format("mgr_%s_vs_%s_tracker_%d",
0574                           x[ix].Data(),
0575                           y[iy].Data(),
0576                           _canvas_index),  // name
0577           //LateXstyle(x[ix]) + TString(" vs. ") + LateXstyle(y[iy]) + TString(" for Tracker") // graph title
0578           TString(";") + LateXstyle(x[ix]) + " /" + _units[x[ix]]          // x axis title
0579               + TString(";") + LateXstyle(y[iy]) + " /" + _units[y[iy]]);  // y axis title
0580
0581       mgraphs[ix][iy][7] = new TMultiGraph(
0582           TString::Format("mgr_%s_vs_%s_pixel_%d",
0583                           x[ix].Data(),
0584                           y[iy].Data(),
0585                           _canvas_index),  // name
0586           //LateXstyle(x[ix]) + TString(" vs. ") + LateXstyle(y[iy]) + TString(" for Tracker") // graph title
0587           TString(";") + LateXstyle(x[ix]) + " /" + _units[x[ix]]          // x axis title
0588               + TString(";") + LateXstyle(y[iy]) + " /" + _units[y[iy]]);  // y axis title
0589
0590       /// TRACKER and PIXEL
0591       // fixing ranges and filling TMultiGraph
0592       // for (unsigned short int jgraph = NB_SUBLEVELS*NB_Z_SLICES-1 ; jgraph >= 0 ; --jgraph)
0593       for (unsigned short int jgraph = 0; jgraph < NB_SUBLEVELS * NB_Z_SLICES * n_module_types; jgraph++) {
0594         unsigned short int igraph =
0595             NB_SUBLEVELS * NB_Z_SLICES * n_module_types - jgraph -
0596             1;  // reverse counting for humane readability (one of the sublevel takes much more place than the others)
0597
0598 #ifdef DEBUG
0599         std::cout << __FILE__ << ":" << __LINE__ << ":Info: writing TGraph to file" << std::endl;
0600 #endif
0601         // write into root file
0602         if (_write)
0603           graphs[ix][iy][igraph]->Write();
0604         if (graphs[ix][iy][igraph]->GetN() == 0) {
0605 #ifdef TALKATIVE
0606           std::cout << __FILE__ << ":" << __LINE__ << ":Info: " << graphs[ix][iy][igraph]->GetName() << " is empty."
0607                     << std::endl;
0608 #endif
0609           continue;
0610         }
0611 #ifdef DEBUG
0612         std::cout << __FILE__ << ":" << __LINE__ << ":Info: cloning, coloring and adding TGraph "
0613                   << _sublevel_names[igraph % NB_SUBLEVELS] << (igraph >= NB_SUBLEVELS ? "(z<0)" : "(z>0)")
0614                   << " to global TMultiGraph" << std::endl;
0615 #endif
0616         // clone to prevent any injure on the graph
0617         TGraph *gr = (TGraph *)graphs[ix][iy][igraph]->Clone();
0618         // color
0619         gr->SetMarkerColor(COLOR_CODE(igraph % NB_SUBLEVELS));
0620         mgraphs[ix][iy][0]->Add(gr, "P");  //, (mgraphs[ix][iy][0]->GetListOfGraphs()==0?"AP":"P"));
0621
0622         if (igraph % NB_SUBLEVELS == 0 || igraph % NB_SUBLEVELS == 1)
0623           mgraphs[ix][iy][7]->Add(gr, "P");  // Add BPIX (0) and FPIX (1) to pixel plot
0624       }
0625
0626       /// SUBLEVELS (1..6)
0627       for (unsigned int isublevel = 1; isublevel <= NB_SUBLEVELS; isublevel++) {
0628 #ifdef DEBUG
0629         std::cout << __FILE__ << ":" << __LINE__ << ":Info: cloning, coloring and adding TGraph "
0630                   << _sublevel_names[isublevel - 1] << " to sublevel TMultiGraph" << std::endl;
0631 #endif
0632         mgraphs[ix][iy][isublevel] =
0633             new TMultiGraph(TString::Format("%s_vs_%s_%s_%d",
0634                                             x[ix].Data(),
0635                                             y[iy].Data(),
0636                                             _sublevel_names[isublevel - 1].Data(),
0637                                             _canvas_index),  // name
0638                             //LateXstyle(x[ix]) + TString(" vs. ") + LateXstyle(y[iy]) + TString(" for ") +
0639                             _sublevel_names[isublevel - 1]                                   // graph title
0640                                 + TString(";") + LateXstyle(x[ix]) + " /" + _units[x[ix]]    // x axis title
0641                                 + TString(";") + LateXstyle(y[iy]) + " /" + _units[y[iy]]);  // y axis title
0642
0643         graphs[ix][iy][isublevel - 1]->SetMarkerColor(kBlack);
0644         graphs[ix][iy][NB_SUBLEVELS + isublevel - 1]->SetMarkerColor(kRed);
0645         graphs[ix][iy][2 * NB_SUBLEVELS + isublevel - 1]->SetMarkerColor(kGray + 1);
0646         graphs[ix][iy][3 * NB_SUBLEVELS + isublevel - 1]->SetMarkerColor(kRed - 7);
0647         graphs[ix][iy][4 * NB_SUBLEVELS + isublevel - 1]->SetMarkerColor(kGray + 1);
0648         graphs[ix][iy][5 * NB_SUBLEVELS + isublevel - 1]->SetMarkerColor(kRed - 7);
0649         if (graphs[ix][iy][isublevel - 1]->GetN() > 0)
0650           mgraphs[ix][iy][isublevel]->Add(
0651               graphs[ix][iy][isublevel - 1],
0652               "P");  //(mgraphs[ix][iy][isublevel-1]->GetListOfGraphs()==0?"AP":"P")); // z>0
0653 #ifdef TALKATIVE
0654         else
0655           std::cout << __FILE__ << ":" << __LINE__
0656                     << ":Info: graphs[ix][iy][isublevel-1]=" << graphs[ix][iy][isublevel - 1]->GetName()
0657                     << " is empty -> not added into " << mgraphs[ix][iy][isublevel]->GetName() << std::endl;
0658 #endif
0659         if (graphs[ix][iy][NB_SUBLEVELS + isublevel - 1]->GetN() > 0)
0660           mgraphs[ix][iy][isublevel]->Add(
0661               graphs[ix][iy][NB_SUBLEVELS + isublevel - 1],
0662               "P");  //(mgraphs[ix][iy][isublevel-1]->GetListOfGraphs()==0?"AP":"P")); // z<0
0663 #ifdef TALKATIVE
0664         else
0665           std::cout << __FILE__ << ":" << __LINE__ << ":Info: graphs[ix][iy][NB_SUBLEVEL+isublevel-1]="
0666                     << graphs[ix][iy][NB_Z_SLICES + isublevel - 1]->GetName() << " is empty -> not added into "
0667                     << mgraphs[ix][iy][isublevel]->GetName() << std::endl;
0668 #endif
0669 #if NB_Z_SLICES != 2
0670         std::cout << __FILE__ << ":" << __LINE__ << ":Error: color code incomplete for Z slices..." << std::endl;
0671 #endif
0672         if (_module_plot_option == "all") {
0673           if (graphs[ix][iy][2 * NB_SUBLEVELS + isublevel - 1]->GetN() > 0)
0674             mgraphs[ix][iy][isublevel]->Add(graphs[ix][iy][2 * NB_SUBLEVELS + isublevel - 1], "P");
0675           if (graphs[ix][iy][3 * NB_SUBLEVELS + isublevel - 1]->GetN() > 0)
0676             mgraphs[ix][iy][isublevel]->Add(graphs[ix][iy][3 * NB_SUBLEVELS + isublevel - 1], "P");
0677           if (graphs[ix][iy][4 * NB_SUBLEVELS + isublevel - 1]->GetN() > 0)
0678             mgraphs[ix][iy][isublevel]->Add(graphs[ix][iy][4 * NB_SUBLEVELS + isublevel - 1], "P");
0679           if (graphs[ix][iy][5 * NB_SUBLEVELS + isublevel - 1]->GetN() > 0)
0680             mgraphs[ix][iy][isublevel]->Add(graphs[ix][iy][5 * NB_SUBLEVELS + isublevel - 1], "P");
0681         }
0682         if (_module_plot_option == "list") {
0683           if (graphs[ix][iy][2 * NB_SUBLEVELS + isublevel - 1]->GetN() > 0)
0684             mgraphs[ix][iy][isublevel]->Add(graphs[ix][iy][2 * NB_SUBLEVELS + isublevel - 1], "P");
0685           if (graphs[ix][iy][3 * NB_SUBLEVELS + isublevel - 1]->GetN() > 0)
0686             mgraphs[ix][iy][isublevel]->Add(graphs[ix][iy][3 * NB_SUBLEVELS + isublevel - 1], "P");
0687         }
0688       }
0689
0690       // fixing ranges, saving, and drawing of TMultiGraph (tracker AND sublevels AND pixel, i.e. 2+NB_SUBLEVELS objects)
0691       // the individual canvases are saved, but the global are just drawn and will be saved later
0692       for (unsigned short int imgr = 0; imgr <= NB_SUBLEVELS + 1; imgr++) {
0693 #ifdef DEBUG
0694         std::cout << __FILE__ << ":" << __LINE__ << ":Info: treating individual canvases." << std::endl;
0695 #endif
0696         // drawing into individual canvas and printing it (including a legend for the tracker canvas)
0697         c[ix][iy][imgr] = new TCanvas(
0698             TString::Format("c_%s_vs_%s_%s_%d",
0699                             x[ix].Data(),
0700                             y[iy].Data(),
0701                             imgr == 0 ? "tracker" : (imgr == 7 ? "pixel" : _sublevel_names[imgr - 1].Data()),
0702                             _canvas_index),
0703             TString::Format("%s vs. %s at %s level",
0704                             x[ix].Data(),
0705                             y[iy].Data(),
0706                             imgr == 0 ? "tracker" : (imgr == 7 ? "pixel" : _sublevel_names[imgr - 1].Data())),
0707             _window_width,
0708             _window_height);
0709         c[ix][iy][imgr]->SetGrid(_grid_x, _grid_y);  // grid
0710
0711         if (mgraphs[ix][iy][imgr]->GetListOfGraphs() != nullptr) {
0712           if (CFG.count(min_branch.Data()) > 0) {
0713             mgraphs[ix][iy][imgr]->SetMinimum(CFG.get<float>(min_branch.Data()));
0714           }
0715           if (CFG.count(max_branch.Data()) > 0) {
0716             mgraphs[ix][iy][imgr]->SetMaximum(CFG.get<float>(max_branch.Data()));
0717           }
0718           mgraphs[ix][iy][imgr]->Draw("A");
0719         }
0720         if (imgr == 0 && _legend)
0721           legend->Draw();  // only for the tracker
0722         if (_print && !_print_only_global)
0723           c[ix][iy][imgr]->Print(
0724               _output_directory + mgraphs[ix][iy][imgr]->GetName() + ExtensionFromPrintOption(_print_option),
0725               _print_option);
0726
0727         // writing into root file
0728         if (_write)
0729           mgraphs[ix][iy][imgr]->Write();
0730
0731         // drawing into global canvas
0732         c_global[imgr]->cd(INDEX_IN_GLOBAL_CANVAS(ix, iy));
0733         c_global[imgr]->GetPad(INDEX_IN_GLOBAL_CANVAS(ix, iy))->SetFillStyle(4000);         //  make the pad transparent
0734         c_global[imgr]->GetPad(INDEX_IN_GLOBAL_CANVAS(ix, iy))->SetGrid(_grid_x, _grid_y);  // grid
0735         if (mgraphs[ix][iy][imgr]->GetListOfGraphs() != nullptr) {
0736           if (CFG.count(min_branch.Data()) > 0) {
0737             mgraphs[ix][iy][imgr]->SetMinimum(CFG.get<float>(min_branch.Data()));
0738           }
0739           if (CFG.count(max_branch.Data()) > 0) {
0740             mgraphs[ix][iy][imgr]->SetMaximum(CFG.get<float>(max_branch.Data()));
0741           }
0742           mgraphs[ix][iy][imgr]->Draw("A");
0743         }
0744         // printing will be performed after customisation (e.g. legend or title) just after the loops on ix and iy
0745       }
0746     }  // end of loop on y
0747   }    // end of loop on x
0748
0749   // CUSTOMISATION
0750   gStyle->SetOptTitle(0);  // otherwise, the title is repeated in every pad of the global canvases
0751                            // -> instead, we will write it in the upper part in a TPaveText or in a TLegend
0752   for (unsigned int ic = 0; ic <= NB_SUBLEVELS + 1; ic++) {
0753     c_global[ic]->Draw();
0754
0755     // setting legend to tracker canvases
0756     if (!_legend)
0757       break;
0758     TCanvas *c_temp = (TCanvas *)c_global[ic]->Clone(c_global[ic]->GetTitle() + TString("_sub"));
0759     c_temp->Draw();
0760     c_global[ic] = new TCanvas(
0761         c_temp->GetName() + TString("_final"), c_temp->GetTitle(), c_temp->GetWindowWidth(), c_temp->GetWindowHeight());
0762     c_global[ic]->Draw();
0763     TPad *p_up = new TPad(TString("legend_") + c_temp->GetName(),
0764                           "",
0765                           0.,
0766                           0.9,
0767                           1.,
0768                           1.,  // relative position
0769                           -1,
0770                           0,
0771                           0),  // display options
0772         *p_down = new TPad(TString("main_") + c_temp->GetName(), "", 0., 0., 1., 0.9, -1, 0, 0);
0773     // in the lower part, draw the plots
0774     p_down->Draw();
0775     p_down->cd();
0776     c_temp->DrawClonePad();
0777     c_global[ic]->cd();
0778     // in the upper part, pimp the canvas :p
0779     p_up->Draw();
0780     p_up->cd();
0781     if (ic == 0)  // tracker
0782     {
0783       TLegend *global_legend = MakeLegend(.05, .1, .7, .8, NB_SUBLEVELS);  //, "brNDC");
0784       global_legend->Draw();
0785       TPaveText *pt_geom = new TPaveText(.75, .1, .95, .8, "NB");
0786       pt_geom->SetFillColor(0);
0787       pt_geom->SetTextSize(0.25);
0788       pt_geom->AddText(TString("x: ") + _reference_name);
0789       pt_geom->AddText(TString("y: ") + _alignment_name + TString(" - ") + _reference_name);
0790       pt_geom->Draw();
0791     } else if (ic == 7)  // pixel
0792     {
0793       TLegend *global_legend = MakeLegend(.05, .1, .7, .8, 2);  //, "brNDC");
0794       global_legend->Draw();
0795       TPaveText *pt_geom = new TPaveText(.75, .1, .95, .8, "NB");
0796       pt_geom->SetFillColor(0);
0797       pt_geom->SetTextSize(0.25);
0798       pt_geom->AddText(TString("x: ") + _reference_name);
0799       pt_geom->AddText(TString("y: ") + _alignment_name + TString(" - ") + _reference_name);
0800       pt_geom->Draw();
0801     } else  // sublevels
0802     {
0803       TPaveText *pt = new TPaveText(.05, .1, .7, .8, "NB");
0804       pt->SetFillColor(0);
0805       pt->AddText(_sublevel_names[ic - 1]);
0806       pt->Draw();
0807       TPaveText *pt_geom = new TPaveText(.6, .1, .95, .8, "NB");
0808       pt_geom->SetFillColor(0);
0809       pt_geom->SetTextSize(0.3);
0810       pt_geom->AddText(TString("x: ") + _reference_name);
0811       pt_geom->AddText(TString("y: ") + _alignment_name + TString(" - ") + _reference_name);
0812       pt_geom->Draw();
0813     }
0814     // printing
0815     if (_print)
0816       c_global[ic]->Print(_output_directory + c_global[ic]->GetName() + ExtensionFromPrintOption(_print_option),
0817                           _print_option);
0818     if (_write)
0819       c_global[ic]->Write();
0820   }
0821
0822   // printing global canvases
0823   if (_write)
0824     output->Close();
0825
0826   // Now produce the profile plots if the option is chosen
0827   // Use seperate loops since no seperate plots are produced for different module qualities
0828   if (_make_profile_plots) {
0829     // Fill Content of 2D-hists into 1D-hists for the profile plots
0830     // Loop over all y-bins for a certain x-bin, calculate mean and RMS as entries of the 1D-hists
0831     bool entries = false;
0832     for (unsigned int ix = 0; ix < x.size(); ix++) {
0833       for (unsigned int iy = 0; iy < y.size(); iy++) {
0834         TString min_branch = TString::Format("%s_min", y[iy].Data());
0835         TString max_branch = TString::Format("%s_max", y[iy].Data());
0836         for (unsigned int igraph = 0; igraph < NB_SUBLEVELS * NB_Z_SLICES; igraph++) {
0837           // Declare hists which will be plotted for the profile plots
0838           histos[ix][iy][igraph] =
0839               new TH1F("1Dhist" + x[ix] + y[iy] + _sublevel_names[igraph % NB_SUBLEVELS] +
0840                            TString(igraph % (NB_SUBLEVELS * NB_Z_SLICES) >= NB_SUBLEVELS ? "n" : "p") +
0841                            std::to_string(_canvas_index),
0842                        "",
0843                        histos2D[ix][iy][igraph]->GetXaxis()->GetNbins(),
0844                        _min[x[ix]],
0845                        _max[x[ix]]);
0846           histos[ix][iy][igraph]->SetMarkerColor(COLOR_CODE(igraph));
0847           histos[ix][iy][igraph]->SetLineColor(COLOR_CODE(igraph));
0848           histos[ix][iy][igraph]->StatOverflows(kTRUE);
0849
0850           // Loop over x bins
0851           for (int binx = 0; binx <= histos2D[ix][iy][igraph]->GetXaxis()->GetNbins(); binx++) {
0852             entries = false;
0853             // Declare y-histogram for each x bin
0854             histosYValues[ix][iy][igraph] =
0855                 new TH1F("1Dhist_Y-Values" + x[ix] + y[iy] + _sublevel_names[igraph % NB_SUBLEVELS] +
0856                              TString(igraph % (NB_SUBLEVELS * NB_Z_SLICES) >= NB_SUBLEVELS ? "n" : "p") +
0857                              std::to_string(_canvas_index) + std::to_string(binx),
0858                          "",
0859                          histos2D[ix][iy][igraph]->GetYaxis()->GetNbins(),
0860                          _min[y[iy]],
0861                          _max[y[iy]] + 1.);
0862             histosYValues[ix][iy][igraph]->StatOverflows(kTRUE);
0863             // Loop over y-bins for each x-bin of the 2D histogram and put it into the 1-d y histograms
0864             // Take overflow bin into account
0865             for (int biny = 0; biny <= histos2D[ix][iy][igraph]->GetYaxis()->GetNbins() + 1; biny++) {
0866               if (histos2D[ix][iy][igraph]->GetBinContent(binx, biny) > 1.) {
0867                 histosYValues[ix][iy][igraph]->SetBinContent(biny, histos2D[ix][iy][igraph]->GetBinContent(binx, biny));
0868                 entries = true;
0869               }
0870             }
0871             if (entries) {
0872               histos[ix][iy][igraph]->SetBinContent(binx, histosYValues[ix][iy][igraph]->GetMean());
0873               histos[ix][iy][igraph]->SetBinError(binx, histosYValues[ix][iy][igraph]->GetRMS());
0874             } else
0875               histos[ix][iy][igraph]->SetBinContent(binx, -999999.);
0876           }
0877         }
0878
0879         // Customize and print the histograms
0880
0881         /// TRACKER
0882         // fixing ranges and draw profile plot histos
0883
0884         c_hist[ix][iy][0] =
0885             new TCanvas(TString::Format("c_hist_%s_vs_%s_tracker_%d", x[ix].Data(), y[iy].Data(), _canvas_index),
0886                         TString::Format("Profile plot %s vs. %s at tracker level", x[ix].Data(), y[iy].Data()),
0887                         _window_width,
0888                         _window_height);
0889         c_hist[ix][iy][0]->SetGrid(_grid_x, _grid_y);  // grid
0890         // Draw the frame that will contain the histograms
0891         // One needs to specify the binning and title
0892         c_hist[ix][iy][0]->GetPad(0)->DrawFrame(
0893             _min[x[ix]],
0894             CFG.count(min_branch.Data()) > 0 ? CFG.get<float>(min_branch.Data()) : _min[y[iy]],
0895             _max[x[ix]],
0896             CFG.count(max_branch.Data()) > 0 ? CFG.get<float>(max_branch.Data()) : _max[y[iy]],
0897             TString(";") + LateXstyle(x[ix]) + " /" + _units[x[ix]] + TString(";") + LateXstyle(y[iy]) + " /" +
0898                 _units[y[iy]]);
0899         if (_legend)
0900           legend->Draw("same");
0901
0902         for (unsigned short int jgraph = 0; jgraph < NB_SUBLEVELS * NB_Z_SLICES; jgraph++) {
0903           unsigned short int igraph =
0904               NB_SUBLEVELS * NB_Z_SLICES - jgraph -
0905               1;  // reverse counting for humane readability (one of the sublevel takes much more place than the others)
0906
0907           // clone to prevent any injure on the graph
0908           histosTracker[ix][iy][igraph] = (TH1F *)histos[ix][iy][igraph]->Clone();
0909           // color
0910           histosTracker[ix][iy][igraph]->SetMarkerColor(COLOR_CODE(igraph % NB_SUBLEVELS));
0911           histosTracker[ix][iy][igraph]->SetLineColor(COLOR_CODE(igraph % NB_SUBLEVELS));
0912           histosTracker[ix][iy][igraph]->SetMarkerStyle(6);
0913           histosTracker[ix][iy][igraph]->Draw("same pe0");
0914         }
0915
0916         if (_print && !_print_only_global)
0917           c_hist[ix][iy][0]->Print(
0918               _output_directory +
0919                   TString::Format("Profile_plot_%s_vs_%s_tracker_%d", x[ix].Data(), y[iy].Data(), _canvas_index) +
0920                   ExtensionFromPrintOption(_print_option),
0921               _print_option);
0922
0923         //Draw into profile hists global tracker canvas
0924         c_global_hist[0]->cd(INDEX_IN_GLOBAL_CANVAS(ix, iy));
0925         c_global_hist[0]->GetPad(INDEX_IN_GLOBAL_CANVAS(ix, iy))->SetFillStyle(4000);  //  make the pad transparent
0926         c_global_hist[0]->GetPad(INDEX_IN_GLOBAL_CANVAS(ix, iy))->SetGrid(_grid_x, _grid_y);  // grid
0927         c_global_hist[0]
0928             ->GetPad(INDEX_IN_GLOBAL_CANVAS(ix, iy))
0929             ->DrawFrame(_min[x[ix]],
0930                         CFG.count(min_branch.Data()) > 0 ? CFG.get<float>(min_branch.Data()) : _min[y[iy]],
0931                         _max[x[ix]],
0932                         CFG.count(max_branch.Data()) > 0 ? CFG.get<float>(max_branch.Data()) : _max[y[iy]],
0933                         TString(";") + LateXstyle(x[ix]) + " /" + _units[x[ix]] + TString(";") + LateXstyle(y[iy]) +
0934                             " /" + _units[y[iy]]);
0935
0936         for (unsigned short int jgraph = 0; jgraph < NB_SUBLEVELS * NB_Z_SLICES; jgraph++) {
0937           unsigned short int igraph =
0938               NB_SUBLEVELS * NB_Z_SLICES - jgraph -
0939               1;  // reverse counting for humane readability (one of the sublevel takes much more place than the others)
0940           histosTracker[ix][iy][igraph]->Draw("same pe0");
0941         }
0942
0943         /// PIXEL
0944         // fixing ranges and draw profile plot histos
0945
0946         c_hist[ix][iy][7] =
0947             new TCanvas(TString::Format("c_hist_%s_vs_%s_pixel_%d", x[ix].Data(), y[iy].Data(), _canvas_index),
0948                         TString::Format("Profile plot %s vs. %s at pixel level", x[ix].Data(), y[iy].Data()),
0949                         _window_width,
0950                         _window_height);
0951         c_hist[ix][iy][7]->SetGrid(_grid_x, _grid_y);  // grid
0952         // Draw the frame that will contain the histograms
0953         // One needs to specify the binning and title
0954         c_hist[ix][iy][7]->GetPad(0)->DrawFrame(
0955             _min[x[ix]],
0956             CFG.count(min_branch.Data()) > 0 ? CFG.get<float>(min_branch.Data()) : _min[y[iy]],
0957             _max[x[ix]],
0958             CFG.count(max_branch.Data()) > 0 ? CFG.get<float>(max_branch.Data()) : _max[y[iy]],
0959             TString(";") + LateXstyle(x[ix]) + " /" + _units[x[ix]] + TString(";") + LateXstyle(y[iy]) + " /" +
0960                 _units[y[iy]]);
0961         if (_legend)
0962           legend->Draw("same");
0963
0964         for (unsigned short int jgraph = 0; jgraph < NB_SUBLEVELS * NB_Z_SLICES; jgraph++) {
0965           unsigned short int igraph =
0966               NB_SUBLEVELS * NB_Z_SLICES - jgraph -
0967               1;  // reverse counting for humane readability (one of the sublevel takes much more place than the others)
0968
0969           if (igraph % NB_SUBLEVELS == 0 || igraph % NB_SUBLEVELS == 1)  //Only BPIX and FPIX
0970           {
0971             // clone to prevent any injure on the graph
0972             histosTracker[ix][iy][igraph] = (TH1F *)histos[ix][iy][igraph]->Clone();
0973             // color
0974             histosTracker[ix][iy][igraph]->SetMarkerColor(COLOR_CODE(igraph % NB_SUBLEVELS));
0975             histosTracker[ix][iy][igraph]->SetLineColor(COLOR_CODE(igraph % NB_SUBLEVELS));
0976             histosTracker[ix][iy][igraph]->SetMarkerStyle(6);
0977             histosTracker[ix][iy][igraph]->Draw("same pe0");
0978           }
0979         }
0980
0981         if (_print && !_print_only_global)
0982           c_hist[ix][iy][7]->Print(
0983               _output_directory +
0984                   TString::Format("Profile_plot_%s_vs_%s_pixel_%d", x[ix].Data(), y[iy].Data(), _canvas_index) +
0985                   ExtensionFromPrintOption(_print_option),
0986               _print_option);
0987
0988         //Draw into profile hists global tracker canvas
0989         c_global_hist[7]->cd(INDEX_IN_GLOBAL_CANVAS(ix, iy));
0990         c_global_hist[7]->GetPad(INDEX_IN_GLOBAL_CANVAS(ix, iy))->SetFillStyle(4000);  //  make the pad transparent
0991         c_global_hist[7]->GetPad(INDEX_IN_GLOBAL_CANVAS(ix, iy))->SetGrid(_grid_x, _grid_y);  // grid
0992         c_global_hist[7]
0993             ->GetPad(INDEX_IN_GLOBAL_CANVAS(ix, iy))
0994             ->DrawFrame(_min[x[ix]],
0995                         CFG.count(min_branch.Data()) > 0 ? CFG.get<float>(min_branch.Data()) : _min[y[iy]],
0996                         _max[x[ix]],
0997                         CFG.count(max_branch.Data()) > 0 ? CFG.get<float>(max_branch.Data()) : _max[y[iy]],
0998                         TString(";") + LateXstyle(x[ix]) + " /" + _units[x[ix]] + TString(";") + LateXstyle(y[iy]) +
0999                             " /" + _units[y[iy]]);
1000
1001         for (unsigned short int jgraph = 0; jgraph < NB_SUBLEVELS * NB_Z_SLICES; jgraph++) {
1002           unsigned short int igraph =
1003               NB_SUBLEVELS * NB_Z_SLICES - jgraph -
1004               1;  // reverse counting for humane readability (one of the sublevel takes much more place than the others)
1005           histosTracker[ix][iy][igraph]->Draw("same pe0");
1006         }
1007         // printing will be performed after customisation (e.g. legend or title) just after the loops on ix and iy
1008         /// SUBLEVELS (1..6)
1009         for (unsigned int isublevel = 1; isublevel <= NB_SUBLEVELS; isublevel++) {
1010           // Draw and print profile histograms
1011           c_hist[ix][iy][isublevel] =
1012               new TCanvas(TString::Format("c_hist_%s_vs_%s_%s_%d",
1013                                           x[ix].Data(),
1014                                           y[iy].Data(),
1015                                           isublevel == 0 ? "tracker" : _sublevel_names[isublevel - 1].Data(),
1016                                           _canvas_index),
1017                           TString::Format("Profile plot %s vs. %s at %s level",
1018                                           x[ix].Data(),
1019                                           y[iy].Data(),
1020                                           isublevel == 0 ? "tracker" : _sublevel_names[isublevel - 1].Data()),
1021                           _window_width,
1022                           _window_height);
1023           c_hist[ix][iy][isublevel]->SetGrid(_grid_x, _grid_y);  // grid
1024           c_hist[ix][iy][isublevel]->GetPad(0)->DrawFrame(
1025               _min[x[ix]],
1026               CFG.count(min_branch.Data()) > 0 ? CFG.get<float>(min_branch.Data()) : _min[y[iy]],
1027               _max[x[ix]],
1028               CFG.count(max_branch.Data()) > 0 ? CFG.get<float>(max_branch.Data()) : _max[y[iy]],
1029               TString(";") + LateXstyle(x[ix]) + " /" + _units[x[ix]] + TString(";") + LateXstyle(y[iy]) + " /" +
1030                   _units[y[iy]]);
1031
1032           histos[ix][iy][isublevel - 1]->SetMarkerColor(kBlack);
1033           histos[ix][iy][isublevel - 1]->SetLineColor(kBlack);
1034           histos[ix][iy][NB_SUBLEVELS + isublevel - 1]->SetMarkerColor(kRed);
1035           histos[ix][iy][NB_SUBLEVELS + isublevel - 1]->SetLineColor(kRed);
1036
1037           histos[ix][iy][isublevel - 1]->Draw("same pe0");
1038           histos[ix][iy][NB_SUBLEVELS + isublevel - 1]->Draw("same pe0");
1039
1040           if (_print && !_print_only_global)
1041             c_hist[ix][iy][isublevel]->Print(_output_directory +
1042                                                  TString::Format("Profile_plot_%s_vs_%s_%s_%d",
1043                                                                  x[ix].Data(),
1044                                                                  y[iy].Data(),
1045                                                                  _sublevel_names[isublevel - 1].Data(),
1046                                                                  _canvas_index) +
1047                                                  ExtensionFromPrintOption(_print_option),
1048                                              _print_option);
1049
1050           // draw into global canvas
1051           // printing will be performed after customisation (e.g. legend or title) just after the loops on ix and iy
1052           c_global_hist[isublevel]->cd(INDEX_IN_GLOBAL_CANVAS(ix, iy));
1053           c_global_hist[isublevel]
1054               ->GetPad(INDEX_IN_GLOBAL_CANVAS(ix, iy))
1055               ->SetFillStyle(4000);  //  make the pad transparent
1056           c_global_hist[isublevel]->GetPad(INDEX_IN_GLOBAL_CANVAS(ix, iy))->SetGrid(_grid_x, _grid_y);  // grid
1057           c_global_hist[isublevel]
1058               ->GetPad(INDEX_IN_GLOBAL_CANVAS(ix, iy))
1059               ->DrawFrame(_min[x[ix]],
1060                           CFG.count(min_branch.Data()) > 0 ? CFG.get<float>(min_branch.Data()) : _min[y[iy]],
1061                           _max[x[ix]],
1062                           CFG.count(max_branch.Data()) > 0 ? CFG.get<float>(max_branch.Data()) : _max[y[iy]],
1063                           TString(";") + LateXstyle(x[ix]) + " /" + _units[x[ix]] + TString(";") + LateXstyle(y[iy]) +
1064                               " /" + _units[y[iy]]);
1065
1066           histos[ix][iy][isublevel - 1]->Draw("same pe0");
1067           histos[ix][iy][NB_SUBLEVELS + isublevel - 1]->Draw("same pe0");
1068         }
1069
1070       }  // end of loop on y
1071     }    // end of loop on x
1072
1073     // CUSTOMISATION
1074     gStyle->SetOptTitle(0);  // otherwise, the title is repeated in every pad of the global canvases
1075                              // -> instead, we will write it in the upper part in a TPaveText or in a TLegend
1076     for (unsigned int ic = 0; ic <= NB_SUBLEVELS; ic++) {
1077       // setting legend to tracker canvases
1078       if (!_legend)
1079         break;
1080
1081       // setting legend to tracker canvases
1082       if (!_legend)
1083         break;
1084       TCanvas *c_temp_hist = (TCanvas *)c_global_hist[ic]->Clone(c_global_hist[ic]->GetTitle() + TString("_sub"));
1085       c_temp_hist->Draw();
1086       c_global_hist[ic] = new TCanvas(c_temp_hist->GetName() + TString("_final"),
1087                                       c_temp_hist->GetTitle(),
1088                                       c_temp_hist->GetWindowWidth(),
1089                                       c_temp_hist->GetWindowHeight());
1090       c_global_hist[ic]->Draw();
1091       TPad *p_up = new TPad(TString("legend_") + c_temp_hist->GetName(),
1092                             "",
1093                             0.,
1094                             0.9,
1095                             1.,
1096                             1.,  // relative position
1097                             -1,
1098                             0,
1099                             0),  // display options
1100           *p_down = new TPad(TString("main_") + c_temp_hist->GetName(), "", 0., 0., 1., 0.9, -1, 0, 0);
1101       // in the lower part, draw the plots
1102       p_down->Draw();
1103       p_down->cd();
1104       c_temp_hist->DrawClonePad();
1105       c_global_hist[ic]->cd();
1106       // in the upper part, pimp the canvas :p
1107       p_up->Draw();
1108       p_up->cd();
1109       if (ic == 0)  // tracker
1110       {
1111         TLegend *global_legend = MakeLegend(.05, .1, .7, .8, NB_SUBLEVELS);  //, "brNDC");
1112         global_legend->Draw();
1113         TPaveText *pt_geom = new TPaveText(.75, .1, .95, .8, "NB");
1114         pt_geom->SetFillColor(0);
1115         pt_geom->SetTextSize(0.25);
1116         pt_geom->AddText(TString("x: ") + _reference_name);
1117         pt_geom->AddText(TString("y: ") + _alignment_name + TString(" - ") + _reference_name);
1118         pt_geom->Draw();
1119       } else if (ic == 7)  // pixel
1120       {
1121         TLegend *global_legend = MakeLegend(.05, .1, .7, .8, 2);  //, "brNDC");
1122         global_legend->Draw();
1123         TPaveText *pt_geom = new TPaveText(.75, .1, .95, .8, "NB");
1124         pt_geom->SetFillColor(0);
1125         pt_geom->SetTextSize(0.25);
1126         pt_geom->AddText(TString("x: ") + _reference_name);
1127         pt_geom->AddText(TString("y: ") + _alignment_name + TString(" - ") + _reference_name);
1128         pt_geom->Draw();
1129       } else  // sublevels
1130       {
1131         TPaveText *pt = new TPaveText(.05, .1, .7, .8, "NB");
1132         pt->SetFillColor(0);
1133         pt->AddText(_sublevel_names[ic - 1]);
1134         pt->Draw();
1135         TPaveText *pt_geom = new TPaveText(.6, .1, .95, .8, "NB");
1136         pt_geom->SetFillColor(0);
1137         pt_geom->SetTextSize(0.3);
1138         pt_geom->AddText(TString("x: ") + _reference_name);
1139         pt_geom->AddText(TString("y: ") + _alignment_name + TString(" - ") + _reference_name);
1140         pt_geom->Draw();
1141       }
1142       // printing
1143       if (_print)
1144         c_global_hist[ic]->Print(
1145             _output_directory + c_global_hist[ic]->GetName() + ExtensionFromPrintOption(_print_option), _print_option);
1146     }
1147   }
1148
1149 #ifdef TALKATIVE
1150   std::cout << __FILE__ << ":" << __LINE__ << ":Info: End of MakePlots method" << std::endl;
1151 #endif
1152 }
1153
1154 // Make additional table for the mean/RMS values of differences
1155 void GeometryComparisonPlotter::MakeTables(
1156     std::vector<TString> x,  // axes to combine to plot
1157     std::vector<TString> y,  // only requires the differences (y values in the plots) and ranges
1158     pt::ptree CFG)           // property tree storing the min and max values under <var>_min and <var>_max
1159 {
1160   /// -1) check that only existing branches are called
1161   // (we use a macro to avoid copy/paste)
1162 #define CHECK_BRANCHES(branchname_vector)                                                       \
1163   for (unsigned int i = 0; i < branchname_vector.size(); i++) {                                 \
1164     if (branch_f.find(branchname_vector[i]) == branch_f.end()) {                                \
1165       std::cout << __FILE__ << ":" << __LINE__ << ":Error: The branch " << branchname_vector[i] \
1166                 << " is not recognised." << std::endl;                                          \
1167       return;                                                                                   \
1168     }                                                                                           \
1169   }
1170   CHECK_BRANCHES(x);
1171   CHECK_BRANCHES(y);
1172
1173   const unsigned int nentries = data->GetEntries();
1174
1175 #ifdef TALKATIVE
1176   std::cout << __FILE__ << ":" << __LINE__ << ":Info: ";
1177   INSIDE_VECTOR(x);
1178   std::cout << std::endl;
1179   std::cout << __FILE__ << ":" << __LINE__ << ":Info: ";
1180   INSIDE_VECTOR(y);
1181   std::cout << std::endl;
1182 #endif
1183
1184   /// 0) min and max values
1185   // the max and min of the graphs are computed from the tree if they have not been manually input yet
1186   // (we use a macro to avoid copy/paste)
1187 #define LIMITS(axes_vector)                                                          \
1188   for (unsigned int i = 0; i < axes_vector.size(); i++) {                            \
1189     if (_SF.find(axes_vector[i]) == _SF.end())                                       \
1190       _SF[axes_vector[i]] = 1.;                                                      \
1191     if (_min.find(axes_vector[i]) == _min.end())                                     \
1192       _min[axes_vector[i]] = _SF[axes_vector[i]] * data->GetMinimum(axes_vector[i]); \
1193     if (_max.find(axes_vector[i]) == _max.end())                                     \
1194       _max[axes_vector[i]] = _SF[axes_vector[i]] * data->GetMaximum(axes_vector[i]); \
1195   }
1196   LIMITS(x);
1197   LIMITS(y);
1198
1199 #ifdef TALKATIVE
1200   CHECK_MAP_CONTENT(_min, float);
1201   CHECK_MAP_CONTENT(_max, float);
1202   CHECK_MAP_CONTENT(_SF, float);
1203 #endif
1204
1205   /// 1) declare histograms
1206   // the idea is to produce tables of the differences and the absolute positions containing mean and RMS values
1207   // for the different subdetectors - 0..5=different sublevels.
1208   // Values for each endcap detector are to be split in +/-z, for the barrel detectors in +/- x (half barrels)
1209   // Since it is easier to handle in the loops, all subdetectors will be split in
1210   // 4 parts at first: (+/-x)X(+/-z)
1211   // This means that 2*2*6 histograms will be filled during the loop on the TTree
1212   // Pairs of histograms need to be combined afterwards again
1213   // Histograms 0-5 are at +x and +z, 6-11 at +x and -z, 12-17 at -x and +z, and 18-23 at -x and -z
1214   //
1215   // Two version of the table containing the differences are produced. Once using Gaussian fits (more stable
1216   // vs single outliers but perform poorly if the distributions are non-Gaussian) and once using
1217   // the mean and RMS of the histograms (more stable but outliers have a strong impact on the RMS).
1218   // For the absolute positions, only mean+RMS are used since the detector layout is not Gaussian
1219   // (structures due to layers/rings etc)
1220 #ifndef NB_SUBLEVELS
1221 #define NB_SUBLEVELS 6
1222 #endif
1223 #define NB_Z_SLICES 2
1224 #define NB_X_SLICES 2
1225
1226   TH1F *histosx[x.size()][NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES];
1227   float meanValuex[x.size()][NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES];
1228   float RMSx[x.size()][NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES];
1229
1230   TH1F *histos[y.size()][NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES];
1231   TF1 *gausFit[y.size()][NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES];
1232   float meanValue[y.size()][NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES];
1233   float meanValueGaussian[y.size()][NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES];
1234   float RMS[y.size()][NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES];
1235   float RMSGaussian[y.size()][NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES];
1236
1237   for (unsigned int iy = 0; iy < y.size(); iy++) {
1238     for (unsigned int ihist = 0; ihist < NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES; ihist++) {
1239       // Create and correctly name a histogram for each subdetector*Z_Slice*X_Slice
1240       histos[iy][ihist] = new TH1F("hist" + y[iy] + _sublevel_names[ihist % NB_SUBLEVELS] +
1241                                        TString(ihist % (NB_SUBLEVELS * NB_Z_SLICES) >= NB_SUBLEVELS ? "zn" : "zp") +
1242                                        TString(ihist >= NB_SUBLEVELS * NB_Z_SLICES ? "xn" : "xp"),
1243                                    "",
1244                                    1000,
1245                                    _min[y[iy]],
1246                                    _max[y[iy]] + 1.);
1247       histos[iy][ihist]->StatOverflows(kTRUE);
1248     }
1249   }
1250
1251   for (unsigned int ix = 0; ix < x.size(); ix++) {
1252     for (unsigned int ihist = 0; ihist < NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES; ihist++) {
1253       // Create and correctly name a histogram for each subdetector*Z_Slice*ModuleType
1254       histosx[ix][ihist] = new TH1F("histx" + x[ix] + _sublevel_names[ihist % NB_SUBLEVELS] +
1255                                         TString(ihist % (NB_SUBLEVELS * NB_Z_SLICES) >= NB_SUBLEVELS ? "zn" : "zp") +
1256                                         TString(ihist >= NB_SUBLEVELS * NB_Z_SLICES ? "xn" : "xp"),
1257                                     "",
1258                                     1000,
1259                                     _min[x[ix]],
1260                                     _max[x[ix]] + 1.);
1261       histosx[ix][ihist]->StatOverflows(kTRUE);
1262     }
1263   }
1264
1265 #ifdef DEBUG
1266   std::cout << __FILE__ << ":" << __LINE__ << ":Info: Creation of the TH1F[" << y.size() << "]["
1267             << NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES << "] ended." << std::endl;
1268 #endif
1269
1270   /// 2) loop on the TTree data
1271 #ifdef DEBUG
1272   std::cout << __FILE__ << ":" << __LINE__ << ":Info: Looping on the TTree" << std::endl;
1273 #endif
1274 #ifdef TALKATIVE
1275   unsigned int progress = 0;
1276   std::cout << __FILE__ << ":" << __LINE__ << ":Info: 0%" << std::endl;
1277 #endif
1278   for (unsigned int ientry = 0; ientry < nentries; ientry++) {
1279 #ifdef TALKATIVE
1280     if (10 * ientry / nentries != progress) {
1281       progress = 10 * ientry / nentries;
1282       std::cout << __FILE__ << ":" << __LINE__ << ":Info: " << 10 * progress << "%" << std::endl;
1283     }
1284 #endif
1285     // load current tree entry
1286     data->GetEntry(ientry);
1287
1288     // CUTS on entry
1289     if (branch_i["level"] != _levelCut)
1290       continue;
1291     if (!_1dModule && branch_i["detDim"] == 1)
1292       continue;
1293     if (!_2dModule && branch_i["detDim"] == 2)
1294       continue;
1295
1296     for (unsigned int iy = 0; iy < y.size(); iy++) {
1297       if (branch_i["sublevel"] < 1 || branch_i["sublevel"] > NB_SUBLEVELS)
1298         continue;
1299       if (_SF[y[iy]] * branch_f[y[iy]] > _max[y[iy]] || _SF[y[iy]] * branch_f[y[iy]] < _min[y[iy]]) {
1300         //#ifdef DEBUG
1301         //                    std::cout << "branch_f[y[iy]]=" << branch_f[y[iy]] << std::endl;
1302         //#endif
1303         continue;
1304       }
1305
1306       // FILLING HISTOGRAMS
1307
1308       // histogram for all modules
1309       const short int ihisto = (branch_i["sublevel"] - 1) + (branch_f["z"] >= 0 ? 0 : NB_SUBLEVELS) +
1310                                (branch_f["x"] >= 0 ? 0 : NB_SUBLEVELS * NB_Z_SLICES);
1311
1312       if (_module_plot_option == "all")
1313         histos[iy][ihisto]->Fill(_SF[y[iy]] * branch_f[y[iy]]);
1314
1315       // Only good modules
1316       else if (_module_plot_option == "good" && branch_i["badModuleQuality"] == 0)
1317         histos[iy][ihisto]->Fill(_SF[y[iy]] * branch_f[y[iy]]);
1318
1319       // Only good modules and those in the list
1320       else if (_module_plot_option == "list" && (branch_i["inModuleList"] == 1 || branch_i["badModuleQuality"] == 0))
1321         histos[iy][ihisto]->Fill(_SF[y[iy]] * branch_f[y[iy]]);
1322     }
1323
1324     for (unsigned int ix = 0; ix < x.size(); ix++) {
1325       if (branch_i["sublevel"] < 1 || branch_i["sublevel"] > NB_SUBLEVELS)
1326         continue;
1327       if (_SF[x[ix]] * branch_f[x[ix]] > _max[x[ix]] || _SF[x[ix]] * branch_f[x[ix]] < _min[x[ix]]) {
1328         //#ifdef DEBUG
1329         //                    std::cout << "branch_f[y[iy]]=" << branch_f[y[iy]] << std::endl;
1330         //#endif
1331         continue;
1332       }
1333
1334       // FILLING HISTOGRAMS
1335
1336       // histogram for all modules
1337       const short int ihistosx = (branch_i["sublevel"] - 1) + (branch_f["z"] >= 0 ? 0 : NB_SUBLEVELS) +
1338                                  (branch_f["x"] >= 0 ? 0 : NB_SUBLEVELS * NB_Z_SLICES);
1339
1340       if (_module_plot_option == "all")
1341         histosx[ix][ihistosx]->Fill(_SF[x[ix]] * branch_f[x[ix]]);
1342
1343       // Only good modules
1344       else if (_module_plot_option == "good" && branch_i["badModuleQuality"] == 0)
1345         histosx[ix][ihistosx]->Fill(_SF[x[ix]] * branch_f[x[ix]]);
1346
1347       // Only good modules and those in the list
1348       else if (_module_plot_option == "list" && (branch_i["inModuleList"] == 1 || branch_i["badModuleQuality"] == 0))
1349         histosx[ix][ihistosx]->Fill(_SF[x[ix]] * branch_f[x[ix]]);
1350     }
1351   }
1352 #ifdef TALKATIVE
1353   std::cout << __FILE__ << ":" << __LINE__ << ":Info: 100%\tLoop ended" << std::endl;
1354 #endif
1355
1356   //~ TString rangeLabel = "";
1357   // Calculate mean and standard deviation for each histogram
1358   for (unsigned int iy = 0; iy < y.size(); iy++) {
1359     for (unsigned int ihist = 0; ihist < NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES; ihist++) {
1360       // combine +/-z histograms for barrel detectors
1361       if (ihist % (NB_SUBLEVELS * NB_Z_SLICES) == 0 || ihist % (NB_SUBLEVELS * NB_Z_SLICES) == 2 ||
1362           ihist % (NB_SUBLEVELS * NB_Z_SLICES) == 4) {
1363         histos[iy][ihist]->Add(histos[iy][ihist + NB_SUBLEVELS]);
1364       }
1365       // combine +/-x histograms for endcap detectors (only used for half shells in barrel)
1366       if (ihist < NB_SUBLEVELS * NB_Z_SLICES &&
1367           (ihist % NB_SUBLEVELS == 1 || ihist % NB_SUBLEVELS == 3 || ihist % NB_SUBLEVELS == 5)) {
1368         histos[iy][ihist]->Add(histos[iy][ihist + NB_SUBLEVELS * NB_Z_SLICES]);
1369       }
1370       meanValue[iy][ihist] = histos[iy][ihist]->GetMean();
1371       RMS[iy][ihist] = histos[iy][ihist]->GetRMS();
1372
1373       histos[iy][ihist]->Fit("gaus");
1374       gausFit[iy][ihist] = histos[iy][ihist]->GetFunction("gaus");
1375       meanValueGaussian[iy][ihist] = gausFit[iy][ihist]->GetParameter(1);
1376       RMSGaussian[iy][ihist] = gausFit[iy][ihist]->GetParameter(2);
1377     }
1378   }
1379
1380   for (unsigned int ix = 0; ix < x.size(); ix++) {
1381     for (unsigned int ihist = 0; ihist < NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES; ihist++) {
1382       // combine +/-z histograms for barrel detectors
1383       if (ihist % (NB_SUBLEVELS * NB_Z_SLICES) == 0 || ihist % (NB_SUBLEVELS * NB_Z_SLICES) == 2 ||
1384           ihist % (NB_SUBLEVELS * NB_Z_SLICES) == 4) {
1385         histosx[ix][ihist]->Add(histosx[ix][ihist + NB_SUBLEVELS]);
1386       }
1387       // combine +/-x histograms for endcap detectors (only used for half shells in barrel)
1388       if (ihist < NB_SUBLEVELS * NB_Z_SLICES &&
1389           (ihist % NB_SUBLEVELS == 1 || ihist % NB_SUBLEVELS == 3 || ihist % NB_SUBLEVELS == 5)) {
1390         histosx[ix][ihist]->Add(histosx[ix][ihist + NB_SUBLEVELS * NB_Z_SLICES]);
1391       }
1392       meanValuex[ix][ihist] = histosx[ix][ihist]->GetMean();
1393       RMSx[ix][ihist] = histosx[ix][ihist]->GetRMS();
1394     }
1395   }
1396
1397   TString tableFileName, tableCaption, tableAlign, tableHeadline;
1398   TString PXBpLine, PXBmLine, PXFpLine, PXFmLine, TIBpLine, TIBmLine, TOBpLine, TOBmLine, TIDpLine, TIDmLine, TECpLine,
1399       TECmLine;
1400
1401   // table using mean and RMS, round to integers in µm etc.
1402   tableFileName = "table_differences.tex";
1403   if (_module_plot_option == "all")
1404     tableCaption = "Means and standard deviations of " + _alignment_name + " - " + _reference_name +
1405                    " for each subdetector, all modules used.";
1406   else if (_module_plot_option == "good")
1407     tableCaption = "Means and standard deviations of " + _alignment_name + " - " + _reference_name +
1408                    " for each subdetector, only good modules used.";
1409   else if (_module_plot_option == "list")
1410     tableCaption = "Means and standard deviations of " + _alignment_name + " - " + _reference_name +
1411                    " for each subdetector, good modules and those in given list used.";
1412
1413   WriteTable(y, NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES, meanValue, RMS, "0", tableCaption, tableFileName);
1414
1415   //~ // table using  Gaussian fit, round to integers in µm etc.
1416   tableFileName = "table_differences_Gaussian.tex";
1417   if (_module_plot_option == "all")
1418     tableCaption = "Means and standard deviations for Gaussian fit of " + _alignment_name + " - " + _reference_name +
1419                    " for each subdetector, all modules used.";
1420   else if (_module_plot_option == "good")
1421     tableCaption = "Means and standard deviations for Gaussian fit of " + _alignment_name + " - " + _reference_name +
1422                    " for each subdetector, only good modules used.";
1423   else if (_module_plot_option == "list")
1424     tableCaption = "Means and standard deviations for Gaussian fit of " + _alignment_name + " - " + _reference_name +
1425                    " for each subdetector, good modules and those in given list used.";
1426
1427   WriteTable(
1428       y, NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES, meanValueGaussian, RMSGaussian, "0", tableCaption, tableFileName);
1429
1430   // Table for the mean positions on the x-axis, round to 3 digits in cm etc.
1431   tableFileName = "table_meanPos.tex";
1432
1433   if (_module_plot_option == "all")
1434     tableCaption = "Mean positions and standard deviations in " + _reference_name +
1435                    " geometry for each subdetector, all modules used.";
1436   else if (_module_plot_option == "good")
1437     tableCaption = "Mean positions and standard deviations in " + _reference_name +
1438                    " geometry for each subdetector, only good modules used.";
1439   else if (_module_plot_option == "list")
1440     tableCaption = "Mean positions and standard deviations in " + _reference_name +
1441                    " geometry for each subdetector, good modules and those in given list used.";
1442
1443   WriteTable(x, NB_SUBLEVELS * NB_Z_SLICES * NB_X_SLICES, meanValuex, RMSx, "3", tableCaption, tableFileName);
1444
1445 #ifdef TALKATIVE
1446   std::cout << __FILE__ << ":" << __LINE__ << ":Info: End of MakeLegends method" << std::endl;
1447 #endif
1448 }
1449
1450 // OPTION METHODS
1451 void GeometryComparisonPlotter::SetPrint(const bool kPrint) { _print = kPrint; }
1452 void GeometryComparisonPlotter::SetLegend(const bool kLegend) { _legend = kLegend; }
1453 void GeometryComparisonPlotter::SetWrite(const bool kWrite) { _write = kWrite; }
1454 void GeometryComparisonPlotter::Set1dModule(const bool k1dModule) { _1dModule = k1dModule; }
1455 void GeometryComparisonPlotter::Set2dModule(const bool k2dModule) { _2dModule = k2dModule; }
1456 void GeometryComparisonPlotter::SetLevelCut(const int kLevelCut) { _levelCut = kLevelCut; }
1457 void GeometryComparisonPlotter::SetBatchMode(const bool kBatchMode) { _batchMode = kBatchMode; }
1458 void GeometryComparisonPlotter::SetGrid(const int kGridX, const int kGridY) {
1459   _grid_x = kGridX;
1460   _grid_y = kGridY;
1461 }
1462 void GeometryComparisonPlotter::SetBranchMax(const TString branchname, const float max) { _max[branchname] = max; }
1463 void GeometryComparisonPlotter::SetBranchMin(const TString branchname, const float min) { _min[branchname] = min; }
1464 void GeometryComparisonPlotter::SetBranchSF(const TString branchname, const float SF) { _SF[branchname] = SF; }
1465 void GeometryComparisonPlotter::SetBranchUnits(const TString branchname, const TString units) {
1466   _units[branchname] = units;
1467 }
1468 void GeometryComparisonPlotter::SetPrintOption(const Option_t *print_option) { _print_option = print_option; }
1469 void GeometryComparisonPlotter::SetCanvasSize(const int window_width, const int window_height) {
1470   _window_width = window_width;
1471   _window_height = window_height;
1472 }
1473 void GeometryComparisonPlotter::SetOutputFileName(const TString name) { _output_filename = name; }
1474 void GeometryComparisonPlotter::SetOutputDirectoryName(const TString name) {
1475   _output_directory = name + TString(name.EndsWith("/") ? "" : "/");
1476 }
1477
1478 // PRIVATE METHODS
1479 TString GeometryComparisonPlotter::LateXstyle(TString word) {
1480   word.ToLower();
1481   if (word.BeginsWith("d"))
1482     word.ReplaceAll("d", "#Delta");
1483   if (word == TString("rdphi"))
1484     word = "r#Delta#phi";  // TO DO: find something less ad hoc...
1485   else if (word.EndsWith("phi"))
1486     word.ReplaceAll("phi", "#phi");
1487   else if (word.EndsWith("alpha"))
1488     word.ReplaceAll("alpha", "#alpha");
1489   else if (word.EndsWith("beta"))
1490     word.ReplaceAll("beta", "#beta");
1491   else if (word.EndsWith("gamma"))
1492     word.ReplaceAll("gamma", "#gamma");
1493   else if (word.EndsWith("eta"))
1494     word.ReplaceAll("eta", "#eta");
1495   return word;
1496 }
1497
1498 TString GeometryComparisonPlotter::LateXstyleTable(TString word) {
1499   word.ToLower();
1500   if (word.BeginsWith("d"))
1501     word.ReplaceAll("d", "$\\Delta$");
1502   if (word == TString("rdphi"))
1503     word = "r$\\Delta\\phi$";  // TO DO: find something less ad hoc...
1504   else if (word.EndsWith("phi"))
1505     word.ReplaceAll("phi", "$\\phi$");
1506   else if (word.EndsWith("alpha"))
1507     word.ReplaceAll("alpha", "$\\alpha$");
1508   else if (word.EndsWith("beta"))
1509     word.ReplaceAll("beta", "$\\beta$");
1510   else if (word.EndsWith("gamma"))
1511     word.ReplaceAll("gamma", "#$\\gamma$");
1512   else if (word.EndsWith("eta"))
1513     word.ReplaceAll("eta", "$\\eta$");
1514   return word;
1515 }
1516
1517 TString GeometryComparisonPlotter::ExtensionFromPrintOption(TString print_option) {
1518   if (print_option.Contains("pdf"))
1519     return TString(".pdf");
1520   else if (print_option.Contains("eps"))
1521     return TString(".eps");
1522   else if (print_option.Contains("ps"))
1523     return TString(".ps");
1524   else if (print_option.Contains("svg"))
1525     return TString(".svg");
1526   else if (print_option.Contains("tex"))
1527     return TString(".tex");
1528   else if (print_option.Contains("gif"))
1529     return TString(".gif");
1530   else if (print_option.Contains("xpm"))
1531     return TString(".xpm");
1532   else if (print_option.Contains("png"))
1533     return TString(".png");
1534   else if (print_option.Contains("jpg"))
1535     return TString(".jpg");
1536   else if (print_option.Contains("tiff"))
1537     return TString(".tiff");
1538   else if (print_option.Contains("cxx"))
1539     return TString(".cxx");
1540   else if (print_option.Contains("xml"))
1541     return TString(".xml");
1542   else if (print_option.Contains("root"))
1543     return TString(".root");
1544   else {
1545     std::cout << __FILE__ << ":" << __LINE__ << ":Warning: unknown format. Returning .pdf, but possibly wrong..."
1546               << std::endl;
1547     return TString(".pdf");
1548   }
1549 }
1550
1551 TLegend *GeometryComparisonPlotter::MakeLegend(
1552     double x1, double y1, double x2, double y2, int nPlottedSublevels, const TString title) {
1553   TLegend *legend = new TLegend(x1, y1, x2, y2, title.Data(), "NBNDC");
1554   legend->SetNColumns(nPlottedSublevels);
1555   legend->SetFillColor(0);
1556   legend->SetLineColor(0);  // redundant with option
1557   legend->SetLineWidth(0);  // redundant with option
1558   for (int isublevel = 0; isublevel < nPlottedSublevels;
1559        isublevel++)  // nPlottedSublevels is either NB_SUBLEVELS for the tracker or 2 for the pixel
1560   {
1561     TGraph *g = new TGraph(0);
1562     g->SetMarkerColor(COLOR_CODE(isublevel));
1563     g->SetFillColor(COLOR_CODE(isublevel));
1564     g->SetMarkerStyle(kFullSquare);
1565     g->SetMarkerSize(10);
1566     legend->AddEntry(g, _sublevel_names[isublevel], "p");
1567   }
1568   return legend;
1569 }
1570
1571 void GeometryComparisonPlotter::WriteTable(std::vector<TString> x,
1572                                            unsigned int nLevelsTimesSlices,
1573                                            float meanValue[10][24],
1574                                            float RMS[10][24],
1575                                            const TString nDigits,
1576                                            const TString tableCaption,
1577                                            const TString tableFileName) {
1578   std::ofstream output(_output_directory + tableFileName);
1579
1580   TString tableAlign, tableHeadline;
1581   TString PXBpLine, PXBmLine, PXFpLine, PXFmLine, TIBpLine, TIBmLine, TOBpLine, TOBmLine, TIDpLine, TIDmLine, TECpLine,
1582       TECmLine;
1583   char meanChar[x.size()][nLevelsTimesSlices][10];
1584   char RMSChar[x.size()][nLevelsTimesSlices][10];
1585
1586   tableAlign = "l";
1587   tableHeadline = "";
1588   PXBpLine = "PXB x$+$";
1589   PXBmLine = "PXB x$-$";
1590   PXFpLine = "PXF z$+$";
1591   PXFmLine = "PXF z$-$";
1592   TIBpLine = "TIB x$+$";
1593   TIBmLine = "TIB x$-$";
1594   TIDpLine = "TID z$+$";
1595   TIDmLine = "TID z$-$";
1596   TOBpLine = "TOB x$+$";
1597   TOBmLine = "TOB x$-$";
1598   TECpLine = "TEC z$+$";
1599   TECmLine = "TEC z$-$";
1600
1601   for (unsigned int ix = 0; ix < x.size(); ix++) {
1602     for (unsigned int isubDet = 0; isubDet < nLevelsTimesSlices; isubDet++) {
1603       sprintf(meanChar[ix][isubDet], "%." + nDigits + "f", meanValue[ix][isubDet]);
1604       sprintf(RMSChar[ix][isubDet], "%." + nDigits + "f", RMS[ix][isubDet]);
1605     }
1606     tableAlign += "|c";
1607     tableHeadline += " & " + LateXstyleTable(x[ix]) + " / " + _units[x[ix]].ReplaceAll("#mum", "$\\mu$m");
1608
1609     PXBpLine += " & $";
1610     PXBpLine += meanChar[ix][0];
1611     PXBpLine += "\\pm";
1612     PXBpLine += RMSChar[ix][0];
1613     PXBpLine += " $";
1614     PXBmLine += " & $";
1615     PXBmLine += meanChar[ix][12];
1616     PXBmLine += "\\pm";
1617     PXBmLine += RMSChar[ix][12];
1618     PXBmLine += " $";
1619     PXFpLine += " & $";
1620     PXFpLine += meanChar[ix][1];
1621     PXFpLine += "\\pm";
1622     PXFpLine += RMSChar[ix][1];
1623     PXFpLine += " $";
1624     PXFmLine += " & $";
1625     PXFmLine += meanChar[ix][7];
1626     PXFmLine += "\\pm";
1627     PXFmLine += RMSChar[ix][7];
1628     PXFmLine += " $";
1629     TIBpLine += " & $";
1630     TIBpLine += meanChar[ix][2];
1631     TIBpLine += "\\pm";
1632     TIBpLine += RMSChar[ix][2];
1633     TIBpLine += " $";
1634     TIBmLine += " & $";
1635     TIBmLine += meanChar[ix][14];
1636     TIBmLine += "\\pm";
1637     TIBmLine += RMSChar[ix][14];
1638     TIBmLine += " $";
1639     TIDpLine += " & $";
1640     TIDpLine += meanChar[ix][3];
1641     TIDpLine += "\\pm";
1642     TIDpLine += RMSChar[ix][3];
1643     TIDpLine += " $";
1644     TIDmLine += " & $";
1645     TIDmLine += meanChar[ix][9];
1646     TIDmLine += "\\pm";
1647     TIDmLine += RMSChar[ix][9];
1648     TIDmLine += " $";
1649     TOBpLine += " & $";
1650     TOBpLine += meanChar[ix][4];
1651     TOBpLine += "\\pm";
1652     TOBpLine += RMSChar[ix][4];
1653     TOBpLine += " $";
1654     TOBmLine += " & $";
1655     TOBmLine += meanChar[ix][16];
1656     TOBmLine += "\\pm";
1657     TOBmLine += RMSChar[ix][16];
1658     TOBmLine += " $";
1659     TECpLine += " & $";
1660     TECpLine += meanChar[ix][5];
1661     TECpLine += "\\pm";
1662     TECpLine += RMSChar[ix][5];
1663     TECpLine += " $";
1664     TECmLine += " & $";
1665     TECmLine += meanChar[ix][11];
1666     TECmLine += "\\pm";
1667     TECmLine += RMSChar[ix][11];
1668     TECmLine += " $";
1669   }
1670
1671   // Write the table to the tex file
1672   output << "\\begin{table}" << std::endl;
1673   output << "\\caption{" << tableCaption << "}" << std::endl;
1674   output << "\\begin{tabular}{" << tableAlign << "}" << std::endl;
1675   output << "\\hline" << std::endl;
1676   output << tableHeadline << " \\\\" << std::endl;
1677   output << "\\hline" << std::endl;
1678   output << PXBpLine << " \\\\" << std::endl;
1679   output << PXBmLine << " \\\\" << std::endl;
1680   output << PXFpLine << " \\\\" << std::endl;
1681   output << PXFmLine << " \\\\" << std::endl;
1682   output << TIBpLine << " \\\\" << std::endl;
1683   output << TIBmLine << " \\\\" << std::endl;
1684   output << TIDpLine << " \\\\" << std::endl;
1685   output << TIDmLine << " \\\\" << std::endl;
1686   output << TOBpLine << " \\\\" << std::endl;
1687   output << TOBmLine << " \\\\" << std::endl;
1688   output << TECpLine << " \\\\" << std::endl;
1689   output << TECmLine << " \\\\" << std::endl;
1690   output << "\\hline" << std::endl;
1691   output << "\\end{tabular}" << std::endl;
1692   output << "\\end{table}" << std::endl;
1693 }