CMS 3D CMS Logo

DMRtrends.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <string>
3 #include <vector>
4 #include <algorithm>
5 #include <map>
6 #include <iomanip>
7 #include <experimental/filesystem>
8 #include "TPad.h"
9 #include "TCanvas.h"
10 #include "TGraph.h"
11 #include "TGraphErrors.h"
12 #include "TMultiGraph.h"
13 #include "TH1.h"
14 #include "THStack.h"
15 #include "TROOT.h"
16 #include "TFile.h"
17 #include "TColor.h"
18 #include "TLegend.h"
19 #include "TLegendEntry.h"
20 #include "TMath.h"
21 #include "TRegexp.h"
22 #include "TPaveLabel.h"
23 #include "TPaveText.h"
24 #include "TStyle.h"
25 #include "TLine.h"
26 
27 using namespace std;
28 namespace fs = std::experimental::filesystem;
29 
33 #define DUMMY -999.
34 
37 #define lumiFactor 1000.
38 
41 #define DMRFactor 10000.
42 
55 struct Point {
56  float run, mu, sigma, muplus, muminus, sigmaplus, sigmaminus;
57 
61  Point (float Run = DUMMY, float y1 = DUMMY,float y2 = DUMMY, float y3 = DUMMY, float y4 = DUMMY, float y5 = DUMMY, float y6 = DUMMY) :
62  run(Run), mu(y1), sigma(y2), muplus(y3), muminus(y5), sigmaplus(y4), sigmaminus(y6)
63  {}
64 
68  Point (float Run, TH1 * histo, TH1 * histoplus, TH1 * histominus) :
69  Point(Run, histo->GetMean(), histo->GetMeanError(),
70  histoplus->GetMean(), histoplus->GetMeanError(),
71  histominus->GetMean(), histominus->GetMeanError())
72  {}
73 
77  Point (float Run, TH1 * histo) :
78  Point(Run, histo->GetMean(), histo->GetMeanError())
79  {}
80 
81  Point& operator=(const Point &p){ run = p.run; mu = p.mu; muplus = p.muplus; muminus = p.muminus; sigma = p.sigma; sigmaplus = p.sigmaplus; sigmaminus = p.sigmaminus; return *this;}
82 
83  float GetRun () const { return run ; }
84  float GetMu () const { return DMRFactor*mu ; }
85  float GetMuPlus () const { return DMRFactor*muplus ; }
86  float GetMuMinus () const { return DMRFactor*muminus ; }
87  float GetSigma () const { return DMRFactor*sigma ; }
88  float GetSigmaPlus () const { return DMRFactor*sigmaplus ; }
89  float GetSigmaMinus () const { return DMRFactor*sigmaminus; }
90  float GetDeltaMu () const { if(muplus==DUMMY&&muminus==DUMMY) return DUMMY;
91  else return DMRFactor*(muplus - muminus); }
92  float GetSigmaDeltaMu () const { if(sigmaplus==DUMMY&&sigmaminus==DUMMY) return DUMMY;
93  else return DMRFactor*hypot(sigmaplus,sigmaminus); };
94 };
95 
96 
100 
101 TString getName (TString structure, int layer, TString geometry);
102 TH1F* ConvertToHist(TGraphErrors *g);
103 vector<int> runlistfromlumifile(TString Year="2018");
104 bool checkrunlist(vector<int> runs, vector<int> IOVlist={}, TString Year="2018");
105 TString lumifileperyear(TString Year="2018", string RunOrIOV="IOV");
106 void scalebylumi(TGraphErrors *g, TString Year="2018", double min=0.);
107 double getintegratedlumiuptorun(int run, TString Year="2018", double min=0.);
108 void PixelUpdateLines(TCanvas *c, TString Year="2018", bool showlumi=false, vector<int>pixelupdateruns={314881, 316758, 317527, 318228, 320377});
109 void PlotDMRTrends(vector<string>labels={"MB"}, TString Year="2018", string myValidation="/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/results/acardini/DMRs/", vector<string> geometries={"GT","SG", "MP pix LBL","PIX HLS+ML STR fix"}, vector<Color_t> colours={kBlue, kRed, kGreen, kCyan}, TString outputdir="/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/alignmentObjects/acardini/DMRsTrends/", bool pixelupdate=false, vector<int> pixelupdateruns={314881, 316758, 317527, 318228, 320377}, bool showlumi=false);
110 void compileDMRTrends(vector<int> IOVlist, vector<string>labels={"MB"}, TString Year="2018", string myValidation="/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/results/acardini/DMRs/", vector<string> geometries={"GT","SG", "MP pix LBL","PIX HLS+ML STR fix"} ,bool showlumi=false, bool FORCE=false);
111 void DMRtrends(vector<int> IOVlist, vector<string>labels={"MB"}, TString Year="2018", string myValidation="/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/results/acardini/DMRs/", vector<string> geometries={"GT","SG", "MP pix LBL","PIX HLS+ML STR fix"}, vector<Color_t> colours={kBlue, kRed, kGreen, kCyan}, TString outputdir="/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/alignmentObjects/acardini/DMRsTrends/", bool pixelupdate=false, vector<int> pixelupdateruns={314881, 316758, 317527, 318228, 320377}, bool showlumi=false, bool FORCE=false);
112 
119 class Geometry {
120  public:
121  vector<Point> points;
122  private:
123  //template<typename T> vector<T> GetQuantity (T (Point::*getter)() const) const {
124  vector<float> GetQuantity (float (Point::*getter)() const) const {
125  vector<float> v;
126  for (Point point: points) {
127  float value = (point.*getter)();
128  v.push_back(value);
129  }
130  return v;
131  }
132  public:
133  TString title;
134  Geometry() : title ("") {}
135  Geometry(TString Title) : title(Title) {}
136  Geometry& operator=(const Geometry &geom){ title = geom.title; points = geom.points; return *this;}
137  void SetTitle(TString Title){ title = Title;}
138  TString GetTitle() { return title; }
139  vector<float> Run () const { return GetQuantity( &Point::GetRun ); }
140  vector<float> Mu () const { return GetQuantity( &Point::GetMu ); }
141  vector<float> MuPlus () const { return GetQuantity( &Point::GetMuPlus ); }
142  vector<float> MuMinus () const { return GetQuantity( &Point::GetMuMinus ); }
143  vector<float> Sigma () const { return GetQuantity( &Point::GetSigma ); }
144  vector<float> SigmaPlus () const { return GetQuantity( &Point::GetSigmaPlus ); }
145  vector<float> SigmaMinus () const { return GetQuantity( &Point::GetSigmaMinus ); }
146  vector<float> DeltaMu () const { return GetQuantity( &Point::GetDeltaMu ); }
147  vector<float> SigmaDeltaMu () const { return GetQuantity( &Point::GetSigmaDeltaMu); }
148  //vector<float> Graph (string variable) const {
149  // };
150 };
151 
153 //struct Layer {
154 // map<string,Geometry> geometries;
155 //};
156 //
157 //struct HLS {
158 // vector<Layer> layers;
159 // map<string,Geometry> geometries;
160 //};
161 
162 
167 TString getName (TString structure, int layer, TString geometry){
168  geometry.ReplaceAll(" ","_");
169  TString name=geometry+"_"+structure;
170  if(layer!=0){
171  if(structure=="TID"||structure=="TEC")name+="_disc";
172  else name+="_layer";
173  name+=layer;
174  }
175 
176  return name;
177 };
179 
183 TString lumifileperyear(TString Year, string RunOrIOV){
184  TString LumiFile=getenv("CMSSW_BASE"); LumiFile +="/src/Alignment/OfflineValidation/data/lumiper";
185  if(RunOrIOV!="run"&&RunOrIOV!="IOV") cout << "ERROR: Please specify \"run\" or \"IOV\" to retrieve the luminosity run by run or for each IOV"<<endl;
186  LumiFile+=RunOrIOV;
187  if(Year!="2017"&&Year!="2018") cout << "ERROR: Only 2017 and 2018 lumi-per-run/IOV are available, please check!"<<endl;
188  LumiFile+=Year;
189  LumiFile+=".txt";
190  return LumiFile;
191 
192 };
193 
198 vector<int> runlistfromlumifile(TString Year){
199  TGraph * scale = new TGraph((lumifileperyear(Year,"run")).Data());
200  double *xscale = scale->GetX();
201  size_t N = scale->GetN();
202  vector<int> runs;
203  for(size_t i=0;i<N;i++)runs.push_back(xscale[i]);
204  return runs;
205 }
206 
211 bool checkrunlist(vector<int> runs,vector<int> IOVlist, TString Year){
212  vector<int> runlist = runlistfromlumifile(Year);
213  vector<int> missingruns; //runs for which the luminosity is not found
214  vector<int> lostruns; //IOVs for which the DMR were not found
215  bool problemfound=false;
216  for(int run : runs){
217  if(find(runlist.begin(),runlist.end(),run)==runlist.end()){
218  problemfound=true;
219  missingruns.push_back(run);
220  }
221  }
222  if(!IOVlist.empty()) for(int IOV : IOVlist){
223  if(find(runs.begin(),runs.end(),IOV)==runs.end()){
224  problemfound=true;
225  lostruns.push_back(IOV);
226  }
227  }
228  std::sort(missingruns.begin(),missingruns.end());
229  if(problemfound){
230  if(!lostruns.empty()){
231  cout << "WARNING: some IOVs where not found among the list of available DMRs" << endl << "List of missing IOVs:" << endl;
232  for (int lostrun : lostruns) cout << to_string(lostrun) << " ";
233  cout << endl;
234  }
235  if(!missingruns.empty()){
236  cout << "WARNING: some runs are missing in the run/luminosity txt file" << endl << "List of missing runs:" << endl;
237  for (int missingrun : missingruns) cout << to_string(missingrun) << " ";
238  cout << endl;
239  }
240  }
241  return problemfound;
242 
243 }
244 
249 void DMRtrends(vector<int> IOVlist,vector<string> labels, TString Year, string myValidation, vector<string> geometries, vector<Color_t> colours, TString outputdir, bool pixelupdate, vector<int> pixelupdateruns, bool showlumi, bool FORCE){
250  compileDMRTrends(IOVlist, labels, Year, myValidation, geometries, showlumi, FORCE);
251  cout<< "Begin plotting"<<endl;
252  PlotDMRTrends(labels, Year, myValidation, geometries, colours, outputdir, pixelupdate, pixelupdateruns, showlumi);
253 
254 };
255 
260 void compileDMRTrends(vector<int> IOVlist, vector<string> labels, TString Year, string myValidation, vector<string> geometries, bool showlumi, bool FORCE){
261  gROOT->SetBatch();
262  vector<int>RunNumbers;
263  vector<TString> filenames;
264  TRegexp regexp("[0-9][0-9][0-9][0-9][0-9][0-9]");
265  for (const auto & entry : fs::recursive_directory_iterator(myValidation)){
266  bool found_all_labels = true;
267  for(string label : labels){
268  if(entry.path().string().find(label)==std::string::npos) found_all_labels = false;
269  }
270  if ((entry.path().string().find("ExtendedOfflineValidation_Images/OfflineValidationSummary.root")!=std::string::npos)&&found_all_labels){
271  if(fs::is_empty(entry.path())) cout << "ERROR: Empty file " << entry.path() << endl;
272  else{
273  TString filename(entry.path().string());
274  filenames.push_back(filename);
275  TString runstring(filename(regexp));
276  if(runstring.IsFloat()){
277  int runN=runstring.Atoi();
278  RunNumbers.push_back(runN);
279  }
280  }
281  }
282  }
283  if(checkrunlist(RunNumbers,IOVlist,Year)&&showlumi&&!FORCE){
284  cout << "Please fix the run/luminosities file!" << endl;
285  exit(EXIT_FAILURE);
286  }
287 
288  vector<TString> structures { "BPIX", "BPIX_y", "FPIX", "FPIX_y", "TIB", "TID", "TOB", "TEC"};
289 
290  const map<TString,int> nlayers{ {"BPIX", 4}, {"FPIX", 3}, {"TIB", 4}, {"TID", 3}, {"TOB", 6}, {"TEC", 9} };
291 
292 
293  map<pair<pair<TString,int>,TString>,Geometry> mappoints; // pair = (structure, layer), geometry
294 
295  std::sort(filenames.begin(),filenames.end());//order the files in alphabetical order
296  for (TString filename: filenames){
297  int runN;
298  TString runstring(filename(regexp));
299  if(runstring.IsFloat()){
300  runN=runstring.Atoi();
301  }
302  else{
303  cout << "ERROR: run number not retrieved for file " << filename << endl;
304  continue;
305  }
306 
307  TFile * f = new TFile(filename, "READ");
308 
309  for (TString& structure: structures) {
310  TString structname = structure;
311  structname.ReplaceAll("_y", "");
312  size_t layersnumber=nlayers.at(structname);
313  for (size_t layer=0; layer<=layersnumber;layer++){
314  for (string geometry: geometries) {
315  TString name = getName(structure, layer, geometry);
316  TH1F *histo = dynamic_cast<TH1F*>(f->Get( name));
317  //Geometry *geom =nullptr;
318  Point * point = nullptr;
319  // Three possibilities:
320  // - All histograms are produced correctly
321  // - Only the non-split histograms are produced
322  // - No histogram is produced correctly
323  // FORCE means that the Point is not added to the points collection in the chosen geometry for that structure
324  // If FORCE is not enabled a default value for the Point is used (-9999) which will appear in the plots
325  if(!histo){
326  cout << "Run" << runN << " Histogram: " << name << " not found" << endl;
327  if(FORCE)continue;
328  point= new Point(runN);
329  }else if(structure!="TID"&&structure!="TEC"){
330  TH1F *histoplus = dynamic_cast<TH1F*>(f->Get((name+"_plus")));
331  TH1F *histominus = dynamic_cast<TH1F*>(f->Get((name+"_minus")));
332  if(!histoplus||!histominus){
333  cout << "Run" << runN << " Histogram: " << name << " plus or minus not found" << endl;
334  if(FORCE)continue;
335  point= new Point(runN, histo);
336  }else point= new Point(runN, histo, histoplus, histominus);
337 
338  }else point= new Point(runN, histo);
339  mappoints[make_pair(make_pair(structure,layer),geometry)].points.push_back(*point);
340  }
341  }
342  }
343  f->Close();
344  }
345  TString outname=myValidation+"DMRtrends";
346  for(TString label : labels){outname+="_"; outname+=label;}
347  outname+=".root";
348  TFile * fout = TFile::Open(outname, "RECREATE");
349  for (TString& structure: structures) {
350  TString structname = structure;
351  structname.ReplaceAll("_y", "");
352  size_t layersnumber=nlayers.at(structname);
353  for (size_t layer=0; layer<=layersnumber;layer++){
354  for (string geometry : geometries) {
355  TString name = getName(structure, layer, geometry);
356  Geometry geom = mappoints[make_pair(make_pair(structure,layer),geometry)];
357  using Trend = vector<float> (Geometry::*)() const;
360  vector<TString> variables {"mu", "sigma", "muplus", "sigmaplus", "muminus", "sigmaminus", "deltamu", "sigmadeltamu" };
361  vector<float> runs = geom.Run();
362  size_t n = runs.size();
363  vector<float> emptyvec;
364  for(size_t i=0; i < runs.size(); i++)emptyvec.push_back(0.);
365  for(size_t iVar=0; iVar < variables.size(); iVar++){
366  Trend trend=trends.at(iVar);
367  TGraphErrors *g = new TGraphErrors(n, runs.data(), (geom.*trend)().data(), emptyvec.data(), emptyvec.data());
368  g->SetTitle(geometry.c_str());
369  g->Write(name+"_"+variables.at(iVar));
370  }
371  vector<pair<Trend,Trend>> trendspair {make_pair(&Geometry::Mu, &Geometry::Sigma), make_pair(&Geometry::MuPlus, &Geometry::SigmaPlus),
372  make_pair(&Geometry::MuMinus, &Geometry::SigmaMinus), make_pair(&Geometry::DeltaMu, &Geometry::SigmaDeltaMu) };
373  vector<pair<TString,TString>> variablepairs {make_pair("mu", "sigma"), make_pair("muplus", "sigmaplus"), make_pair("muminus", "sigmaminus"), make_pair("deltamu", "sigmadeltamu") };
374  for(size_t iVar=0; iVar < variablepairs.size(); iVar++){
375  Trend meantrend=trendspair.at(iVar).first;
376  Trend sigmatrend=trendspair.at(iVar).second;
377  TGraphErrors *g = new TGraphErrors(n, runs.data(), (geom.*meantrend)().data(), emptyvec.data(), (geom.*sigmatrend)().data());
378  g->SetTitle(geometry.c_str());
379  TString graphname = name+"_"+variablepairs.at(iVar).first;
380  graphname+=variablepairs.at(iVar).second;
381  g->Write(graphname);
382  }
383  }
384  }
385  }
386  fout->Close();
387 
388 }
389 
393 void PixelUpdateLines(TCanvas *c, TString Year, bool showlumi, vector<int>pixelupdateruns){
394  vector<TPaveText*> labels;
395  double lastlumi=0.;
396  c->cd();
397  size_t index = 0;
398  for(int pixelupdaterun : pixelupdateruns){
399  double lumi=0.;
400  if(showlumi)lumi=getintegratedlumiuptorun(pixelupdaterun,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.
401  else lumi=pixelupdaterun;
402  TLine *line = new TLine (lumi,c->GetUymin(),lumi,c->GetUymax());
403  line->SetLineColor(kBlue);
404  line->SetLineStyle(9);
405  line->Draw();
406  //Due to the way the coordinates within the Canvas are set, the following steps are required to draw the TPaveText:
407  // Compute the gPad coordinates in TRUE normalized space (NDC)
408  int ix1;
409  int ix2;
410  int iw = gPad->GetWw();
411  int ih = gPad->GetWh();
412  double x1p,y1p,x2p,y2p;
413  gPad->GetPadPar(x1p,y1p,x2p,y2p);
414  ix1 = (Int_t)(iw*x1p);
415  ix2 = (Int_t)(iw*x2p);
416  double wndc = TMath::Min(1.,(double)iw/(double)ih);
417  double rw = wndc/(double)iw;
418  double x1ndc = (double)ix1*rw;
419  double x2ndc = (double)ix2*rw;
420  // Ratios to convert user space in TRUE normalized space (NDC)
421  double rx1,ry1,rx2,ry2;
422  gPad->GetRange(rx1,ry1,rx2,ry2);
423  double rx = (x2ndc-x1ndc)/(rx2-rx1);
424  double _sx;
425  // Left limit of the TPaveText
426  _sx = rx*(lumi-rx1)+x1ndc;
427  // To avoid an overlap between the TPaveText a vertical shift is done when the IOVs are too close
428  if( _sx < lastlumi ){
429  index++;
430  }else index=0;
431  TPaveText *box= new TPaveText(_sx+0.0015,0.86-index*0.04,_sx+0.035,0.89-index*0.04,"blNDC");
432  box->SetFillColor(10);
433  box->SetBorderSize(1);
434  box->SetLineColor(kBlack);
435  TText *textRun = box->AddText(Form("%i",int(pixelupdaterun)));
436  textRun->SetTextSize(0.025);
437  labels.push_back(box);
438  lastlumi=_sx+0.035;
439  }
440  //Drawing in a separate loop to ensure that the labels are drawn on top of the lines
441  for(auto label: labels){
442  label->Draw("same");
443  }
444  c->Update();
445 }
446 
451 double getintegratedlumiuptorun(int run, TString Year, double min){
452  TGraph * scale = new TGraph((lumifileperyear(Year,"run")).Data());
453  int Nscale=scale->GetN();
454  double *xscale=scale->GetX();
455  double *yscale=scale->GetY();
456 
457 
458  double lumi=min;
459  int index=-1;
460  for(int j=0;j<Nscale;j++){
461  lumi+=yscale[j];
462  if(run>=xscale[j]){
463  index=j;
464  continue;
465  }
466  }
467  lumi=min;
468  for(int j=0;j<index;j++)lumi+=yscale[j]/lumiFactor;
469 
470  return lumi;
471 
472 }
477 void scalebylumi(TGraphErrors *g, TString Year, double min){
478  size_t N=g->GetN();
479  vector<double> x,y,xerr,yerr;
480 
481  TGraph * scale = new TGraph((lumifileperyear(Year,"IOV")).Data());
482  size_t Nscale=scale->GetN();
483  double *xscale=scale->GetX();
484  double *yscale=scale->GetY();
485 
486  size_t i=0;
487  while(i<N){
488  double run,yvalue;
489  g->GetPoint(i,run,yvalue);
490  size_t index=-1;
491  for(size_t j=0;j<Nscale;j++){
492  if(run==xscale[j]){
493  index=j;
494  continue;
495  }else if(run>xscale[j]) continue;
496  }
497  if(yscale[index]==0||index<0.){
498  N=N-1;
499  g->RemovePoint(i);
500  }else{
501  double xvalue=min;
502  for(size_t j=0;j<index;j++)xvalue+=yscale[j]/lumiFactor;
503  x .push_back(xvalue+(yscale[index]/(lumiFactor*2.)));
504  if(yvalue<=DUMMY){
505  y.push_back(DUMMY);
506  yerr.push_back(0.);
507  }else{
508  y.push_back(yvalue);
509  yerr.push_back(g->GetErrorY(i));
510  }
511  xerr.push_back(yscale[index]/(lumiFactor*2.));
512  i=i+1;
513  }
514 
515  }
516  g->GetHistogram()->Delete();
517  g->SetHistogram(nullptr);
518  for(size_t i=0;i<N;i++){ g->SetPoint(i, x.at(i),y.at(i)); g->SetPointError(i, xerr.at(i),yerr.at(i));}
519 
520 }
521 
526 TH1F *ConvertToHist(TGraphErrors *g){
527  size_t N=g->GetN();
528  double* x=g->GetX();
529  double* y=g->GetY();
530  double* xerr=g->GetEX();
531  vector<float> bins;
532  bins.push_back(x[0]-xerr[0]);
533  for(size_t i=1;i<N;i++){
534  if((x[i-1]+xerr[i-1]) > (bins.back() + pow(10,-6))) bins.push_back(x[i-1]+xerr[i-1]);
535  if((x[i]-xerr[i]) > (bins.back() + pow(10,-6))) bins.push_back(x[i]-xerr[i]);
536 
537  }
538  bins.push_back(x[N-1]+xerr[N-1]);
539  TString histoname="histo_";
540  histoname+=g->GetName();
541  TH1F *histo = new TH1F(histoname,g->GetTitle(),bins.size()-1,bins.data());
542  for(size_t i=0;i<N;i++){
543  histo->Fill(x[i],y[i]);
544  histo->SetBinError(histo->FindBin(x[i]),pow(10,-6));
545  }
546  return histo;
547 }
548 
549 
554 void PlotDMRTrends(vector<string> labels, TString Year, string myValidation, vector<string> geometries, vector<Color_t> colours, TString outputdir, bool pixelupdate, vector<int> pixelupdateruns, bool showlumi){
556  checkrunlist(pixelupdateruns,{},Year);
557  vector<TString> structures { "BPIX", "BPIX_y", "FPIX", "FPIX_y", "TIB", "TID", "TOB", "TEC"};
558 
559  const map<TString,int> nlayers{ {"BPIX", 4}, {"FPIX", 3}, {"TIB", 4}, {"TID", 3}, {"TOB", 6}, {"TEC", 9} };
560 
561 
562  TString filename=myValidation+"DMRtrends";
563  for(TString label : labels){ filename+="_"; filename+=label;}
564  filename+=".root";
565  TFile *in= new TFile(filename);
566  for (TString& structure: structures) {
567  TString structname = structure;
568  structname.ReplaceAll("_y", "");
569  int layersnumber=nlayers.at(structname);
570  for (int layer=0; layer<=layersnumber;layer++){
571  vector<TString> variables {"mu", "sigma", "muplus", "sigmaplus", "muminus", "sigmaminus", "deltamu", "sigmadeltamu", "musigma", "muplussigmaplus", "muminussigmaminus", "deltamusigmadeltamu"};
572  vector<string> YaxisNames { "#mu [#mum]", "#sigma_{#mu} [#mum]", "#mu outward [#mum]", "#sigma_{#mu outward} [#mum]", "#mu inward [#mum]", "#sigma_{#mu inward} [#mum]", "#Delta#mu [#mum]", "#sigma_{#Delta#mu} [#mum]", "#mu [#mum]", "#mu outward [#mum]", "#mu inward [#mum]", "#Delta#mu [#mum]",};
573  //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
574  for(size_t i=8; i < variables.size(); i++){
575  TString variable= variables.at(i);
576  TCanvas * c = new TCanvas("dummy","",2000,800);
577  vector<Color_t>::iterator colour = colours.begin();
578 
579  TMultiGraph *mg = new TMultiGraph(structure,structure);
580  THStack *mh = new THStack(structure,structure);
581  size_t igeom=0;
582  for (string geometry: geometries) {
583  TString name = getName(structure, layer, geometry);
584  TGraphErrors *g = dynamic_cast<TGraphErrors*> (in->Get(name+"_"+variables.at(i)));
585  g->SetName(name+"_"+variables.at(i));
586  if(i>=8){
587  g->SetLineWidth(1);
588  g->SetLineColor(*colour);
589  g->SetFillColorAlpha(*colour,0.2);
590  }
591  vector<vector<double>> vectors;
592  //if(showlumi&&i<8)scalebylumi(dynamic_cast<TGraph*>(g));
593  if(showlumi)scalebylumi(g,Year);
594  g->SetLineColor(*colour);
595  g->SetMarkerColor(*colour);
596  TH1F *h = ConvertToHist(g);
597  h->SetLineColor(*colour);
598  h->SetMarkerColor(*colour);
599  h->SetMarkerSize(0);
600  h->SetLineWidth(1);
601 
602  if(i<8){
603  mg->Add(g,"PL");
604  mh->Add(h,"E");
605  }else{
606  mg->Add(g,"2");
607  mh->Add(h,"E");
608  }
609  ++colour;
610  ++igeom;
611  }
612  if(i<8){
613  mg->Draw("a");
614  }else{
615  mg->Draw("a2");
616  }
617  double max=6;
618  double min=-4;
619  double range=max-min;
620  if(((variable=="sigma"||variable=="sigmaplus"||variable=="sigmaminus"||variable=="sigmadeltamu")&&range>=2)){
621  mg->SetMaximum(4);
622  mg->SetMinimum(-2);
623  }else{
624  mg->SetMaximum(max+range*0.1);
625  mg->SetMinimum(min-range*0.3);
626  }
627 
628  char* Ytitle= (char *)YaxisNames.at(i).c_str();
629  mg->GetYaxis()->SetTitle(Ytitle);
630  mg->GetXaxis()->SetTitle(showlumi ? "Integrated lumi [1/fb]" : "IOV number");
631  mg->GetXaxis()->CenterTitle(true);
632  mg->GetYaxis()->CenterTitle(true);
633  mg->GetYaxis()->SetTitleOffset(.5);
634  mg->GetYaxis()->SetTitleSize(.05);
635  mg->GetXaxis()->SetTitleSize(.04);
636  if(showlumi) mg->GetXaxis()->SetLimits(0.,mg->GetXaxis()->GetXmax());
637  gStyle->SetOptTitle(0); // TODO
638  gPad->SetTickx();
639  gPad->SetTicky();
640  //c->SetLeftMargin(0.11);
641 
642  c->Update();
643 
644  //gStyle->SetLegendBorderSize(0);
645  gStyle->SetLegendTextSize(0.025);
646 
647  //TLegend *legend = c->BuildLegend();
648  TLegend *legend = c->BuildLegend(0.3,0.15,0.3,0.15);
649  int Ngeom=geometries.size();
650  legend->SetNColumns(Ngeom);
651  //legend->SetTextSize(0.05);
652  TString structtitle = structure;
653  if(layer!=0){
654  if(structure=="TID"||structure=="TEC"||structure=="FPIX"||structure=="FPIX_y")structtitle+="_disc";
655  else structtitle+="_layer";
656  structtitle+=layer;
657  }
658  legend->SetHeader("#scale[1.2]{#bf{CMS} Work in progress}");
659  //TLegendEntry *header = (TLegendEntry*)legend->GetListOfPrimitives()->First();
660  //header->SetTextSize(.04);
661  legend->AddEntry((TObject*)nullptr,structtitle.Data(),"h");
662 
663  //TLegendEntry *str = (TLegendEntry*)legend->GetListOfPrimitives()->Last();
664  //str->SetTextSize(.03);
665  PixelUpdateLines(c, Year, showlumi, pixelupdateruns);
666 
667  legend->Draw();
668  mh->Draw("nostack same");
669  c->Update();
670  TString structandlayer = getName(structure,layer,"");
671  TString printfile=outputdir;
672  for(TString label : labels){printfile+=label;printfile+="_";}
673  printfile+=variable+structandlayer;
674  c->SaveAs(printfile+".pdf");
675  c->SaveAs(printfile+".eps");
676  c->SaveAs(printfile+".png");
677  c->Destructor();
678  }
679 
680  }
681 
682  }
683  in->Close();
684 }
685 
701 int main (int argc, char * argv[]) {
702  if (argc == 1) {
703 
704  vector<int>IOVlist={290543, 296702, 296966, 297224, 297281, 297429, 297467, 297484, 297494, 297503, 297557, 297599, 297620, 297660, 297670, 298678, 298996, 299062, 299096, 299184, 299327, 299368, 299381, 299443, 299480, 299592, 299594, 299649, 300087, 300155, 300233, 300237, 300280, 300364, 300389, 300399, 300459, 300497, 300515, 300538, 300551, 300574, 300636, 300673, 300780, 300806, 300812, 301046, 301417, 302131, 302573, 302635, 303825, 303998, 304170, 304505, 304672, 305040, 305081};//UL17
705  //vector<int> pixelupdateruns {316758, 317527,317661,317664,318227, 320377};//2018
706  vector<int> pixelupdateruns {290543, 297281, 298653, 299443, 300389, 302131, 303790, 304911, 305898};//2017
707 
708  cout << "WARNING: Running function with arguments specified in DMRtrends.cc" << endl << "If you want to specify the arguments from command line run the macro as follows:" << endl << "DMRtrends labels pathtoDMRs geometriesandcolourspairs outputdirectory showpixelupdate showlumi FORCE" << endl;
709 
710  //Example provided for a currently working set of parameters
711  DMRtrends(IOVlist,{"vUL17","MB"},"2017", "/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/results/acardini/DMRs/EOY17/", {"EOY17","full ML pixel + strip","SG-mp2607","mp2993"}, {kRed, kBlack, kBlue, kGreen}, "/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/alignmentObjects/acardini/DMRsTrends/", true, pixelupdateruns, true, true);
712 
713 
714  return 0;
715  }
716  else if (argc < 11) {
717  cout << "DMRtrends IOVlist labels Year pathtoDMRs geometriesandcolourspairs outputdirectory pixelupdatelist showpixelupdate showlumi FORCE" << endl;
718 
719  return 1;
720  }
721 
722  TString runlist = argv[1],
723  all_labels = argv[2],
724  Year = argv[3],
725  pathtoDMRs = argv[4],
726  geometrieandcolours = argv[5], //name1:title1:color1,name2:title2:color2,name3:title3:color3
727  outputdirectory = argv[6],
728  pixelupdatelist = argv[7];
729  bool showpixelupdate = argv[8],
730  showlumi = argv[9],
731  FORCE = argv[10];
732  TObjArray *labelarray = all_labels.Tokenize(",");
733  vector<string> labels;
734  for(int i=0; i < labelarray->GetEntries(); i++)labels.push_back((string)(labelarray->At(i)->GetName()));
735  TObjArray *IOVarray = runlist.Tokenize(",");
736  vector<int> IOVlist;
737  for(int i=0; i < IOVarray->GetEntries(); i++)IOVlist.push_back(stoi(IOVarray->At(i)->GetName()));
738  vector<int> pixelupdateruns;
739  TObjArray *PIXarray = pixelupdatelist.Tokenize(",");
740  for(int i=0; i < PIXarray->GetEntries(); i++)pixelupdateruns.push_back(stoi(PIXarray->At(i)->GetName()));
741  vector<string> geometries;
742  vector<Color_t> colours;
743  TObjArray *geometrieandcolourspairs = geometrieandcolours.Tokenize(",");
744  for (int i=0; i < geometrieandcolourspairs->GetEntries(); i++) {
745  TObjArray *geomandcolourvec = TString(geometrieandcolourspairs->At(i)->GetName()).Tokenize(":");
746  geometries.push_back(geomandcolourvec->At(0)->GetName());
747  colours.push_back((Color_t)(atoi(geomandcolourvec->At(1)->GetName())));
748  }
749  DMRtrends(IOVlist,labels,Year,pathtoDMRs.Data(),geometries,colours,outputdirectory.Data(),showpixelupdate, pixelupdateruns,showlumi,FORCE);
750 
751 
752  return 0;
753 }
Geometry(TString Title)
Definition: DMRtrends.cc:135
vector< float > SigmaPlus() const
Definition: DMRtrends.cc:144
float mu
Definition: DMRtrends.cc:56
Geometry & operator=(const Geometry &geom)
Definition: DMRtrends.cc:136
#define DUMMY
Definition: DMRtrends.cc:33
vector< float > SigmaDeltaMu() const
Definition: DMRtrends.cc:147
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
TString lumifileperyear(TString Year="2018", string RunOrIOV="IOV")
Definition: DMRtrends.cc:183
vector< Point > points
Definition: DMRtrends.cc:121
void DMRtrends(vector< int > IOVlist, vector< string >labels={"MB"}, TString Year="2018", string myValidation="/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/results/acardini/DMRs/", vector< string > geometries={"GT","SG","MP pix LBL","PIX HLS+ML STR fix"}, vector< Color_t > colours={kBlue, kRed, kGreen, kCyan}, TString outputdir="/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/alignmentObjects/acardini/DMRsTrends/", bool pixelupdate=false, vector< int > pixelupdateruns={314881, 316758, 317527, 318228, 320377}, bool showlumi=false, bool FORCE=false)
Definition: DMRtrends.cc:249
float sigma
Definition: DMRtrends.cc:56
TString title
Definition: DMRtrends.cc:133
Point(float Run, TH1 *histo, TH1 *histoplus, TH1 *histominus)
Definition: DMRtrends.cc:68
double getintegratedlumiuptorun(int run, TString Year="2018", double min=0.)
Definition: DMRtrends.cc:451
T Min(T a, T b)
Definition: MathUtil.h:39
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:20
void SetTitle(TString Title)
Definition: DMRtrends.cc:137
std::pair< double, double > Point
Definition: CaloEllipse.h:18
vector< float > DeltaMu() const
Definition: DMRtrends.cc:146
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
bool checkrunlist(vector< int > runs, vector< int > IOVlist={}, TString Year="2018")
Definition: DMRtrends.cc:211
Point(float Run, TH1 *histo)
Definition: DMRtrends.cc:77
void PlotDMRTrends(vector< string >labels={"MB"}, TString Year="2018", string myValidation="/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/results/acardini/DMRs/", vector< string > geometries={"GT","SG","MP pix LBL","PIX HLS+ML STR fix"}, vector< Color_t > colours={kBlue, kRed, kGreen, kCyan}, TString outputdir="/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/alignmentObjects/acardini/DMRsTrends/", bool pixelupdate=false, vector< int > pixelupdateruns={314881, 316758, 317527, 318228, 320377}, bool showlumi=false)
Definition: DMRtrends.cc:554
vector< float > MuMinus() const
Definition: DMRtrends.cc:142
Point & operator=(const Point &p)
Definition: DMRtrends.cc:81
void scalebylumi(TGraphErrors *g, TString Year="2018", double min=0.)
Definition: DMRtrends.cc:477
char const * label
vector< int > runlistfromlumifile(TString Year="2018")
Definition: DMRtrends.cc:198
Class Geometry Contains vector for fit parameters (mean, sigma, etc.) obtained from multiple IOVs See...
Definition: DMRtrends.cc:119
float GetRun() const
Definition: DMRtrends.cc:83
vector< float > Run() const
Definition: DMRtrends.cc:139
vector< float > Sigma() const
Definition: DMRtrends.cc:143
float run
Definition: DMRtrends.cc:56
float muminus
Definition: DMRtrends.cc:56
void compileDMRTrends(vector< int > IOVlist, vector< string >labels={"MB"}, TString Year="2018", string myValidation="/afs/cern.ch/cms/CAF/CMSALCA/ALCA_TRACKERALIGN/data/commonValidation/results/acardini/DMRs/", vector< string > geometries={"GT","SG","MP pix LBL","PIX HLS+ML STR fix"}, bool showlumi=false, bool FORCE=false)
Definition: DMRtrends.cc:260
double f[11][100]
float GetSigmaMinus() const
Definition: DMRtrends.cc:89
const int mu
Definition: Constants.h:22
float GetMuMinus() const
Definition: DMRtrends.cc:86
Definition: value.py:1
T min(T a, T b)
Definition: MathUtil.h:58
void PixelUpdateLines(TCanvas *c, TString Year="2018", bool showlumi=false, vector< int >pixelupdateruns={314881, 316758, 317527, 318228, 320377})
Definition: DMRtrends.cc:393
vector< float > Mu() const
Definition: DMRtrends.cc:140
float sigmaplus
Definition: DMRtrends.cc:56
TString getName(TString structure, int layer, TString geometry)
Definition: DMRtrends.cc:167
#define lumiFactor
Definition: DMRtrends.cc:37
#define N
Definition: blowfish.cc:9
def is_empty(h)
Definition: utils.py:177
float GetSigmaDeltaMu() const
Definition: DMRtrends.cc:92
gErrorIgnoreLevel
Definition: utils.py:28
float GetMu() const
Definition: DMRtrends.cc:84
float GetSigmaPlus() const
Definition: DMRtrends.cc:88
vector< float > SigmaMinus() const
Definition: DMRtrends.cc:145
Structure Point Contains parameters of Gaussian fits to DMRs.
Definition: DMRtrends.cc:55
vector< float > MuPlus() const
Definition: DMRtrends.cc:141
float GetMuPlus() const
Definition: DMRtrends.cc:85
float sigmaminus
Definition: DMRtrends.cc:56
vector< float > GetQuantity(float(Point::*getter)() const) const
Definition: DMRtrends.cc:124
float GetSigma() const
Definition: DMRtrends.cc:87
#define DMRFactor
Definition: DMRtrends.cc:41
Point(float Run=-999., float y1=-999., float y2=-999., float y3=-999., float y4=-999., float y5=-999., float y6=-999.)
Definition: DMRtrends.cc:61
TString GetTitle()
Definition: DMRtrends.cc:138
float GetDeltaMu() const
Definition: DMRtrends.cc:90
float muplus
Definition: DMRtrends.cc:56
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
int main(int argc, char *argv[])
Definition: DMRtrends.cc:701
TH1F * ConvertToHist(TGraphErrors *g)
Definition: DMRtrends.cc:526