Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:45:39

0001 #include <iostream>
0002 #include <string>
0003 #include <vector>
0004 #include <algorithm>
0005 #include <map>
0006 #include <iomanip>
0007 #include <experimental/filesystem>
0008 #include "TPad.h"
0009 #include "TCanvas.h"
0010 #include "TGraph.h"
0011 #include "TGraphErrors.h"
0012 #include "TMultiGraph.h"
0013 #include "TH1.h"
0014 #include "THStack.h"
0015 #include "TROOT.h"
0016 #include "TFile.h"
0017 #include "TColor.h"
0018 #include "TLegend.h"
0019 #include "TLegendEntry.h"
0020 #include "TMath.h"
0021 #include "TRegexp.h"
0022 #include "TPaveLabel.h"
0023 #include "TPaveText.h"
0024 #include "TStyle.h"
0025 #include "TLine.h"
0026 #include "Alignment/OfflineValidation/plugins/ColorParser.C"
0027 
0028 using namespace std;
0029 namespace fs = std::experimental::filesystem;
0030 
0031 /*!
0032  * \def Dummy value in case a DMR would fail for instance
0033  */
0034 #define DUMMY -999.
0035 /*!
0036  * \def Scale factor value to have luminosity expressed in fb^-1
0037  */
0038 #define lumiFactor 1000.
0039 /*!
0040  * \def Scale factor value to have mean and sigmas expressed in micrometers.
0041  */
0042 #define DMRFactor 10000.
0043 
0044 /*! \struct Point
0045  *  \brief Structure Point
0046  *         Contains parameters of Gaussian fits to DMRs
0047  *  
0048  * @param run:             run number (IOV boundary)
0049  * @param scale:           scale for the measured quantity: cm->μm for DMRs, 1 for normalized residuals
0050  * @param mu:              mu/mean from Gaussian fit to DMR/DrmsNR
0051  * @param sigma:           sigma/standard deviation from Gaussian fit to DMR/DrmsNR
0052  * @param muplus:          mu/mean for the inward pointing modules
0053  * @param muminus:         mu/mean for outward pointing modules
0054  * @param sigmaplus:       sigma/standard for inward pointing modules 
0055  * @param sigmaminus: //!< sigma/standard for outward pointing modules
0056  */
0057 struct Point {
0058   float run, scale, mu, sigma, muplus, muminus, sigmaplus, sigmaminus;
0059 
0060   /*! \fn Point
0061      *  \brief Constructor of structure Point, initialising all members one by one
0062      */
0063   Point(float Run = DUMMY,
0064         float ScaleFactor = DMRFactor,
0065         float y1 = DUMMY,
0066         float y2 = DUMMY,
0067         float y3 = DUMMY,
0068         float y4 = DUMMY,
0069         float y5 = DUMMY,
0070         float y6 = DUMMY)
0071       : run(Run), scale(ScaleFactor), mu(y1), sigma(y2), muplus(y3), muminus(y5), sigmaplus(y4), sigmaminus(y6) {}
0072 
0073   /*! \fn Point
0074      *  \brief Constructor of structure Point, initialising all members from DMRs directly (with split)
0075      */
0076   Point(float Run, float ScaleFactor, TH1 *histo, TH1 *histoplus, TH1 *histominus)
0077       : Point(Run,
0078               ScaleFactor,
0079               histo->GetMean(),
0080               histo->GetMeanError(),
0081               histoplus->GetMean(),
0082               histoplus->GetMeanError(),
0083               histominus->GetMean(),
0084               histominus->GetMeanError()) {}
0085 
0086   /*! \fn Point
0087      *  \brief Constructor of structure Point, initialising all members from DMRs directly (without split)
0088      */
0089   Point(float Run, float ScaleFactor, TH1 *histo) : Point(Run, ScaleFactor, histo->GetMean(), histo->GetMeanError()) {}
0090 
0091   Point &operator=(const Point &p) {
0092     run = p.run;
0093     mu = p.mu;
0094     muplus = p.muplus;
0095     muminus = p.muminus;
0096     sigma = p.sigma;
0097     sigmaplus = p.sigmaplus;
0098     sigmaminus = p.sigmaminus;
0099     return *this;
0100   }
0101 
0102   float GetRun() const { return run; }
0103   float GetMu() const { return scale * mu; }
0104   float GetMuPlus() const { return scale * muplus; }
0105   float GetMuMinus() const { return scale * muminus; }
0106   float GetSigma() const { return scale * sigma; }
0107   float GetSigmaPlus() const { return scale * sigmaplus; }
0108   float GetSigmaMinus() const { return scale * sigmaminus; }
0109   float GetDeltaMu() const {
0110     if (muplus == DUMMY && muminus == DUMMY)
0111       return DUMMY;
0112     else
0113       return scale * (muplus - muminus);
0114   }
0115   float GetSigmaDeltaMu() const {
0116     if (sigmaplus == DUMMY && sigmaminus == DUMMY)
0117       return DUMMY;
0118     else
0119       return scale * hypot(sigmaplus, sigmaminus);
0120   };
0121 };
0122 
0123 ///**************************
0124 ///*  Function declaration  *
0125 ///**************************
0126 
0127 TString getName(TString structure, int layer, TString geometry);
0128 TH1F *ConvertToHist(TGraphErrors *g);
0129 const map<TString, int> numberOfLayers(TString Year = "2018");
0130 vector<int> runlistfromlumifile(TString Year = "2018");
0131 bool checkrunlist(vector<int> runs, vector<int> IOVlist = {}, TString Year = "2018");
0132 TString lumifileperyear(TString Year = "2018", string RunOrIOV = "IOV");
0133 void scalebylumi(TGraphErrors *g, vector<pair<int, double>> lumiIOVpairs);
0134 vector<pair<int, double>> lumiperIOV(vector<int> IOVlist, TString Year = "2018");
0135 double getintegratedlumiuptorun(int run, TString Year = "2018", double min = 0.);
0136 void PixelUpdateLines(TCanvas *c,
0137                       TString Year = "2018",
0138                       bool showlumi = false,
0139                       vector<int> pixelupdateruns = {314881, 316758, 317527, 318228, 320377});
0140 void PlotDMRTrends(
0141     vector<int> IOVlist,
0142     TString Variable = "median",
0143     vector<string> labels = {"MB"},
0144     TString Year = "2018",
0145     string myValidation = "/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/results/acardini/DMRs/",
0146     vector<string> geometries = {"GT", "SG", "MP pix LBL", "PIX HLS+ML STR fix"},
0147     vector<Color_t> colours = {kBlue, kRed, kGreen, kCyan},
0148     TString outputdir =
0149         "/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/alignmentObjects/acardini/DMRsTrends/",
0150     bool pixelupdate = false,
0151     vector<int> pixelupdateruns = {314881, 316758, 317527, 318228, 320377},
0152     bool showlumi = false);
0153 void compileDMRTrends(
0154     vector<int> IOVlist,
0155     TString Variable = "median",
0156     vector<string> labels = {"MB"},
0157     TString Year = "2018",
0158     string myValidation = "/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/results/acardini/DMRs/",
0159     vector<string> geometries = {"GT", "SG", "MP pix LBL", "PIX HLS+ML STR fix"},
0160     bool showlumi = false,
0161     bool FORCE = false);
0162 void DMRtrends(
0163     vector<int> IOVlist,
0164     vector<string> Variables = {"median", "DrmsNR"},
0165     vector<string> labels = {"MB"},
0166     TString Year = "2018",
0167     string myValidation = "/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/results/acardini/DMRs/",
0168     vector<string> geometries = {"GT", "SG", "MP pix LBL", "PIX HLS+ML STR fix"},
0169     vector<Color_t> colours = {kBlue, kRed, kGreen, kCyan},
0170     TString outputdir =
0171         "/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/alignmentObjects/acardini/DMRsTrends/",
0172     bool pixelupdate = false,
0173     vector<int> pixelupdateruns = {314881, 316758, 317527, 318228, 320377},
0174     bool showlumi = false,
0175     bool FORCE = false);
0176 
0177 /*! \class Geometry
0178  *  \brief Class Geometry
0179  *         Contains vector for fit parameters (mean, sigma, etc.) obtained from multiple IOVs
0180  *         See Structure Point for description of the parameters.
0181  */
0182 
0183 class Geometry {
0184 public:
0185   vector<Point> points;
0186 
0187 private:
0188   //template<typename T> vector<T> GetQuantity (T (Point::*getter)() const) const {
0189   vector<float> GetQuantity(float (Point::*getter)() const) const {
0190     vector<float> v;
0191     for (Point point : points) {
0192       float value = (point.*getter)();
0193       v.push_back(value);
0194     }
0195     return v;
0196   }
0197 
0198 public:
0199   TString title;
0200   Geometry() : title("") {}
0201   Geometry(TString Title) : title(Title) {}
0202   Geometry &operator=(const Geometry &geom) {
0203     title = geom.title;
0204     points = geom.points;
0205     return *this;
0206   }
0207   void SetTitle(TString Title) { title = Title; }
0208   TString GetTitle() { return title; }
0209   vector<float> Run() const { return GetQuantity(&Point::GetRun); }
0210   vector<float> Mu() const { return GetQuantity(&Point::GetMu); }
0211   vector<float> MuPlus() const { return GetQuantity(&Point::GetMuPlus); }
0212   vector<float> MuMinus() const { return GetQuantity(&Point::GetMuMinus); }
0213   vector<float> Sigma() const { return GetQuantity(&Point::GetSigma); }
0214   vector<float> SigmaPlus() const { return GetQuantity(&Point::GetSigmaPlus); }
0215   vector<float> SigmaMinus() const { return GetQuantity(&Point::GetSigmaMinus); }
0216   vector<float> DeltaMu() const { return GetQuantity(&Point::GetDeltaMu); }
0217   vector<float> SigmaDeltaMu() const { return GetQuantity(&Point::GetSigmaDeltaMu); }
0218   //vector<float> Graph (string variable) const {
0219   // };
0220 };
0221 
0222 /// DEPRECATED
0223 //struct Layer {
0224 //    map<string,Geometry> geometries;
0225 //};
0226 //
0227 //struct HLS {
0228 //    vector<Layer> layers;
0229 //    map<string,Geometry> geometries;
0230 //};
0231 
0232 /*! \fn getName
0233  *  \brief Function used to get a string containing information on the high level structure, the layer/disc and the geometry.
0234  */
0235 
0236 TString getName(TString structure, int layer, TString geometry) {
0237   geometry.ReplaceAll(" ", "_");
0238   TString name = geometry + "_" + structure;
0239   if (layer != 0) {
0240     if (structure == "TID" || structure == "TEC")
0241       name += "_disc";
0242     else
0243       name += "_layer";
0244     name += layer;
0245   }
0246 
0247   return name;
0248 };
0249 
0250 /*! \fn numberOfLayers
0251  *  \brief Function used to retrieve a map containing the number of layers per subdetector
0252  */
0253 
0254 const map<TString, int> numberOfLayers(TString Year) {
0255   if (Year == "2016")
0256     return {{"BPIX", 3}, {"FPIX", 2}, {"TIB", 4}, {"TID", 3}, {"TOB", 6}, {"TEC", 9}};
0257   else
0258     return {{"BPIX", 4}, {"FPIX", 3}, {"TIB", 4}, {"TID", 3}, {"TOB", 6}, {"TEC", 9}};
0259 }
0260 
0261 /// TO DO: once the information on the luminosity is passed through the root files this method needs to be changed
0262 /*! \fn lumifileperyear
0263  *  \brief Function to retrieve the file with luminosity per run/IOV
0264  *         The use of a lumi-per-IOV file is deprecated, but can still be useful for debugging
0265  */
0266 
0267 TString lumifileperyear(TString Year, string RunOrIOV) {
0268   TString LumiFile = std::getenv("CMSSW_BASE");
0269   LumiFile += "/src/Alignment/OfflineValidation/data/lumiper";
0270   if (RunOrIOV != "run" && RunOrIOV != "IOV") {
0271     cout << "ERROR: Please specify \"run\" or \"IOV\" to retrieve the luminosity run by run or for each IOV" << endl;
0272     exit(EXIT_FAILURE);
0273   }
0274   LumiFile += RunOrIOV;
0275   if (Year != "2016" && Year != "2017" && Year != "2018") {
0276     cout << "ERROR: Only 2016, 2017 and 2018 lumi-per-run files are available, please check!" << endl;
0277     exit(EXIT_FAILURE);
0278   }
0279   LumiFile += Year;
0280   LumiFile += ".txt";
0281   return LumiFile;
0282 };
0283 
0284 /*! \fn runlistfromlumifile
0285  *  \brief Get a vector containing the list of runs for which the luminosity is known.
0286  */
0287 
0288 vector<int> runlistfromlumifile(TString Year) {
0289   TString filename = lumifileperyear(Year, "run");
0290   fs::path path(filename.Data());
0291   if (fs::is_empty(path)) {
0292     cout << "ERROR: Empty file " << path.c_str() << endl;
0293     exit(EXIT_FAILURE);
0294   }
0295   TGraph *scale = new TGraph(filename.Data());
0296   double *xscale = scale->GetX();
0297   size_t N = scale->GetN();
0298   vector<int> runs;
0299   for (size_t i = 0; i < N; i++)
0300     runs.push_back(xscale[i]);
0301   return runs;
0302 }
0303 
0304 /*! \fn checkrunlist
0305  *  \brief Check whether all runs of interest are present in the luminosity per run txt file and whether all IOVs analized have been correctly processed
0306  */
0307 
0308 bool checkrunlist(vector<int> runs, vector<int> IOVlist, TString Year) {
0309   vector<int> runlist = runlistfromlumifile(Year);
0310   vector<int> missingruns;  //runs for which the luminosity is not found
0311   vector<int> lostruns;     //IOVs for which the DMR were not found
0312   bool problemfound = false;
0313   for (int run : IOVlist) {
0314     if (find(runlist.begin(), runlist.end(), run) == runlist.end()) {
0315       problemfound = true;
0316       missingruns.push_back(run);
0317     }
0318   }
0319   if (!IOVlist.empty())
0320     for (int IOV : IOVlist) {
0321       if (find(runs.begin(), runs.end(), IOV) == runs.end()) {
0322         problemfound = true;
0323         lostruns.push_back(IOV);
0324       }
0325     }
0326   std::sort(missingruns.begin(), missingruns.end());
0327   if (problemfound) {
0328     if (!lostruns.empty()) {
0329       cout << "WARNING: some IOVs where not found among the list of available DMRs" << endl
0330            << "List of missing IOVs:" << endl;
0331       for (int lostrun : lostruns)
0332         cout << to_string(lostrun) << " ";
0333       cout << endl;
0334     }
0335     if (!missingruns.empty()) {
0336       cout << "WARNING: some IOVs are missing in the run/luminosity txt file" << endl
0337            << "List of missing runs:" << endl;
0338       for (int missingrun : missingruns)
0339         cout << to_string(missingrun) << " ";
0340       cout << endl;
0341     }
0342   }
0343   return problemfound;
0344 }
0345 
0346 /*! \fn DMRtrends
0347  *  \brief Create and plot the DMR trends.
0348  */
0349 
0350 void DMRtrends(vector<int> IOVlist,
0351                vector<string> Variables,
0352                vector<string> labels,
0353                TString Year,
0354                string myValidation,
0355                vector<string> geometries,
0356                vector<Color_t> colours,
0357                TString outputdir,
0358                bool pixelupdate,
0359                vector<int> pixelupdateruns,
0360                bool showlumi,
0361                bool FORCE) {
0362   fs::path path(outputdir.Data());
0363   if (!(fs::exists(path))) {
0364     cout << "WARNING: Output directory (" << outputdir.Data() << ") not found, it will be created automatically!"
0365          << endl;
0366     //  exit(EXIT_FAILURE);
0367     fs::create_directory(path);
0368     if (!(fs::exists(path))) {
0369       cout << "ERROR: Output directory (" << outputdir.Data() << ") has not been created!" << endl
0370            << "At least the parent directory needs to exist, please check!" << endl;
0371       exit(EXIT_FAILURE);
0372     }
0373   }
0374   for (const auto &Variable : Variables) {
0375     compileDMRTrends(IOVlist, Variable, labels, Year, myValidation, geometries, showlumi, FORCE);
0376     cout << "Begin plotting" << endl;
0377     PlotDMRTrends(IOVlist,
0378                   Variable,
0379                   labels,
0380                   Year,
0381                   myValidation,
0382                   geometries,
0383                   colours,
0384                   outputdir,
0385                   pixelupdate,
0386                   pixelupdateruns,
0387                   showlumi);
0388   }
0389 };
0390 
0391 /*! \fn compileDMRTrends
0392  *  \brief  Create a file where the DMR trends are stored in the form of TGraph.
0393  */
0394 
0395 void compileDMRTrends(vector<int> IOVlist,
0396                       TString Variable,
0397                       vector<string> labels,
0398                       TString Year,
0399                       string myValidation,
0400                       vector<string> geometries,
0401                       bool showlumi,
0402                       bool FORCE) {
0403   gROOT->SetBatch();
0404   vector<int> RunNumbers;
0405   vector<TString> filenames;
0406   TRegexp regexp("[0-9][0-9][0-9][0-9][0-9][0-9]");
0407   for (const auto &entry : fs::recursive_directory_iterator(myValidation)) {
0408     bool found_all_labels = true;
0409     for (const string &label : labels) {
0410       if (entry.path().string().find(label) == std::string::npos)
0411         found_all_labels = false;
0412     }
0413     if ((entry.path().string().find("ExtendedOfflineValidation_Images/OfflineValidationSummary.root") !=
0414          std::string::npos) &&
0415         found_all_labels) {
0416       if (fs::is_empty(entry.path()))
0417         cout << "ERROR: Empty file " << entry.path() << endl;
0418       else {
0419         TString filename(entry.path().string());
0420         filenames.push_back(filename);
0421         TString runstring(filename(regexp));
0422         if (runstring.IsFloat()) {
0423           int runN = runstring.Atoi();
0424           RunNumbers.push_back(runN);
0425         }
0426       }
0427     }
0428   }
0429   if (RunNumbers.empty()) {
0430     cout << "ERROR: No available DMRs found!" << endl
0431          << "Please check that the OfflineValidationSummary.root file is in any directory where the DMR per IOV have "
0432             "been stored!"
0433          << endl;
0434     exit(EXIT_FAILURE);
0435   } else if (checkrunlist(RunNumbers, IOVlist, Year)) {
0436     cout << "Please check the DMRs that have been produced!" << endl;
0437     if (!FORCE)
0438       exit(EXIT_FAILURE);
0439   }
0440 
0441   vector<TString> structures{"BPIX", "BPIX_y", "FPIX", "FPIX_y", "TIB", "TID", "TOB", "TEC"};
0442 
0443   const map<TString, int> nlayers = numberOfLayers(Year);
0444 
0445   float ScaleFactor = DMRFactor;
0446   if (Variable == "DrmsNR")
0447     ScaleFactor = 1;
0448 
0449   map<pair<pair<TString, int>, TString>, Geometry> mappoints;  // pair = (structure, layer), geometry
0450 
0451   std::sort(filenames.begin(), filenames.end());  //order the files in alphabetical order
0452   for (const TString &filename : filenames) {
0453     int runN;
0454     TString runstring(filename(regexp));
0455     if (runstring.IsFloat()) {
0456       runN = runstring.Atoi();
0457     } else {
0458       cout << "ERROR: run number not retrieved for file " << filename << endl;
0459       continue;
0460     }
0461 
0462     TFile *f = new TFile(filename, "READ");
0463 
0464     for (TString &structure : structures) {
0465       TString structname = structure;
0466       structname.ReplaceAll("_y", "");
0467       size_t layersnumber = nlayers.at(structname);
0468       for (size_t layer = 0; layer <= layersnumber; layer++) {
0469         for (const string &geometry : geometries) {
0470           TString name = Variable + "_" + getName(structure, layer, geometry);
0471           TH1F *histo = dynamic_cast<TH1F *>(f->Get(name));
0472           //Geometry *geom =nullptr;
0473           Point *point = nullptr;
0474           // Three possibilities:
0475           //  - All histograms are produced correctly
0476           //  - Only the non-split histograms are produced
0477           //  - No histogram is produced correctly
0478           //  FORCE means that the Point is not added to the points collection in the chosen geometry for that structure
0479           //  If FORCE is not enabled a default value for the Point is used (-9999) which will appear in the plots
0480           if (!histo) {
0481             cout << "Run" << runN << " Histogram: " << name << " not found" << endl;
0482             if (FORCE)
0483               continue;
0484             point = new Point(runN, ScaleFactor);
0485           } else if (structure != "TID" && structure != "TEC") {
0486             TH1F *histoplus = dynamic_cast<TH1F *>(f->Get((name + "_plus")));
0487             TH1F *histominus = dynamic_cast<TH1F *>(f->Get((name + "_minus")));
0488             if (!histoplus || !histominus) {
0489               cout << "Run" << runN << " Histogram: " << name << " plus or minus not found" << endl;
0490               if (FORCE)
0491                 continue;
0492               point = new Point(runN, ScaleFactor, histo);
0493             } else
0494               point = new Point(runN, ScaleFactor, histo, histoplus, histominus);
0495 
0496           } else
0497             point = new Point(runN, ScaleFactor, histo);
0498           mappoints[make_pair(make_pair(structure, layer), geometry)].points.push_back(*point);
0499         }
0500       }
0501     }
0502     f->Close();
0503   }
0504   TString outname = myValidation + "DMRtrends";
0505   for (const auto &label : labels) {
0506     outname += "_";
0507     outname += label;
0508   }
0509   outname += ".root";
0510   cout << outname << endl;
0511   TFile *fout = TFile::Open(outname, "RECREATE");
0512   for (TString &structure : structures) {
0513     TString structname = structure;
0514     structname.ReplaceAll("_y", "");
0515     size_t layersnumber = nlayers.at(structname);
0516     for (size_t layer = 0; layer <= layersnumber; layer++) {
0517       for (const string &geometry : geometries) {
0518         TString name = Variable + "_" + getName(structure, layer, geometry);
0519         Geometry geom = mappoints[make_pair(make_pair(structure, layer), geometry)];
0520         using Trend = vector<float> (Geometry::*)() const;
0521         vector<Trend> trends{&Geometry::Mu,
0522                              &Geometry::Sigma,
0523                              &Geometry::MuPlus,
0524                              &Geometry::SigmaPlus,
0525                              &Geometry::MuMinus,
0526                              &Geometry::SigmaMinus,
0527                              &Geometry::DeltaMu,
0528                              &Geometry::SigmaDeltaMu};
0529         vector<TString> variables{
0530             "mu", "sigma", "muplus", "sigmaplus", "muminus", "sigmaminus", "deltamu", "sigmadeltamu"};
0531         vector<float> runs = geom.Run();
0532         size_t n = runs.size();
0533         vector<float> emptyvec;
0534         for (size_t i = 0; i < runs.size(); i++)
0535           emptyvec.push_back(0.);
0536         for (size_t iVar = 0; iVar < variables.size(); iVar++) {
0537           Trend trend = trends.at(iVar);
0538           TGraphErrors *g = new TGraphErrors(n, runs.data(), (geom.*trend)().data(), emptyvec.data(), emptyvec.data());
0539           g->SetTitle(geometry.c_str());
0540           g->Write(name + "_" + variables.at(iVar));
0541         }
0542         vector<pair<Trend, Trend>> trendspair{make_pair(&Geometry::Mu, &Geometry::Sigma),
0543                                               make_pair(&Geometry::MuPlus, &Geometry::SigmaPlus),
0544                                               make_pair(&Geometry::MuMinus, &Geometry::SigmaMinus),
0545                                               make_pair(&Geometry::DeltaMu, &Geometry::SigmaDeltaMu)};
0546         vector<pair<TString, TString>> variablepairs{make_pair("mu", "sigma"),
0547                                                      make_pair("muplus", "sigmaplus"),
0548                                                      make_pair("muminus", "sigmaminus"),
0549                                                      make_pair("deltamu", "sigmadeltamu")};
0550         for (size_t iVar = 0; iVar < variablepairs.size(); iVar++) {
0551           Trend meantrend = trendspair.at(iVar).first;
0552           Trend sigmatrend = trendspair.at(iVar).second;
0553           TGraphErrors *g = new TGraphErrors(
0554               n, runs.data(), (geom.*meantrend)().data(), emptyvec.data(), (geom.*sigmatrend)().data());
0555           g->SetTitle(geometry.c_str());
0556           TString graphname = name + "_" + variablepairs.at(iVar).first;
0557           graphname += variablepairs.at(iVar).second;
0558           g->Write(graphname);
0559         }
0560       }
0561     }
0562   }
0563   fout->Close();
0564 }
0565 
0566 /*! \fn PixelUpdateLines
0567  *  \brief  Adds to the canvas vertical lines corresponding to the pixelupdateruns
0568  */
0569 void PixelUpdateLines(TCanvas *c, TString Year, bool showlumi, vector<int> pixelupdateruns) {
0570   vector<TPaveText *> labels;
0571   double lastlumi = 0.;
0572   c->cd();
0573   size_t index = 0;
0574   for (int pixelupdaterun : pixelupdateruns) {
0575     double lumi = 0.;
0576     if (showlumi)
0577       lumi = getintegratedlumiuptorun(
0578           pixelupdaterun,
0579           Year);  //The vertical line needs to be drawn at the beginning of the run where the pixel update was implemented, thus only the integrated luminosity up to that run is required.
0580     else
0581       lumi = pixelupdaterun;
0582     TLine *line = new TLine(lumi, c->GetUymin(), lumi, c->GetUymax());
0583     line->SetLineColor(kBlue);
0584     line->SetLineStyle(9);
0585     line->Draw();
0586     //Due to the way the coordinates within the Canvas are set, the following steps are required to draw the TPaveText:
0587     // Compute the gPad coordinates in TRUE normalized space (NDC)
0588     int ix1;
0589     int ix2;
0590     int iw = gPad->GetWw();
0591     int ih = gPad->GetWh();
0592     double x1p, y1p, x2p, y2p;
0593     gPad->GetPadPar(x1p, y1p, x2p, y2p);
0594     ix1 = (Int_t)(iw * x1p);
0595     ix2 = (Int_t)(iw * x2p);
0596     double wndc = TMath::Min(1., (double)iw / (double)ih);
0597     double rw = wndc / (double)iw;
0598     double x1ndc = (double)ix1 * rw;
0599     double x2ndc = (double)ix2 * rw;
0600     // Ratios to convert user space in TRUE normalized space (NDC)
0601     double rx1, ry1, rx2, ry2;
0602     gPad->GetRange(rx1, ry1, rx2, ry2);
0603     double rx = (x2ndc - x1ndc) / (rx2 - rx1);
0604     double _sx;
0605     // Left limit of the TPaveText
0606     _sx = rx * (lumi - rx1) + x1ndc;
0607     // To avoid an overlap between the TPaveText a vertical shift is done when the IOVs are too close
0608     if (_sx < lastlumi) {
0609       index++;
0610     } else
0611       index = 0;
0612     TPaveText *box = new TPaveText(_sx + 0.0015, 0.86 - index * 0.04, _sx + 0.035, 0.89 - index * 0.04, "blNDC");
0613     box->SetFillColor(10);
0614     box->SetBorderSize(1);
0615     box->SetLineColor(kBlack);
0616     TText *textRun = box->AddText(Form("%i", int(pixelupdaterun)));
0617     textRun->SetTextSize(0.025);
0618     labels.push_back(box);
0619     lastlumi = _sx + 0.035;
0620 
0621     gPad->RedrawAxis();
0622   }
0623   //Drawing in a separate loop to ensure that the labels are drawn on top of the lines
0624   for (auto label : labels) {
0625     label->Draw("same");
0626   }
0627   c->Update();
0628 }
0629 
0630 /*! \fn getintegratedlumiuptorun
0631  *  \brief Returns the integrated luminosity up to the run of interest
0632  */
0633 
0634 double getintegratedlumiuptorun(int run, TString Year, double min) {
0635   TGraph *scale = new TGraph((lumifileperyear(Year, "run")).Data());
0636   int Nscale = scale->GetN();
0637   double *xscale = scale->GetX();
0638   double *yscale = scale->GetY();
0639 
0640   double lumi = min;
0641   int index = -1;
0642   for (int j = 0; j < Nscale; j++) {
0643     lumi += yscale[j];
0644     if (run >= xscale[j]) {
0645       index = j;
0646       continue;
0647     }
0648   }
0649   lumi = min;
0650   for (int j = 0; j < index; j++)
0651     lumi += yscale[j] / lumiFactor;
0652 
0653   return lumi;
0654 }
0655 /*! \fn scalebylumi
0656  *  \brief Scale X-axis of the TGraph and the error on that axis according to the integrated luminosity.
0657  */
0658 
0659 void scalebylumi(TGraphErrors *g, vector<pair<int, double>> lumiIOVpairs) {
0660   size_t N = g->GetN();
0661   vector<double> x, y, xerr, yerr;
0662 
0663   //TGraph * scale = new TGraph((lumifileperyear(Year,"IOV")).Data());
0664   size_t Nscale = lumiIOVpairs.size();
0665 
0666   size_t i = 0;
0667   while (i < N) {
0668     double run, yvalue;
0669     g->GetPoint(i, run, yvalue);
0670     size_t index = -1;
0671     for (size_t j = 0; j < Nscale; j++) {
0672       if (run == (lumiIOVpairs.at(j)
0673                       .first)) {  //If the starting run of an IOV is included in the list of IOVs, the index is stored
0674         index = j;
0675         continue;
0676       } else if (run > (lumiIOVpairs.at(j).first))
0677         continue;
0678     }
0679     if (lumiIOVpairs.at(index).second == 0 || index < 0.) {
0680       N = N - 1;
0681       g->RemovePoint(i);
0682     } else {
0683       double xvalue = 0.;
0684       for (size_t j = 0; j < index; j++)
0685         xvalue += lumiIOVpairs.at(j).second / lumiFactor;
0686       x.push_back(xvalue + (lumiIOVpairs.at(index).second / (lumiFactor * 2.)));
0687       if (yvalue <= DUMMY) {
0688         y.push_back(DUMMY);
0689         yerr.push_back(0.);
0690       } else {
0691         y.push_back(yvalue);
0692         yerr.push_back(g->GetErrorY(i));
0693       }
0694       xerr.push_back(lumiIOVpairs.at(index).second / (lumiFactor * 2.));
0695       i = i + 1;
0696     }
0697   }
0698   g->GetHistogram()->Delete();
0699   g->SetHistogram(nullptr);
0700   for (size_t i = 0; i < N; i++) {
0701     g->SetPoint(i, x.at(i), y.at(i));
0702     g->SetPointError(i, xerr.at(i), yerr.at(i));
0703   }
0704 }
0705 
0706 /*! \fn lumiperIOV
0707  *  \brief Retrieve luminosity per IOV
0708  */
0709 
0710 vector<pair<int, double>> lumiperIOV(vector<int> IOVlist, TString Year) {
0711   size_t N = IOVlist.size();
0712   vector<pair<int, double>> lumiperIOV;
0713   TGraph *scale = new TGraph((lumifileperyear(Year, "run")).Data());
0714   size_t Nscale = scale->GetN();
0715   double *xscale = scale->GetX();
0716   double *yscale = scale->GetY();
0717 
0718   size_t i = 0;
0719   size_t index = 0;
0720   while (i <= N) {
0721     double run = 0;
0722     double lumi = 0.;
0723     if (i != N)
0724       run = IOVlist.at(i);
0725     else
0726       run = 0;
0727     for (size_t j = index; j < Nscale; j++) {
0728       if (run == xscale[j]) {
0729         index = j;
0730         break;
0731       } else
0732         lumi += yscale[j];
0733     }
0734     if (i == 0)
0735       lumiperIOV.push_back(make_pair(0, lumi));
0736     else
0737       lumiperIOV.push_back(make_pair(IOVlist.at(i - 1), lumi));
0738     ++i;
0739   }
0740   //for debugging:
0741   double lumiInput = 0;
0742   double lumiOutput = 0.;
0743   for (size_t j = 0; j < Nscale; j++)
0744     lumiInput += yscale[j];
0745   //cout << "Total lumi: " << lumiInput <<endl;
0746   for (size_t j = 0; j < lumiperIOV.size(); j++)
0747     lumiOutput += lumiperIOV.at(j).second;
0748   //cout << "Total lumi saved for IOVs: " << lumiOutput <<endl;
0749   if (abs(lumiInput - lumiOutput) > 0.5) {
0750     cout << "ERROR: luminosity retrieved for IOVs does not match the one for the runs" << endl
0751          << "Please check that all IOV first runs are part of the run-per-lumi file!" << endl;
0752     exit(EXIT_FAILURE);
0753   }
0754   return lumiperIOV;
0755 }
0756 
0757 /*! \fn ConvertToHist
0758  *  \brief A TH1F is constructed using the points and the errors collected in the TGraphErrors
0759  */
0760 
0761 TH1F *ConvertToHist(TGraphErrors *g) {
0762   size_t N = g->GetN();
0763   double *x = g->GetX();
0764   double *y = g->GetY();
0765   double *xerr = g->GetEX();
0766   vector<float> bins;
0767   bins.push_back(x[0] - xerr[0]);
0768   for (size_t i = 1; i < N; i++) {
0769     if ((x[i - 1] + xerr[i - 1]) > (bins.back() + pow(10, -6)))
0770       bins.push_back(x[i - 1] + xerr[i - 1]);
0771     if ((x[i] - xerr[i]) > (bins.back() + pow(10, -6)))
0772       bins.push_back(x[i] - xerr[i]);
0773   }
0774   bins.push_back(x[N - 1] + xerr[N - 1]);
0775   TString histoname = "histo_";
0776   histoname += g->GetName();
0777   TH1F *histo = new TH1F(histoname, g->GetTitle(), bins.size() - 1, bins.data());
0778   for (size_t i = 0; i < N; i++) {
0779     histo->Fill(x[i], y[i]);
0780     histo->SetBinError(histo->FindBin(x[i]), pow(10, -6));
0781   }
0782   return histo;
0783 }
0784 
0785 /*! \fn PlotDMRTrends
0786  *  \brief Plot the DMR trends.
0787  */
0788 
0789 void PlotDMRTrends(vector<int> IOVlist,
0790                    TString Variable,
0791                    vector<string> labels,
0792                    TString Year,
0793                    string myValidation,
0794                    vector<string> geometries,
0795                    vector<Color_t> colours,
0796                    TString outputdir,
0797                    bool pixelupdate,
0798                    vector<int> pixelupdateruns,
0799                    bool showlumi) {
0800   gErrorIgnoreLevel = kWarning;
0801   checkrunlist(pixelupdateruns, {}, Year);
0802   vector<TString> structures{"BPIX", "BPIX_y", "FPIX", "FPIX_y", "TIB", "TID", "TOB", "TEC"};
0803 
0804   const map<TString, int> nlayers = numberOfLayers(Year);
0805 
0806   vector<pair<int, double>> lumiIOVpairs;
0807   if (showlumi)
0808     lumiIOVpairs = lumiperIOV(IOVlist, Year);
0809 
0810   TString filename = myValidation + "DMRtrends";
0811   for (const auto &label : labels) {
0812     filename += "_";
0813     filename += label;
0814   }
0815   filename += ".root";
0816   cout << filename << endl;
0817   TFile *in = new TFile(filename);
0818   for (TString &structure : structures) {
0819     TString structname = structure;
0820     structname.ReplaceAll("_y", "");
0821     int layersnumber = nlayers.at(structname);
0822     for (int layer = 0; layer <= layersnumber; layer++) {
0823       vector<TString> variables{"mu",
0824                                 "sigma",
0825                                 "muplus",
0826                                 "sigmaplus",
0827                                 "muminus",
0828                                 "sigmaminus",
0829                                 "deltamu",
0830                                 "sigmadeltamu",
0831                                 "musigma",
0832                                 "muplussigmaplus",
0833                                 "muminussigmaminus",
0834                                 "deltamusigmadeltamu"};
0835       vector<string> YaxisNames{
0836           "#mu [#mum]",
0837           "#sigma_{#mu} [#mum]",
0838           "#mu outward [#mum]",
0839           "#sigma_{#mu outward} [#mum]",
0840           "#mu inward [#mum]",
0841           "#sigma_{#mu inward} [#mum]",
0842           "#Delta#mu [#mum]",
0843           "#sigma_{#Delta#mu} [#mum]",
0844           "#mu [#mum]",
0845           "#mu outward [#mum]",
0846           "#mu inward [#mum]",
0847           "#Delta#mu [#mum]",
0848       };
0849       if (Variable == "DrmsNR")
0850         YaxisNames = {
0851             "RMS(x'_{pred}-x'_{hit} /#sigma)",
0852             "#sigma_{RMS(x'_{pred}-x'_{hit} /#sigma)}",
0853             "RMS(x'_{pred}-x'_{hit} /#sigma) outward",
0854             "#sigma_{#mu outward}",
0855             "RMS(x'_{pred}-x'_{hit} /#sigma) inward",
0856             "#sigma_{RMS(x'_{pred}-x'_{hit} /#sigma) inward}",
0857             "#DeltaRMS(x'_{pred}-x'_{hit} /#sigma)",
0858             "#sigma_{#DeltaRMS(x'_{pred}-x'_{hit} /#sigma)}",
0859             "RMS(x'_{pred}-x'_{hit} /#sigma)",
0860             "RMS(x'_{pred}-x'_{hit} /#sigma) outward",
0861             "RMS(x'_{pred}-x'_{hit} /#sigma) inward",
0862             "#DeltaRMS(x'_{pred}-x'_{hit} /#sigma)",
0863         };
0864       //For debugging purposes we still might want to have a look at plots for a variable without errors, once ready for the PR those variables will be removed and the iterator will start from 0
0865       for (size_t i = 0; i < variables.size(); i++) {
0866         TString variable = variables.at(i);
0867         if (variable.Contains("plus") || variable.Contains("minus") || variable.Contains("delta")) {
0868           if (structname == "TEC" || structname == "TID")
0869             continue;  //Lorentz drift cannot appear in TEC and TID. These structures are skipped when looking at outward and inward pointing modules.
0870         }
0871         TCanvas *c = new TCanvas("dummy", "", 2000, 800);
0872 
0873         vector<Color_t>::iterator colour = colours.begin();
0874 
0875         TMultiGraph *mg = new TMultiGraph(structure, structure);
0876         THStack *mh = new THStack(structure, structure);
0877         size_t igeom = 0;
0878         for (const string &geometry : geometries) {
0879           TString name = Variable + "_" + getName(structure, layer, geometry);
0880           TGraphErrors *g = dynamic_cast<TGraphErrors *>(in->Get(name + "_" + variables.at(i)));
0881           g->SetName(name + "_" + variables.at(i));
0882           if (i >= 8) {
0883             g->SetLineWidth(1);
0884             g->SetLineColor(*colour);
0885             g->SetFillColorAlpha(*colour, 0.2);
0886           }
0887           vector<vector<double>> vectors;
0888           //if(showlumi&&i<8)scalebylumi(dynamic_cast<TGraph*>(g));
0889           if (showlumi)
0890             scalebylumi(g, lumiIOVpairs);
0891           g->SetLineColor(*colour);
0892           g->SetMarkerColor(*colour);
0893           TH1F *h = ConvertToHist(g);
0894           h->SetLineColor(*colour);
0895           h->SetMarkerColor(*colour);
0896           h->SetMarkerSize(0);
0897           h->SetLineWidth(1);
0898 
0899           if (i < 8) {
0900             mg->Add(g, "PL");
0901             mh->Add(h, "E");
0902           } else {
0903             mg->Add(g, "2");
0904             mh->Add(h, "E");
0905           }
0906           ++colour;
0907           ++igeom;
0908         }
0909 
0910         gStyle->SetOptTitle(0);
0911         gStyle->SetPadLeftMargin(0.08);
0912         gStyle->SetPadRightMargin(0.05);
0913         gPad->SetTickx();
0914         gPad->SetTicky();
0915         gStyle->SetLegendTextSize(0.025);
0916 
0917         double max = 6;
0918         double min = -4;
0919         if (Variable == "DrmsNR") {
0920           if (variable.Contains("delta")) {
0921             max = 0.15;
0922             min = -0.1;
0923           } else {
0924             max = 1.2;
0925             min = 0.6;
0926           }
0927         }
0928         double range = max - min;
0929 
0930         if (((variable == "sigma" || variable == "sigmaplus" || variable == "sigmaminus" ||
0931               variable == "sigmadeltamu") &&
0932              range >= 2)) {
0933           mg->SetMaximum(4);
0934           mg->SetMinimum(-2);
0935         } else {
0936           mg->SetMaximum(max + range * 0.1);
0937           mg->SetMinimum(min - range * 0.3);
0938         }
0939 
0940         if (i < 8) {
0941           mg->Draw("a");
0942         } else {
0943           mg->Draw("a2");
0944         }
0945 
0946         char *Ytitle = (char *)YaxisNames.at(i).c_str();
0947         mg->GetYaxis()->SetTitle(Ytitle);
0948         mg->GetXaxis()->SetTitle(showlumi ? "Integrated lumi [1/fb]" : "IOV number");
0949         mg->GetXaxis()->CenterTitle(true);
0950         mg->GetYaxis()->CenterTitle(true);
0951         mg->GetYaxis()->SetTitleOffset(.5);
0952         mg->GetYaxis()->SetTitleSize(.05);
0953         mg->GetXaxis()->SetTitleSize(.04);
0954         if (showlumi)
0955           mg->GetXaxis()->SetLimits(0., mg->GetXaxis()->GetXmax());
0956 
0957         c->Update();
0958 
0959         TLegend *legend = c->BuildLegend();
0960         // TLegend *legend = c->BuildLegend(0.15,0.18,0.15,0.18);
0961         int Ngeom = geometries.size();
0962         if (Ngeom >= 6)
0963           legend->SetNColumns(2);
0964         else if (Ngeom >= 9)
0965           legend->SetNColumns(3);
0966         else
0967           legend->SetNColumns(1);
0968         TString structtitle = "#bf{";
0969         if (structure.Contains("PIX") && !(structure.Contains("_y")))
0970           structtitle += structure + " (x)";
0971         else if (structure.Contains("_y")) {
0972           TString substring(structure(0, 4));
0973           structtitle += substring + " (y)";
0974         } else
0975           structtitle += structure;
0976         if (layer != 0) {
0977           if (structure == "TID" || structure == "TEC" || structure == "FPIX" || structure == "FPIX_y")
0978             structtitle += "  disc ";
0979           else
0980             structtitle += "  layer ";
0981           structtitle += layer;
0982         }
0983         structtitle += "}";
0984         PixelUpdateLines(c, Year, showlumi, pixelupdateruns);
0985 
0986         TPaveText *CMSworkInProgress = new TPaveText(
0987             0, mg->GetYaxis()->GetXmax() + range * 0.02, 2.5, mg->GetYaxis()->GetXmax() + range * 0.12, "nb");
0988         CMSworkInProgress->AddText("#scale[1.1]{CMS} #bf{Internal}");
0989         CMSworkInProgress->SetTextAlign(12);
0990         CMSworkInProgress->SetTextSize(0.04);
0991         CMSworkInProgress->SetFillColor(10);
0992         CMSworkInProgress->Draw();
0993         TPaveText *TopRightCorner = new TPaveText(0.95 * (mg->GetXaxis()->GetXmax()),
0994                                                   mg->GetYaxis()->GetXmax() + range * 0.02,
0995                                                   (mg->GetXaxis()->GetXmax()),
0996                                                   mg->GetYaxis()->GetXmax() + range * 0.12,
0997                                                   "nb");
0998         TopRightCorner->AddText(Year + " pp collisions");
0999         TopRightCorner->SetTextAlign(32);
1000         TopRightCorner->SetTextSize(0.04);
1001         TopRightCorner->SetFillColor(10);
1002         TopRightCorner->Draw();
1003         TPaveText *structlabel = new TPaveText(0.95 * (mg->GetXaxis()->GetXmax()),
1004                                                mg->GetYaxis()->GetXmin() + range * 0.02,
1005                                                0.99 * (mg->GetXaxis()->GetXmax()),
1006                                                mg->GetYaxis()->GetXmin() + range * 0.12,
1007                                                "nb");
1008         structlabel->AddText(structtitle.Data());
1009         structlabel->SetTextAlign(32);
1010         structlabel->SetTextSize(0.04);
1011         structlabel->SetFillColor(10);
1012         structlabel->Draw();
1013 
1014         legend->Draw();
1015         mh->Draw("nostack same");
1016         c->Update();
1017         TString structandlayer = getName(structure, layer, "");
1018         TString printfile = outputdir;
1019         if (!(outputdir.EndsWith("/")))
1020           outputdir += "/";
1021         for (const auto &label : labels) {
1022           printfile += label;
1023           printfile += "_";
1024         }
1025         printfile += Variable;
1026         printfile += "_";
1027         printfile += variable + structandlayer;
1028         c->SaveAs(printfile + ".pdf");
1029         c->SaveAs(printfile + ".eps");
1030         c->SaveAs(printfile + ".png");
1031         c->Destructor();
1032       }
1033     }
1034   }
1035   in->Close();
1036 }
1037 
1038 /*! \fn main
1039  *  \brief main function: if no arguments are specified a default list of arguments is used, otherwise a total of 9 arguments are required:
1040  * @param IOVlist:                 string containing the list of IOVs separated by a ","
1041  * @param labels:                  string containing labels that must be part of the input files
1042  * @param Year:                    string containing the year of the studied runs (needed to retrieve the lumi-per-run file)
1043  * @param pathtoDMRs:              string containing the path to the directory where the DMRs are stored
1044  * @param geometrieandcolours:     string containing the list of geometries and colors in the following way name1:color1,name2:color2 etc.
1045  * @param outputdirectory:         string containing the output directory for the plots
1046  * @param pixelupdatelist:         string containing the list of pixelupdates separated by a ","
1047  * @param showpixelupdate:         boolean that if set to true will allow to plot vertical lines in the canvas corresponding to the pixel updates
1048  * @param showlumi:                boolean, if set to false the trends will be presented in function of the run (IOV) number, if set to true the luminosity is used on the x axis
1049  * @param FORCE:              //!< boolean, if set to true the plots will be made regardless of possible errors.
1050  *                                 Eventual errors while running the code will be ignored and just warnings will appear in the output.
1051  */
1052 
1053 int main(int argc, char *argv[]) {
1054   if (argc == 1) {
1055     vector<int> IOVlist = {314881, 315257, 315488, 315489, 315506, 316239, 316271, 316361, 316363, 316378, 316456,
1056                            316470, 316505, 316569, 316665, 316758, 317080, 317182, 317212, 317295, 317339, 317382,
1057                            317438, 317527, 317661, 317664, 318712, 319337, 319460, 320841, 320854, 320856, 320888,
1058                            320916, 320933, 320980, 321009, 321119, 321134, 321164, 321261, 321294, 321310, 321393,
1059                            321397, 321431, 321461, 321710, 321735, 321773, 321774, 321778, 321820, 321831, 321880,
1060                            321960, 322014, 322510, 322603, 323232, 323423, 323472, 323475, 323693, 323794, 323976,
1061                            324202, 324206, 324245, 324729, 324764, 324840, 324999, 325097, 325110};
1062     vector<int> pixelupdateruns{316758, 317527, 317661, 317664, 318227, 320377};  //2018
1063     //  vector<int> pixelupdateruns {290543, 297281, 298653, 299443, 300389, 301046, 302131, 303790, 303998, 304911, 305898};//2017
1064 
1065     cout << "WARNING: Running function with arguments specified in DMRtrends.cc" << endl
1066          << "If you want to specify the arguments from command line run the macro as follows:" << endl
1067          << "DMRtrends labels pathtoDMRs geometriesandcolourspairs outputdirectory showpixelupdate showlumi FORCE"
1068          << endl;
1069 
1070     //PLEASE READ: for debugging purposes please keep at least one example that works commented.
1071     //             Error messages are still a W.I.P. and having a working example available is useful for debugging.
1072     //Example provided for a currently working set of parameters:
1073     DMRtrends(IOVlist,
1074               {"median", "DrmsNR"},
1075               {"minbias"},
1076               "2018",
1077               "/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/results/adewit/UL18_DMR_MB_eoj_v3/",
1078               {"1-step hybrid", "mid18 rereco", "counter twist"},
1079               {kBlue, kBlack, kGreen + 3},
1080               "/afs/cern.ch/user/a/acardini/commonValidation/results/acardini/DMRs/DMRTrends/test/",
1081               true,
1082               pixelupdateruns,
1083               true,
1084               true);
1085     return 0;
1086   } else if (argc < 12) {
1087     cout << "DMRtrends IOVlist labels Year pathtoDMRs geometriesandcolourspairs outputdirectory pixelupdatelist "
1088             "showpixelupdate showlumi FORCE"
1089          << endl;
1090 
1091     return 1;
1092   }
1093 
1094   TString runlist = argv[1], all_variables = argv[2], all_labels = argv[3], Year = argv[4], pathtoDMRs = argv[5],
1095           geometrieandcolours = argv[6],  //name1:title1:color1,name2:title2:color2,name3:title3:color3
1096       outputdirectory = argv[7], pixelupdatelist = argv[8];
1097   bool showpixelupdate = argv[9], showlumi = argv[10], FORCE = argv[11];
1098   TObjArray *vararray = all_variables.Tokenize(",");
1099   vector<string> Variables;
1100   Variables.reserve(vararray->GetEntries());
1101   for (int i = 0; i < vararray->GetEntries(); i++)
1102     Variables.push_back((string)(vararray->At(i)->GetName()));
1103   TObjArray *labelarray = all_labels.Tokenize(",");
1104   vector<string> labels;
1105   labels.reserve(labelarray->GetEntries());
1106   for (int i = 0; i < labelarray->GetEntries(); i++)
1107     labels.push_back((string)(labelarray->At(i)->GetName()));
1108   TObjArray *IOVarray = runlist.Tokenize(",");
1109   vector<int> IOVlist;
1110   IOVlist.reserve(IOVarray->GetEntries());
1111   for (int i = 0; i < IOVarray->GetEntries(); i++)
1112     IOVlist.push_back(stoi(IOVarray->At(i)->GetName()));
1113   vector<int> pixelupdateruns;
1114   TObjArray *PIXarray = pixelupdatelist.Tokenize(",");
1115   pixelupdateruns.reserve(PIXarray->GetEntries());
1116   for (int i = 0; i < PIXarray->GetEntries(); i++)
1117     pixelupdateruns.push_back(stoi(PIXarray->At(i)->GetName()));
1118   vector<string> geometries;
1119   //TO DO: the color is not taken correctly from command line
1120   vector<Color_t> colours;
1121   TObjArray *geometrieandcolourspairs = geometrieandcolours.Tokenize(",");
1122   for (int i = 0; i < geometrieandcolourspairs->GetEntries(); i++) {
1123     TObjArray *geomandcolourvec = TString(geometrieandcolourspairs->At(i)->GetName()).Tokenize(":");
1124     geometries.push_back(geomandcolourvec->At(0)->GetName());
1125     colours.push_back(ColorParser(geomandcolourvec->At(1)->GetName()));
1126   }
1127   DMRtrends(IOVlist,
1128             Variables,
1129             labels,
1130             Year,
1131             pathtoDMRs.Data(),
1132             geometries,
1133             colours,
1134             outputdirectory.Data(),
1135             showpixelupdate,
1136             pixelupdateruns,
1137             showlumi,
1138             FORCE);
1139 
1140   return 0;
1141 }