OfflineValidation/bin/DMRtrends.cc

0001 #include <cstdlib>
0002 #include <string>
0003 #include <tuple>
0004 #include <iostream>
0005 #include <numeric>
0006 #include <functional>
0007 #include <unistd.h>
0008
0009 #include "TFile.h"
0010 #include "TGraph.h"
0011 #include "TH1.h"
0012
0013 #include "exceptions.h"
0014 #include "toolbox.h"
0015 #include "Options.h"
0016
0017 #include "FWCore/ParameterSet/interface/FileInPath.h"
0018 #include "boost/filesystem.hpp"
0019 #include "boost/algorithm/string.hpp"
0020 #include "boost/property_tree/ptree.hpp"
0021 #include "boost/property_tree/json_parser.hpp"
0022 #include "boost/optional.hpp"
0023
0024 #include "TString.h"
0025 #include "TColor.h"
0026
0027 #include "Alignment/OfflineValidation/interface/PrepareDMRTrends.h"
0028 #include "Alignment/OfflineValidation/interface/Trend.h"
0029
0030 using namespace std;
0031 using namespace AllInOneConfig;
0032 namespace fs = boost::filesystem;
0033 namespace bc = boost::container;
0034
0035 static const char *bold = "\e[1m", *normal = "\e[0m";
0036 static const float defaultConvertScale = 1000.;
0037 static const int startRun2016 = 272930;
0038 static const int endRun2018 = 325175;
0039
0040 namespace pt = boost::property_tree;
0041
0042 int trends(int argc, char *argv[]) {
0043   // parse the command line
0044
0045   Options options;
0046   options.helper(argc, argv);
0047   options.parser(argc, argv);
0048
0049   //Read in AllInOne json config
0050   pt::ptree main_tree;
0051   pt::read_json(options.config, main_tree);
0052
0053   pt::ptree alignments = main_tree.get_child("alignments");
0054   pt::ptree validation = main_tree.get_child("validation");
0055   pt::ptree style = main_tree.get_child("style");
0056
0057   //Read all configure variables and set default for missing keys
0058   string outputdir = main_tree.get<string>("output");
0059   bool FORCE = validation.count("FORCE") ? validation.get<bool>("FORCE") : false;
0060   string year = validation.count("year") ? validation.get<string>("year") : "Run2";
0061   TString lumiInputFile = style.get_child("trends").count("lumiInputFile")
0062                               ? style.get_child("trends").get<string>("lumiInputFile")
0063                               : "Alignment/OfflineValidation/data/lumiPerRun_Run2.txt";
0064   fs::path lumiFile = lumiInputFile.Data();
0065   edm::FileInPath fip = edm::FileInPath(lumiFile.string());
0066   fs::path pathToLumiFile = "";
0067   if (!fs::exists(lumiFile)) {
0068     pathToLumiFile = fip.fullPath();
0069   } else {
0070     pathToLumiFile = lumiFile;
0071   }
0072   if (!fs::exists(pathToLumiFile)) {
0073     cout << "ERROR: lumi-per-run file (" << lumiFile.string().data() << ") not found!" << endl
0074          << "Please check!" << endl;
0075     exit(EXIT_FAILURE);
0076   } else {
0077     cout << "Found lumi-per-run file: " << pathToLumiFile.string().data() << endl;
0078   }
0079   if (!lumiInputFile.Contains(year)) {
0080     cout << "ERROR: lumi-per-run file must contain (" << year.data() << ")!" << endl << "Please check!" << endl;
0081     exit(EXIT_FAILURE);
0082   }
0083
0084   string lumiAxisType = "recorded";
0085   if (lumiInputFile.Contains("delivered"))
0086     lumiAxisType = "delivered";
0087
0088   std::cout << Form("NOTE: using %s luminosity!", lumiAxisType.data()) << std::endl;
0089
0090   vector<int> IOVlist;
0091   vector<string> inputFiles;
0092   for (auto const &childTree : validation.get_child("IOV")) {
0093     int iov = childTree.second.get_value<int>();
0094     IOVlist.push_back(iov);
0095     TString mergeFile = validation.get<string>("mergeFile");
0096     string input = Form("%s/OfflineValidationSummary.root", mergeFile.ReplaceAll("{}", to_string(iov)).Data());
0097     inputFiles.push_back(input);
0098   }
0099
0100   string labels_to_add = "";
0101   if (validation.count("labels")) {
0102     for (auto const &label : validation.get_child("labels")) {
0103       labels_to_add += "_";
0104       labels_to_add += label.second.get_value<string>();
0105     }
0106   }
0107
0108   fs::path pname = Form("%s/DMRtrends%s.root", outputdir.data(), labels_to_add.data());
0109
0110   vector<TString> structures{"BPIX", "BPIX_y", "FPIX", "FPIX_y", "TIB", "TID", "TOB", "TEC"};
0111
0112   map<TString, int> nlayers{{"BPIX", 4}, {"FPIX", 3}, {"TIB", 4}, {"TID", 3}, {"TOB", 6}, {"TEC", 9}};
0113   if (year == "2016")
0114     nlayers = {{"BPIX", 3}, {"FPIX", 2}, {"TIB", 4}, {"TID", 3}, {"TOB", 6}, {"TEC", 9}};
0115
0116   PrepareDMRTrends prepareTrends(pname.c_str(), alignments);
0117   if (validation.count("Variables")) {
0118     for (auto const &Variable : validation.get_child("Variables")) {
0119       prepareTrends.compileDMRTrends(
0120           IOVlist, Variable.second.get_value<string>(), inputFiles, structures, nlayers, FORCE);
0121     }
0122   } else
0123     prepareTrends.compileDMRTrends(IOVlist, "median", inputFiles, structures, nlayers, FORCE);
0124
0125   assert(fs::exists(pname));
0126
0127   float convertUnit = style.get_child("trends").count("convertUnit")
0128                           ? style.get_child("trends").get<float>("convertUnit")
0129                           : defaultConvertScale;
0130   int firstRun = validation.count("firstRun") ? validation.get<int>("firstRun") : startRun2016;
0131   int lastRun = validation.count("lastRun") ? validation.get<int>("lastRun") : endRun2018;
0132
0133   const Run2Lumi GetLumi(pathToLumiFile.string().data(), firstRun, lastRun, convertUnit);
0134
0135   auto f = TFile::Open(pname.c_str());
0136
0137   for (auto const &Variable : validation.get_child("Variables")) {
0138     vector<tuple<TString, TString, float, float>> DMRs{{"mu", "#mu [#mum]", -6, 6},
0139                                                        {"sigma", "#sigma_{#mu} [#mum]", -15, 15},
0140                                                        {"muplus", "#mu outward [#mum]", -6, 6},
0141                                                        {"sigmaplus", "#sigma_{#mu outward} [#mum]", -15, 15},
0142                                                        {"muminus", "#mu inward [#mum]", -6, 6},
0143                                                        {"sigmaminus", "#sigma_{#mu inward} [#mum]", -15, 15},
0144                                                        {"deltamu", "#Delta#mu [#mum]", -15, 15},
0145                                                        {"sigmadeltamu", "#sigma_{#Delta#mu} [#mum]", -15, 15},
0146                                                        {"musigma", "#mu [#mum]", -6, 6},
0147                                                        {"muplussigmaplus", "#mu outward [#mum]", -15, 15},
0148                                                        {"muminussigmaminus", "#mu inward [#mum]", -15, 15},
0149                                                        {"deltamusigmadeltamu", "#Delta#mu [#mum]", -15, 15}};
0150
0151     if (Variable.second.get_value<string>() == "DrmsNR") {
0152       DMRs = {{"mu", "RMS(x'_{pred}-x'_{hit} /#sigma)", -1.2, 1.2},
0153               {"sigma", "#sigma_{RMS(x'_{pred}-x'_{hit} /#sigma)}", -6, 6},
0154               {"muplus", "RMS(x'_{pred}-x'_{hit} /#sigma) outward", -1.2, 1.2},
0155               {"sigmaplus", "#sigma_{#mu outward}", -6, 6},
0156               {"muminus", "RMS(x'_{pred}-x'_{hit} /#sigma) inward", -1.2, 1.2},
0157               {"sigmaminus", "#sigma_{RMS(x'_{pred}-x'_{hit} /#sigma) inward}", -6, 6},
0158               {"deltamu", "#DeltaRMS(x'_{pred}-x'_{hit} /#sigma)", -0.15, 0.15},
0159               {"sigmadeltamu", "#sigma_{#DeltaRMS(x'_{pred}-x'_{hit} /#sigma)}", -6, 6},
0160               {"musigma", "RMS(x'_{pred}-x'_{hit} /#sigma)", -1.2, 1.2},
0161               {"muplussigmaplus", "RMS(x'_{pred}-x'_{hit} /#sigma) outward", -1.2, 1.2},
0162               {"muminussigmaminus", "RMS(x'_{pred}-x'_{hit} /#sigma) inward", -1.2, 1.2},
0163               {"deltamusigmadeltamu", "#DeltaRMS(x'_{pred}-x'_{hit} /#sigma)", -0.15, 0.15}};
0164     }
0165
0166     for (const auto &structure : structures) {
0167       TString structname = structure;
0168       structname.ReplaceAll("_y", "");
0169       size_t layersnumber = nlayers.at(structname);
0170       for (Size_t layer = 0; layer <= layersnumber; layer++) {
0171         TString structtitle = "";
0172         if (structure.Contains("PIX") && !(structure.Contains("_y")))
0173           structtitle = structure + " (x)";
0174         else if (structure.Contains("_y")) {
0175           TString substring(structure(0, 4));
0176           structtitle = substring + " (y)";
0177         } else
0178           structtitle = structure;
0179         if (layer != 0) {
0180           if (structure == "TID" || structure == "TEC" || structure == "FPIX" || structure == "FPIX_y")
0181             structtitle += "  disc ";
0182           else
0183             structtitle += "  layer ";
0184           structtitle += layer;
0185         }
0186
0187         TString structandlayer = structure;
0188         if (layer != 0) {
0189           if (structure == "TID" || structure == "TEC")
0190             structandlayer += "_disc";
0191           else
0192             structandlayer += "_layer";
0193           structandlayer += layer;
0194         }
0195
0196         for (auto &DMR : DMRs) {
0197           auto name = get<0>(DMR), ytitle = get<1>(DMR);
0198
0199           if (name.Contains("plus") || name.Contains("minus") || name.Contains("delta")) {
0200             if (structname == "TEC" || structname == "TID")
0201               continue;  //Lorentz drift cannot appear in TEC and TID. These structures are skipped when looking at outward and inward pointing modules.
0202           }
0203
0204           cout << bold << name << normal << endl;
0205
0206           float ymin = get<2>(DMR), ymax = get<3>(DMR);
0207           Trend trend(Form("%s_%s_%s", Variable.second.get_value<string>().data(), structandlayer.Data(), name.Data()),
0208                       outputdir.data(),
0209                       ytitle,
0210                       ytitle,
0211                       ymin,
0212                       ymax,
0213                       style,
0214                       GetLumi,
0215                       lumiAxisType.data());
0216           trend.lgd.SetHeader(structtitle);
0217
0218           for (auto const &alignment : alignments) {
0219             bool fullRange = true;
0220             if (style.get_child("trends").count("earlyStops")) {
0221               for (auto const &earlyStop : style.get_child("trends.earlyStops")) {
0222                 if (earlyStop.second.get_value<string>() == alignment.first)
0223                   fullRange = false;
0224               }
0225             }
0226
0227             TString gtitle = alignment.second.get<string>("title");
0228             TString gname = Form("%s_%s_%s_%s",
0229                                  Variable.second.get_value<string>().data(),
0230                                  gtitle.Data(),
0231                                  structandlayer.Data(),
0232                                  name.Data());
0233             gname.ReplaceAll(" ", "_");
0234             auto g = Get<TGraphErrors>(gname);
0235             assert(g != nullptr);
0236             g->SetTitle(gtitle);  // for the legend
0237             g->SetMarkerSize(0.6);
0238             int color = alignment.second.get<int>("color");
0239             int style = floor(alignment.second.get<double>("style") / 100.);
0240             g->SetFillColorAlpha(color, 0.2);
0241             g->SetMarkerColor(color);
0242             g->SetMarkerStyle(style);
0243             g->SetLineColor(kWhite);
0244             trend(g, "P2", "pf", fullRange);
0245           }
0246         }
0247       }
0248     }
0249   }
0250
0251   f->Close();
0252   cout << bold << "Done" << normal << endl;
0253
0254   return EXIT_SUCCESS;
0255 }
0256
0257 #ifndef DOXYGEN_SHOULD_SKIP_THIS
0258 int main(int argc, char *argv[]) { return exceptions<trends>(argc, argv); }
0259 #endif