Plot.cxx
Go to the documentation of this file.
1 #include "CherenkovCalc.h"
2 #include <vector>
3 
4 #include "TLatex.h"
5 #include "TCanvas.h"
6 #include "TH1.h"
7 #include "TH2.h"
8 #include "TH3.h"
9 #include "TLegend.h"
10 #include "TError.h"
11 #include "TStyle.h"
12 #include "TGaxis.h"
13 #include "TF1.h"
14 
15 /* Plotting code for LightTuning */
16 
17 // Base text size for drawing plots
18 static const double tsize = 0.05;
19 
20 static const int datacolor = kBlack;
21 static const int mccolor = kRed;
22 static const int flatpadbins = 1;
23 
24 // Plot the distributions as 2D color plots and save them to PDF files. Four
25 // plots are generated for each call: data and MC in X and Y. The top of the
26 // page includes 'title', as does the name of the output file.
27 void TwoDplots(CherenkovCalc & calc, const char * const title)
28 {
29  // Suppress excess printing. NOTE: persists outside this function because
30  // that's how ROOT is. I don't think I care in this case, but it's messy.
31  gErrorIgnoreLevel = kFatal;
32 
33  gStyle->SetOptStat(11);
34  gStyle->SetStatX(0.6);
35 
36  TH3D* const mc_x = (TH3D*)calc.mcX()->Clone("MCX");
37  TH3D* const mc_y = (TH3D*)calc.mcY()->Clone("MCY");
38  TH3D* const data_x = (TH3D*)calc.dataX()->Clone("DATAX");
39  TH3D* const data_y = (TH3D*)calc.dataY()->Clone("DATAY");
40 
41  TCanvas* c = new TCanvas("Fitter", "Fitter");
42 
43  TLatex ltitle(0.5, 0.995, title);
44  ltitle.SetNDC();
45  ltitle.SetTextAlign(23);
46  ltitle.SetTextFont(42);
47  ltitle.Draw();
48 
49  c->SetLogz();
50 
51  c->Divide(2,2, 0.01, 0.05);
52 
53  TH3D * const colhists[4] = { mc_x, mc_y, data_x, data_y };
54 
55  const int NRGBs = 7, NCont = 20;
56  gStyle->SetNumberContours(NCont);
57  Double_t stops[NRGBs] = { 0.00, 0.08, 0.25, 0.46, 0.60, 0.85, 1.00 };
58  Double_t red[NRGBs] = { 0.70, 0.00, 0.00, 0.00, 1.00, 1.00, 0.33 };
59  Double_t green[NRGBs] = { 0.70, 0.70, 0.30, 0.40, 1.00, 0.00, 0.00 };
60  Double_t blue[NRGBs] = { 0.70, 1.00, 1.00, 0.00, 0.00, 0.00, 0.00 };
61  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
62 
63  for(int i = 0; i < 4; i++){
64  c->cd(i+1);
65  gPad->SetTopMargin(0.155);
66  gPad->SetRightMargin(0.16);
67  colhists[i]->Project3D("yx")->Draw("colz");
68  }
69 
70  c->SaveAs(Form("%s.pdf", title));
71 
72  delete mc_x;
73  delete mc_y;
74  delete data_x;
75  delete data_y;
76  delete c;
77 }
78 
79 // Flatten the input 2D histogram 'h3', returning the result in 'h1'. The
80 // flattened histogram is normalized by column, i.e. the integral of each
81 // column (set of bins in the output) is adjusted to 1. Set 'max' to the
82 // highest value of any non-overflow bin. Overflow bins are not included in
83 // the flattened histogram. This is only used for drawing.
84 static TH1D * flattenTH3(const TH3D * const h3, double & max)
85 {
86  #ifdef VERBOSEFLATTENTH3
87  printf("Flattening %s\n", h3->GetName());
88  #endif
89 
90  const int nx = h3->GetNbinsX(), ny = h3->GetNbinsY(), nz = h3->GetNbinsZ();
91 
92  std::vector<double> bins;
93  {
94  double curr = 0; // current position along the flattened axis
95  for(int k = 1; k <= nz; k++){
96  for(int i = 1; i <= nx; i++){
97  for(int j = 1; j <= ny; j++){
98  bins.push_back(curr);
99  // since we're flattening onto Y (energy)
100  curr += h3->GetYaxis()->GetBinWidth(j);
101  }
102 
103  // Add a pad after each group, except at the very end
104  if(!(i == nx && k == nz)){
105  for(int j = 1; j <= flatpadbins; j++){
106  bins.push_back(curr);
107  curr += h3->GetYaxis()->GetBinWidth(h3->GetNbinsY()/2);
108  }
109  }
110  }
111  }
112  bins.push_back(curr);
113  }
114 
115  TH1D * h1 = new TH1D("h1", "", bins.size()-1, &bins[0]);
116 
117  max = 0;
118 
119  const double stdw = h3->GetYaxis()->GetBinWidth(h3->GetNbinsY()/2);
120 
121  for(int k = 1; k <= nz; k++){
122  for(int i = 1; i <= nx; i++){
123  double integral = 0;
124  for(int j = 1; j <= ny; j++) integral += h3->GetBinContent(i, j, k);
125 
126  for(int j = 1; j <= ny; j++){
127  #ifdef VERBOSEFLATTENTH3
128  printf("%3d %3d %3d: %7.0f\n", i, j, k, h3->GetBinContent(i, j, k));
129  #endif
130  const double val = integral == 0?0:h3->GetBinContent(i, j, k)/integral;
131  const double wcorr = stdw/h3->GetYaxis()->GetBinWidth(j);
132  if(val*wcorr > max) max = val*wcorr;
133  const int xzbin = (i-1) + (k-1)*nx;
134  h1->SetBinContent(xzbin*(ny+flatpadbins) + j, wcorr*val);
135  h1->SetBinError (xzbin*(ny+flatpadbins) + j,
136  integral == 0?0:wcorr*h3->GetBinError(i, j, k)/integral);
137  }
138  }
139  }
140 
141  // Hackery to avoid drawing artifacts. Remember that this histogram is
142  // only used for drawing.
143  for(int k = 1; k <= nz; k++){
144  for(int i = 1; i <= nx; i++){
145  for(int j = 0; j < flatpadbins; j++){
146  const int bin = (i-1)*(ny+flatpadbins) + j + ny + 1;
147  h1->SetBinContent(bin,
148  j+1 < flatpadbins? h1->GetBinContent(bin-1)
149  : h1->GetBinContent(bin+1));
150  h1->SetBinError(bin, 1);
151  }
152  }
153  }
154  return h1;
155 }
156 
157 // Converts the input histogram so that each bin is zero if it was zero
158 // and one if it wasn't. This is a hack to workaround the fact that my own
159 // flattenTH3 changes normalizations of the mask histograms, which is
160 // undesirable for plotting. The input histogram is modified in place and
161 // also returned.
162 static TH1D * maskify(TH1D * mask)
163 {
164  for(int i = 1; i <= mask->GetNbinsX(); i++)
165  mask->SetBinContent(i, !!mask->GetBinContent(i));
166  return mask;
167 }
168 
169 // Plot the distributions as 1D histograms by flatting them in column-major
170 // (distance bin-major) order. Save them to PDF files. Two plots are
171 // generated for each call: data and MC overlaid in X and Y. The top of the
172 // page includes 'title', as does the name of the output file.
173 void Flatplots(CherenkovCalc & calc, const char * const title1,
174  const char * const title2)
175 {
176  // See above comment on gErrorIgnoreLevel
177  gErrorIgnoreLevel = kFatal;
178 
179  gStyle->SetOptStat(0);
180 
181  double max[4] = { 0 };
182  double dummy = 0;
183 
184  const int nx = calc.mcX()->GetNbinsX(), ny = calc.mcX()->GetNbinsY(),
185  nz = calc.mcX()->GetNbinsZ();
186 
187  TH1D * flat_mc[2] = { flattenTH3(calc. mcX(), max[0]),
188  flattenTH3(calc. mcY(), max[1]) };
189  TH1D * flat_data[2] = { flattenTH3(calc.dataX(), max[2]),
190  flattenTH3(calc.dataY(), max[3]) };
191  TH1D * flat_mask[2] = { maskify(flattenTH3(calc.invmaskX(), dummy)),
192  maskify(flattenTH3(calc.invmaskY(), dummy)) };
193 
194  const double Max = std::max(std::max(max[0],max[1]), std::max(max[2],max[3]));
195 
196  for(int i = 0; i < 2; i++){
197  TCanvas * c = new TCanvas("flat", "flat");
198 
199  const double divheight = 0.3;
200  const double textratio = (1-divheight)/divheight;
201 
202  TPad * toppad = new TPad("top", "top", 0, divheight, 1, 1);
203  TPad * botpad = new TPad("liv", "liv", 0, 0, 1, divheight);
204 
205  const double leftmargin = 0.075, rightmargin = 0.02;
206 
207  toppad->SetRightMargin(rightmargin);
208  toppad->SetLeftMargin(leftmargin);
209  toppad->SetTopMargin(0.07);
210  toppad->SetBottomMargin(0);
211 
212  botpad->SetTopMargin(0);
213  botpad->SetBottomMargin(0.085/divheight);
214  botpad->SetRightMargin(rightmargin);
215  botpad->SetLeftMargin(leftmargin);
216 
217  botpad->Draw();
218  toppad->Draw();
219 
220  toppad->cd();
221 
222  flat_mc[i]->SetLineColor(mccolor);
223  flat_mc[i]->SetMarkerColor(mccolor);
224  flat_mc[i]->SetLineWidth(nz>1?1:2);
225 
226  flat_data[i]->SetLineColor(datacolor);
227  flat_data[i]->SetMarkerColor(datacolor);
228 
229  TH1D * firsthist = flat_mc[i];
230 
231  const double maxy = Max*1.14;
232 
233  firsthist->GetYaxis()->SetRangeUser(0.0001, maxy);
234  firsthist->GetXaxis()->SetTitle("PE/length in cell (cm^{-1})");
235  firsthist->GetYaxis()->SetTitle("Probability normalized by distance bin");
236  firsthist->GetXaxis()->SetTitleSize(tsize);
237  firsthist->GetXaxis()->SetLabelSize(0);
238  firsthist->GetYaxis()->SetTitleSize(tsize);
239  firsthist->GetYaxis()->SetLabelSize(tsize);
240  firsthist->GetXaxis()->CenterTitle();
241  firsthist->GetYaxis()->CenterTitle();
242  firsthist->GetYaxis()->SetTickSize(nz > 1?0.005:0.01);
243  firsthist->GetXaxis()->SetTickSize(0);
244  firsthist->GetYaxis()->SetDecimals();
245  firsthist->GetYaxis()->SetTitleOffset(leftmargin*10.);
246 
247  firsthist->GetXaxis()->SetNdivisions(-(100*ny+nx));
248 
249  flat_mc[i]->Draw("histe");
250  flat_data[i]->Draw("samee");
251 
252  const double ylo = calc.mcX()->GetYaxis()->GetBinLowEdge(1);
253  const double yhi = calc.mcX()->GetYaxis()->GetBinLowEdge(calc.mcX()->GetNbinsY()+1);
254 
255  TH1D * borders = (TH1D *)firsthist->Clone("borders");
256  borders->Reset();
257  borders->SetFillStyle(1001);
258  borders->SetFillColor(kBlack);
259  borders->SetLineColor(kBlack);
260  borders->SetLineWidth(0);
261 
262  for(int z = 1; z <= nz; z++){
263  for(int x = 1; x <= nx; x++){
264  const int xzbin = (x-1) + (z-1)*nx + 1;
265  const double text_x_pos =
266  firsthist->GetBinCenter( xzbin *(ny+flatpadbins)+1)/2 +
267  firsthist->GetBinCenter((xzbin-1)*(ny+flatpadbins)+1)/2;
268  TLatex * l = new TLatex(text_x_pos, maxy*0.83,
269  Form("#splitline{l: %.0f}{to %.0fcm}",
270  calc.mcX()->GetXaxis()->GetBinLowEdge(x),
271  calc.mcX()->GetXaxis()->GetBinLowEdge(x+1)));
272  l->SetTextFont(42);
273  l->SetTextAlign(22);
274  l->SetTextSize(tsize*0.5);
275  l->SetTextAngle(270);
276  l->Draw();
277 
278  if(nz > 1 && x == 1){
279  TLatex * l = new TLatex(text_x_pos, maxy*0.94,
280  Form("#splitline{w: %.1f}{to %.1fm}",
281  calc.mcX()->GetZaxis()->GetBinLowEdge(z)/100,
282  calc.mcX()->GetZaxis()->GetBinLowEdge(z+1)/100));
283  l->SetTextFont(42);
284  l->SetTextAlign(22);
285  l->SetTextSize(tsize*0.5);
286  l->SetTextAngle(270);
287  l->Draw();
288  }
289 
290  for(int i = 0; i < flatpadbins; i++)
291  borders->SetBinContent((ny+flatpadbins)*xzbin+i-flatpadbins+1, maxy);
292  }
293  }
294 
295  gStyle->SetGridColor(kBlack);
296 
297  TLatex ltitle(0.5, 0.993,
298  Form("%s %c %s - black is data", title1, i?'y':'x', title2));
299  ltitle.SetNDC();
300  ltitle.SetTextAlign(23);
301  ltitle.SetTextFont(42);
302  ltitle.Draw();
303 
304  flat_mask[i]->Scale(maxy);
305  flat_mask[i]->SetFillStyle(1001);
306  flat_mask[i]->SetFillColorAlpha(kGray, 0.5);
307  flat_mask[i]->SetLineColorAlpha(kGray, 0.5);
308  flat_mask[i]->Draw("histsame");
309 
310  borders->Draw("samehist");
311 
312  ////////////////////////////////////////////////////////////////////
313 
314  TH1D * flat_d = (TH1D*)flat_data[i]->Clone("flat_d");
315 
316  flat_d->Divide(flat_mc[i]);
317 
318  botpad->cd();
319  flat_d->Draw("e");
320  flat_d->GetXaxis()->SetLabelSize(textratio*tsize);
321  flat_d->GetYaxis()->SetLabelSize(textratio*tsize);
322  flat_d->GetXaxis()->SetTitleSize(textratio*tsize);
323  flat_d->GetYaxis()->SetTitleSize(textratio*tsize);
324  flat_d->GetXaxis()->SetLabelSize(0);
325  flat_d->GetXaxis()->SetTitle("PE/length in cell (cm^{-1})");
326  flat_d->GetYaxis()->SetTitle("Ratio");
327  flat_d->GetXaxis()->CenterTitle();
328  flat_d->GetYaxis()->CenterTitle();
329  flat_d->GetYaxis()->SetTitleOffset(leftmargin*10./textratio);
330  flat_d->GetYaxis()->SetTickSize(0.01);
331  flat_d->GetYaxis()->SetNdivisions(508);
332  flat_d->GetXaxis()->SetTickSize(0);
333 
334  TF1 perfect("perfect", "1", 0, firsthist->GetBinLowEdge(firsthist->GetNbinsX()+1));
335 
336  perfect.SetLineColor(kBlack);
337  perfect.SetLineWidth(1);
338  perfect.Draw("same");
339 
340  const double bot_ylo = 0.5, bot_yhi = 1.5;
341  flat_d->GetYaxis()->SetRangeUser(bot_ylo, bot_yhi);
342 
343  for(int k = 0; k < nz; k++){
344  for(int i = 0; i < nx; i++){
345  const int xzbin = i + nx*k;
346  TGaxis * a = new TGaxis(
347  firsthist->GetBinLowEdge((ny+flatpadbins)*xzbin+1) , bot_ylo,
348  firsthist->GetBinLowEdge((ny+flatpadbins)*xzbin+1+ny), bot_ylo,
349  ylo, yhi);
350  a->SetNdivisions(4);
351  a->SetLabelFont(42);
352  a->SetLabelSize(tsize*textratio*0.6666666);
353  a->SetTickSize(3);
354  a->SetTickLength(3);
355  a->Draw();
356  }
357  }
358 
359  TH1D * bot_borders = (TH1D *)flat_d->Clone("bot_borders");
360  bot_borders->Reset();
361  bot_borders->SetFillStyle(1001);
362  bot_borders->SetFillColor(kBlack);
363  bot_borders->SetLineColor(kBlack);
364  bot_borders->SetLineWidth(0);
365 
366  for(int z = 1; z <= nz; z++){
367  for(int x = 1; x <= nx; x++){
368  for(int i = 0; i < flatpadbins; i++){
369  const int xzbin = x + nx*(z-1);
370  bot_borders->SetBinContent((ny+flatpadbins)*xzbin+i-flatpadbins+1, bot_yhi);
371  }
372  }
373  }
374 
375  bot_borders->Draw("samehist");
376 
377  TH1D * bot_mask = (TH1D *)flat_mask[i]->Clone("bot_mask");
378  bot_mask->Scale(bot_yhi/maxy);
379  bot_mask->Draw("same hist");
380 
381  c->SaveAs(Form("%s_%c_%s.pdf", title1, i?'y':'x', title2));
382  delete flat_d;
383  delete toppad;
384  delete botpad;
385  delete c;
386  }
387  delete flat_mc[0];
388  delete flat_mc[1];
389  delete flat_data[0];
390  delete flat_data[1];
391 }
T max(const caf::Proxy< T > &a, T b)
static const int datacolor
Definition: Plot.cxx:20
enum BeamMode kRed
static const int flatpadbins
Definition: Plot.cxx:22
int NCont
Definition: rootlogon.py:78
double maxy
TH1F * h3
Definition: berger.C:36
osc::OscCalcDumb calc
static const double tsize
Definition: Plot.cxx:18
tuple blue
Definition: rootlogon.py:65
const double a
printf("%d Experimental points found\n", nlines)
const double j
Definition: BetheBloch.cxx:29
const Binning bins
float bin[41]
Definition: plottest35.C:14
z
Definition: test.py:28
static TH1D * maskify(TH1D *mask)
Definition: Plot.cxx:162
TH1F * h1
TH3D * dataX()
static const int mccolor
Definition: Plot.cxx:21
int NRGBs
Definition: rootlogon.py:77
static TH1D * flattenTH3(const TH3D *const h3, double &max)
Definition: Plot.cxx:84
void Flatplots(CherenkovCalc &calc, const char *const title1, const char *const title2)
Definition: Plot.cxx:173
void TwoDplots(CherenkovCalc &calc, const char *const title)
Definition: Plot.cxx:27
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:66
TH3D * dataY()