LOGO

RestFrames  v1.0.1
RestFrames HEP Event Analysis Software Library
HistPlot.cc
Go to the documentation of this file.
1 // RestFrames: particle physics event analysis library
3 // --------------------------------------------------------------------
4 // Copyright (c) 2014-2016, Christopher Rogan
14 // This file is part of RestFrames.
15 //
16 // RestFrames is free software; you can redistribute it and/or modify
17 // it under the terms of the GNU General Public License as published by
18 // the Free Software Foundation; either version 2 of the License, or
19 // (at your option) any later version.
20 //
21 // RestFrames is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with RestFrames. If not, see <http://www.gnu.org/licenses/>.
29 
30 #include <TFile.h>
31 #include <TLatex.h>
32 #include <TLegend.h>
33 
34 #include "RestFrames/HistPlot.hh"
37 
38 namespace RestFrames {
39 
40  HistPlot::HistPlot(const std::string& sname, const std::string& stitle)
41  : RFPlot(sname, stitle)
42  {
43  SetPlotLabel("#bf{#it{RestFrames}} Event Generation");
44  SetPlotTitle(GetTitle());
45  SetScaleLabel("a. u.");
46  m_Scale = 1.;
47  m_SetScale = false;
48  m_Rebin = 4;
49  m_Log.SetSource("HistPlot "+GetName());
50  SetStyle();
51  }
52 
53  HistPlot::~HistPlot(){
54  int Nv = m_Vars.GetN();
55  for(int i = 0; i < Nv; i++)
56  delete &m_Vars[i];
57  int Nc = m_Cats.GetN();
58  for(int i = 0; i < Nc; i++)
59  delete &m_Cats[i];
60  Clear();
61  }
62 
64  int N = m_1DHists.size();
65  for(int i = 0; i < N; i++)
66  delete m_1DHists[i];
67  m_1DHists.clear();
68  N = m_2DHists.size();
69  for(int i = 0; i < N; i++)
70  delete m_2DHists[i];
71  m_2DHists.clear();
72 
73  m_HistToVar.clear();
74  m_HistToVars.clear();
75  m_CatToHist1D.clear();
76  m_CatToHist2D.clear();
77  m_Plot1D_Var.clear();
78  m_Plot1D_Cats.clear();
79  m_Plot1D_Color.clear();
80  m_Plot2D_Vars.clear();
81  m_Plot2D_Cat.clear();
82  m_Plot2D_Color.clear();
83  RFPlot::Clear();
84  }
85 
86  HistPlotVar const& HistPlot::GetNewVar(const std::string& name,
87  const std::string& title,
88  double minval, double maxval,
89  const std::string& unit){
90  HistPlotVar* var = new HistPlotVar(name,title,minval,maxval,unit);
91  m_Vars += *var;
92  return *var;
93  }
94 
95  HistPlotCategory const& HistPlot::GetNewCategory(const std::string& name,
96  const std::string& title){
97  HistPlotCategory* cat = new HistPlotCategory(name, title);
98  m_Cats += *cat;
99  return *cat;
100  }
101 
102  void HistPlot::AddPlot(const HistPlotVar& var,
103  RFList<const HistPlotCategory> cats,
104  bool invert_colors){
105  if(!m_Vars.Contains(var))
106  return;
107 
108  int Ncat = cats.GetN();
109  for(int i = 0; i < Ncat; i++){
110  if(m_Cats.Contains(cats[i]))
111  continue;
112  Ncat--;
113  cats -= cats[i];
114  }
115 
116  if(Ncat == 0){
117  const HistPlotCategory* empty = &HistPlotCategory::Empty();
118  if(m_CatToHist1D.count(empty) <= 0)
119  m_CatToHist1D[empty] = std::vector<TH1D*>();
120 
121  int Nhist = m_CatToHist1D[empty].size();
122  bool exists = false;
123  for(int i = 0; i < Nhist; i++){
124  if(*m_HistToVar[m_CatToHist1D[empty][i]] == var){
125  exists = true;
126  break;
127  }
128  }
129  if(!exists){
130  std::string name = GetUniqueName(var.GetName()+"_"+GetName());
131  TH1D* hist = new TH1D(name.c_str(),name.c_str(),
132  256,var.GetMin(),var.GetMax());
133  m_HistToVar[hist] = &var;
134  m_CatToHist1D[empty].push_back(hist);
135  m_1DHists.push_back(hist);
136  }
137  } else {
138  for(int c = 0; c < Ncat; c++){
139  if(m_CatToHist1D.count(&cats[c]) <= 0)
140  m_CatToHist1D[&cats[c]] = std::vector<TH1D*>();
141 
142  int Nhist = m_CatToHist1D[&cats[c]].size();
143  bool exists = false;
144  for(int i = 0; i < Nhist; i++){
145  if(*m_HistToVar[m_CatToHist1D[&cats[c]][i]] == var){
146  exists = true;
147  break;
148  }
149  }
150  if(!exists){
151  std::string name = GetUniqueName(var.GetName()+"_"+cats[c].GetName()+"_"+GetName());
152  TH1D* hist = new TH1D(name.c_str(),name.c_str(),
153  256,var.GetMin(),var.GetMax());
154  m_HistToVar[hist] = &var;
155  m_CatToHist1D[&cats[c]].push_back(hist);
156  m_1DHists.push_back(hist);
157  }
158  }
159  }
160  m_Plot1D_Var.push_back(&var);
161  m_Plot1D_Cats.push_back(cats);
162  m_Plot1D_Color.push_back(invert_colors);
163  }
164 
165  void HistPlot::AddPlot(const HistPlotVar& varX, const HistPlotVar& varY,
166  RFList<const HistPlotCategory> cats,
167  bool invert_colors){
168  if(!m_Vars.Contains(varX) ||
169  !m_Vars.Contains(varY))
170  return;
171 
172  int Ncat = cats.GetN();
173  for(int i = 0; i < Ncat; i++){
174  if(m_Cats.Contains(cats[i]))
175  continue;
176  Ncat--;
177  cats -= cats[i];
178  }
179 
180  if(Ncat == 0){
181  const HistPlotCategory* empty = &HistPlotCategory::Empty();
182  if(m_CatToHist2D.count(empty) <= 0)
183  m_CatToHist2D[empty] = std::vector<TH2D*>();
184 
185  int Nhist = m_CatToHist2D[empty].size();
186  bool exists = false;
187  for(int i = 0; i < Nhist; i++){
188  if(*m_HistToVars[m_CatToHist2D[empty][i]].first == varX &&
189  *m_HistToVars[m_CatToHist2D[empty][i]].second == varY){
190  exists = true;
191  break;
192  }
193  }
194  if(!exists){
195  std::string name = GetUniqueName(varX.GetName()+"_v_"+varY.GetName()+"_"+GetName());
196  TH2D* hist = new TH2D(name.c_str(),name.c_str(),
197  128,varX.GetMin(),varX.GetMax(),
198  128,varY.GetMin(),varY.GetMax());
199 
200  m_HistToVars[hist] =
201  std::pair<const HistPlotVar*,const HistPlotVar*>(&varX,&varY);
202  m_CatToHist2D[empty].push_back(hist);
203  m_2DHists.push_back(hist);
204  m_Plot2D_Vars.push_back(std::pair<const HistPlotVar*,
205  const HistPlotVar*>(&varX,&varY));
206  m_Plot2D_Cat.push_back(empty);
207  m_Plot2D_Color.push_back(invert_colors);
208  }
209  } else {
210  for(int c = 0; c < Ncat; c++){
211  if(m_CatToHist2D.count(&cats[c]) <= 0)
212  m_CatToHist2D[&cats[c]] = std::vector<TH2D*>();
213 
214  int Nhist = m_CatToHist2D[&cats[c]].size();
215  bool exists = false;
216  for(int i = 0; i < Nhist; i++){
217  if(*m_HistToVars[m_CatToHist2D[&cats[c]][i]].first == varX &&
218  *m_HistToVars[m_CatToHist2D[&cats[c]][i]].second == varY){
219  exists = true;
220  break;
221  }
222  }
223  if(!exists){
224  std::string name = GetUniqueName(varX.GetName()+"_v_"+varY.GetName()+"_"+
225  cats[c].GetName()+"_"+GetName());
226  TH2D* hist = new TH2D(name.c_str(),name.c_str(),
227  128,varX.GetMin(),varX.GetMax(),
228  128,varY.GetMin(),varY.GetMax());
229  m_HistToVars[hist] =
230  std::pair<const HistPlotVar*,const HistPlotVar*>(&varX,&varY);
231  m_CatToHist2D[&cats[c]].push_back(hist);
232  m_2DHists.push_back(hist);
233  m_Plot2D_Vars.push_back(std::pair<const HistPlotVar*,
234  const HistPlotVar*>(&varX,&varY));
235  m_Plot2D_Cat.push_back(&cats[c]);
236  m_Plot2D_Color.push_back(invert_colors);
237  }
238  }
239  }
240  }
241 
242  void HistPlot::Fill(double weight){
243  const HistPlotCategory* empty = &HistPlotCategory::Empty();
244  int N = m_CatToHist1D[empty].size();
245  for(int i = 0; i < N; i++)
246  m_CatToHist1D[empty][i]->Fill(m_HistToVar[m_CatToHist1D[empty][i]]->GetVal(),
247  weight);
248  N = m_CatToHist2D[empty].size();
249  for(int i = 0; i < N; i++)
250  m_CatToHist2D[empty][i]->Fill(m_HistToVars[m_CatToHist2D[empty][i]].first->GetVal(),
251  m_HistToVars[m_CatToHist2D[empty][i]].second->GetVal(),
252  weight);
253  }
254 
255  void HistPlot::Fill(const HistPlotCategory& cat, double weight){
256  if(!cat){
257  Fill(weight);
258  return;
259  }
260 
261  int N = m_CatToHist1D[&cat].size();
262  for(int i = 0; i < N; i++)
263  m_CatToHist1D[&cat][i]->Fill(m_HistToVar[m_CatToHist1D[&cat][i]]->GetVal(),
264  weight);
265  N = m_CatToHist2D[&cat].size();
266  for(int i = 0; i < N; i++)
267  m_CatToHist2D[&cat][i]->Fill(m_HistToVars[m_CatToHist2D[&cat][i]].first->GetVal(),
268  m_HistToVars[m_CatToHist2D[&cat][i]].second->GetVal(),
269  weight);
270  }
271 
272  void HistPlot::Draw(bool invert_colors){
273  int N1D = m_Plot1D_Var.size();
274  int N2D = m_Plot2D_Vars.size();
275 
276  for(int i = 0; i < N2D; i++)
277  DrawPlot(m_Plot2D_Vars[i], *m_Plot2D_Cat[i],
278  (m_Plot2D_Color[i] || invert_colors));
279  for(int i = 0; i < N1D; i++)
280  DrawPlot(*m_Plot1D_Var[i], m_Plot1D_Cats[i],
281  (m_Plot1D_Color[i] || invert_colors));
282  }
283 
284  void HistPlot::DrawPlot(const HistPlotVar& var,
285  const HistPlotCatList& cats,
286  bool invert_colors){
287  std::vector<TH1D*> hists;
288  int Ncat = cats.GetN();
289  std::string catname = "";
290 
291  if(Ncat == 0){
292  const HistPlotCategory* empty = &HistPlotCategory::Empty();
293  int Nhist = m_CatToHist1D[empty].size();
294  for(int i = 0; i < Nhist; i++){
295  if(*m_HistToVar[m_CatToHist1D[empty][i]] == var){
296  hists.push_back(m_CatToHist1D[empty][i]);
297  break;
298  }
299  }
300  } else {
301  for(int c = 0; c < Ncat; c++){
302  catname += cats[c].GetName() + "_";
303  int Nhist = m_CatToHist1D[&cats[c]].size();
304  for(int i = 0; i < Nhist; i++){
305  if(*m_HistToVar[m_CatToHist1D[&cats[c]][i]] == var){
306  hists.push_back(m_CatToHist1D[&cats[c]][i]);
307  break;
308  }
309  }
310  }
311  }
312 
313  std::string name = GetUniqueName("c_"+var.GetName()+"_"+catname+GetName());
314  TCanvas* can = new TCanvas(name.c_str(),name.c_str(),600,500);
315  can->SetLeftMargin(0.2);
316  can->SetRightMargin(0.05);
317  can->Draw();
318  if(invert_colors){
319  can->SetFillColor(kBlack);
320  can->Modified();
321  }
322  can->SetGridx();
323  can->SetGridy();
324 
325  std::string XLabel = var.GetTitle();
326 
327  std::string ScaleLabel;
328  if(!m_SetScale){
329  ScaleLabel = "#frac{1}{N} #frac{dN}{";
330  ScaleLabel += "d( "+XLabel+" )}";
331  } else {
332  ScaleLabel = m_ScaleLabel;
333  }
334 
335  if(var.GetUnit() != "")
336  XLabel += " "+var.GetUnit();
337 
338  int N = hists.size();
339 
340  int imax = -1;
341  int imin = -1;
342  int imin0 = -1;
343  double hmax = -1.;
344  double hmin = -1.;
345  double hmin0 = -1.;
346  for(int i = 0; i < N; i++){
347  hists[i]->Rebin(m_Rebin);
348  if(!m_SetScale){
349  if(hists[i]->Integral() > 0.)
350  hists[i]->Scale(1./hists[i]->Integral());
351  } else {
352  hists[i]->Scale(m_Scale);
353  }
354  if(hists[i]->GetMaximum() > hmax || imax < 0){
355  hmax = hists[i]->GetMaximum();
356  imax = i;
357  }
358  if(hists[i]->GetMinimum(0.) < hmin || imin < 0){
359  hmin = hists[i]->GetMinimum(0.);
360  imin = i;
361  }
362  if(hists[i]->GetMinimum() < hmin0 || imin0 < 0){
363  hmin0 = hists[i]->GetMinimum();
364  imin0 = i;
365  }
366  }
367 
368  hists[imax]->Draw();
369  hists[imax]->GetXaxis()->SetTitle(XLabel.c_str());
370  hists[imax]->GetXaxis()->SetTitleOffset(1.27);
371  hists[imax]->GetXaxis()->CenterTitle();
372  hists[imax]->GetYaxis()->SetTitle(ScaleLabel.c_str());
373  hists[imax]->GetYaxis()->SetTitleOffset(1.42);
374  hists[imax]->GetYaxis()->CenterTitle();
375  if(hmin0 > 0.)
376  hists[imax]->GetYaxis()->SetRangeUser(0., 1.1*hmax);
377  else
378  hists[imax]->GetYaxis()->SetRangeUser(0.9*hmin, 1.1*hmax);
379 
380  if(invert_colors){
381  hists[imax]->GetXaxis()->SetTitleColor(kWhite);
382  hists[imax]->GetXaxis()->SetLabelColor(kWhite);
383  hists[imax]->GetXaxis()->SetAxisColor(kWhite);
384  hists[imax]->GetYaxis()->SetTitleColor(kWhite);
385  hists[imax]->GetYaxis()->SetLabelColor(kWhite);
386  hists[imax]->GetYaxis()->SetAxisColor(kWhite);
387  }
388 
389  for(int i = N-1; i >= 0; i--){
390  int icolor0 = 7003 + (i%8)*10;
391  int icolor1 = 7000 + (i%8)*10;
392  if(invert_colors){
393  icolor0 = 7000 + (i%8)*10;
394  icolor1 = 7003 + (i%8)*10;
395  }
396  hists[i]->SetFillColor(icolor1);
397  hists[i]->SetFillStyle(3002);
398  hists[i]->SetLineColor(icolor0);
399  hists[i]->SetLineWidth(3);
400  hists[i]->SetMarkerColor(icolor0);
401  hists[i]->SetMarkerSize(0);
402  hists[i]->Draw("same");
403  }
404 
405  TLatex l(0.6,0.943,m_PlotTitle.c_str());
406  l.SetNDC();
407  if(invert_colors)
408  l.SetTextColor(kWhite);
409  l.SetTextSize(0.045);
410  l.SetTextFont(132);
411  l.DrawLatex(0.48+std::max(0.,0.47-l.GetXsize()),0.947,m_PlotTitle.c_str());
412  l.SetTextSize(0.04);
413  l.SetTextFont(42);
414  l.DrawLatex(0.02,0.95,m_PlotLabel.c_str());
415 
416  if(N > 1){
417  TLegend* leg = new TLegend(0.225,std::max(0.5,0.884-double(N)*0.073),0.488,0.884);
418  AddTObject(leg);
419  leg->SetShadowColor(kWhite);
420  leg->SetLineColor(kWhite);
421  leg->SetFillColor(kWhite);
422  if(invert_colors){
423  leg->SetShadowColor(kBlack);
424  leg->SetLineColor(kBlack);
425  leg->SetFillColor(kBlack);
426  leg->SetTextColor(kWhite);
427  }
428  leg->SetTextFont(132);
429  leg->SetTextSize(0.045);
430  for(int i = 0; i < N; i++)
431  leg->AddEntry(hists[i],cats[i].GetTitle().c_str());
432 
433  leg->Draw();
434  }
435  if(N == 1 && !cats[0].IsEmpty()){
436  l.SetTextSize(0.05);
437  l.SetTextFont(132);
438  l.DrawLatex(0.23,0.84,cats[0].GetTitle().c_str());
439  }
440 
441  AddCanvas(can);
442  }
443 
444  void HistPlot::DrawPlot(const std::pair<const HistPlotVar*,
445  const HistPlotVar*>& vars,
446  const HistPlotCategory& cat,
447  bool invert_colors){
448  const HistPlotVar& varX = *vars.first;
449  const HistPlotVar& varY = *vars.second;
450  TH2D* hist = nullptr;
451  std::string catname = "";
452  if(!cat){
453  const HistPlotCategory* empty = &HistPlotCategory::Empty();
454  int Nhist = m_CatToHist2D[empty].size();
455  for(int i = 0; i < Nhist; i++){
456  if(*m_HistToVars[m_CatToHist2D[empty][i]].first == varX &&
457  *m_HistToVars[m_CatToHist2D[empty][i]].second == varY){
458  hist = m_CatToHist2D[empty][i];
459  break;
460  }
461  }
462  } else {
463  catname += cat.GetName() + "_";
464  int Nhist = m_CatToHist2D[&cat].size();
465  for(int i = 0; i < Nhist; i++){
466  if(*m_HistToVars[m_CatToHist2D[&cat][i]].first == varX &&
467  *m_HistToVars[m_CatToHist2D[&cat][i]].second == varY){
468  hist = m_CatToHist2D[&cat][i];
469  break;
470  }
471  }
472  }
473 
474  std::string name = GetUniqueName("c_"+varX.GetName()+"_v_"+varY.GetName()+"_"+catname+GetName());
475  TCanvas* can = new TCanvas(name.c_str(),name.c_str(),600,500);
476  can->Draw();
477  if(invert_colors){
478  can->SetFillColor(kBlack);
479  can->Modified();
480  }
481  can->SetGridx();
482  can->SetGridy();
483  can->SetLogz();
484 
485  std::string XLabel = varX.GetTitle();
486  std::string YLabel = varY.GetTitle();
487 
488  std::string ScaleLabel;
489  if(!m_SetScale){
490  if(hist->Integral() > 0.)
491  hist->Scale(1./hist->Integral());
492  ScaleLabel = "#frac{1}{N} #frac{dN}{";
493  ScaleLabel += "d( "+XLabel+" ) ";
494  ScaleLabel += "d( "+YLabel+" )";
495  ScaleLabel += "}";
496  } else {
497  hist->Scale(m_Scale);
498  ScaleLabel = m_ScaleLabel;
499  }
500 
501  hist->RebinX(m_Rebin);
502  hist->RebinY(m_Rebin);
503 
504  if(varX.GetUnit() != "")
505  XLabel += " "+varX.GetUnit();
506  if(varY.GetUnit() != "")
507  YLabel += " "+varY.GetUnit();
508 
509  SetZPalette(invert_colors);
510 
511  hist->Draw("COLZ");
512  hist->GetXaxis()->SetTitle(XLabel.c_str());
513  hist->GetXaxis()->SetTitleOffset(1.24);
514  hist->GetXaxis()->CenterTitle();
515  hist->GetYaxis()->SetTitle(YLabel.c_str());
516  hist->GetYaxis()->SetTitleOffset(1.11);
517  hist->GetYaxis()->CenterTitle();
518  hist->GetZaxis()->SetTitle(ScaleLabel.c_str());
519  hist->GetZaxis()->SetTitleOffset(1.5);
520  hist->GetZaxis()->CenterTitle();
521  hist->GetZaxis()->SetRangeUser(0.9*hist->GetMinimum(0.0),1.1*hist->GetMaximum());
522  if(invert_colors){
523  hist->GetXaxis()->SetTitleColor(kWhite);
524  hist->GetXaxis()->SetLabelColor(kWhite);
525  hist->GetXaxis()->SetAxisColor(kWhite);
526  hist->GetYaxis()->SetTitleColor(kWhite);
527  hist->GetYaxis()->SetLabelColor(kWhite);
528  hist->GetYaxis()->SetAxisColor(kWhite);
529  hist->GetZaxis()->SetTitleColor(kWhite);
530  hist->GetZaxis()->SetLabelColor(kWhite);
531  hist->GetZaxis()->SetAxisColor(kWhite);
532  }
533  hist->Draw("COLZ");
534 
535  std::string title = m_PlotTitle;
536  if(!cat.IsEmpty())
537  title = cat.GetTitle();
538 
539  TLatex l(0.6,0.943,title.c_str());
540  l.SetNDC();
541  if(invert_colors)
542  l.SetTextColor(kWhite);
543  l.SetTextSize(0.045);
544  l.SetTextFont(132);
545  l.DrawLatex(0.48+std::max(0.,0.32-l.GetXsize()),0.947,title.c_str());
546  l.SetTextSize(0.04);
547  l.SetTextFont(42);
548  l.DrawLatex(0.02,0.95,m_PlotLabel.c_str());
549 
550  AddCanvas(can);
551  }
552 
553  void HistPlot::SetScale(double scale){
554  if(scale <= 0){
555  m_Scale = 1.;
556  m_SetScale = false;
557  m_ScaleLabel = "a. u.";
558  } else{
559  m_Scale = scale;
560  m_SetScale = true;
561  }
562  }
563 
564  void HistPlot::SetScaleLabel(const std::string& label){
565  m_ScaleLabel = label;
566  }
567 
568  void HistPlot::SetPlotLabel(const std::string& label){
569  m_PlotLabel = label;
570  }
571 
572  void HistPlot::SetPlotTitle(const std::string& title){
573  m_PlotTitle = title;
574  }
575 
576  void HistPlot::SetRebin(int rebin){
577  if(rebin > 0)
578  m_Rebin = rebin;
579  }
580 
581  void HistPlot::WriteHist(const std::string& name){
582  TFile* file = new TFile(name.c_str(),"UPDATE");
583  file->mkdir((GetName()+"/hist").c_str());
584  file->cd((GetName()+"/hist").c_str());
585  int N = m_1DHists.size();
586  for(int i = 0; i < N; i++)
587  m_1DHists[i]->Write("",TObject::kOverwrite);
588  N = m_2DHists.size();
589  for(int i = 0; i < N; i++)
590  m_2DHists[i]->Write("",TObject::kOverwrite);
591  file->Close();
592  delete file;
593  }
594 
595 }
std::string GetName() const
Returns object name.
Definition: RFBase.cc:104
std::string GetTitle() const
Returns object title.
Definition: RFBase.cc:108
virtual void Clear()
Clears RFBase of all connections to other objects.
Definition: RFPlot.cc:48
bool IsEmpty() const
Checks whether this is default (empty) instance of class.
Definition: RFBase.cc:84
virtual void Clear()
Clears RFBase of all connections to other objects.
Definition: HistPlot.cc:63