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