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
0033
0034 #define DUMMY -999.
0035
0036
0037
0038 #define lumiFactor 1000.
0039
0040
0041
0042 #define DMRFactor 10000.
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057 struct Point {
0058 float run, scale, mu, sigma, muplus, muminus, sigmaplus, sigmaminus;
0059
0060
0061
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
0074
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
0087
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
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
0178
0179
0180
0181
0182
0183 class Geometry {
0184 public:
0185 vector<Point> points;
0186
0187 private:
0188
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
0219
0220 };
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230
0231
0232
0233
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
0251
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
0262
0263
0264
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
0285
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
0305
0306
0307
0308 bool checkrunlist(vector<int> runs, vector<int> IOVlist, TString Year) {
0309 vector<int> runlist = runlistfromlumifile(Year);
0310 vector<int> missingruns;
0311 vector<int> lostruns;
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
0347
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
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
0392
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;
0450
0451 std::sort(filenames.begin(), filenames.end());
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
0473 Point *point = nullptr;
0474
0475
0476
0477
0478
0479
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
0567
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);
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
0587
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
0601 double rx1, ry1, rx2, ry2;
0602 gPad->GetRange(rx1, ry1, rx2, ry2);
0603 double rx = (x2ndc - x1ndc) / (rx2 - rx1);
0604 double _sx;
0605
0606 _sx = rx * (lumi - rx1) + x1ndc;
0607
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
0624 for (auto label : labels) {
0625 label->Draw("same");
0626 }
0627 c->Update();
0628 }
0629
0630
0631
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
0656
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
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)) {
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
0707
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
0741 double lumiInput = 0;
0742 double lumiOutput = 0.;
0743 for (size_t j = 0; j < Nscale; j++)
0744 lumiInput += yscale[j];
0745
0746 for (size_t j = 0; j < lumiperIOV.size(); j++)
0747 lumiOutput += lumiperIOV.at(j).second;
0748
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
0758
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
0786
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
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;
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
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
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
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
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};
1063
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
1071
1072
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],
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
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 }