make_extrap_figure.C
Go to the documentation of this file.
1 // This macro uses the histograms produced by make_extrap_figure_hists.C
4 #include "CAFAna/Core/OscCurve.h"
9 using namespace ana;
10 
11 #include "OscLib/IOscCalc.h"
12 #include "Utilities/rootlogon.C"
13 
14 #include "TCanvas.h"
15 #include "TColor.h"
16 #include "TFile.h"
17 #include "TGraph.h"
18 #include "TH2.h"
19 #include "TLatex.h"
20 #include "TLegend.h"
21 #include "TLine.h"
22 #include "TPad.h"
23 #include "TStyle.h"
24 
25 #include <iostream>
26 #include <vector>
27 
28 // How far up the canvas the bottom set of plots extends
29 const double kVertSplit = .55;
30 // How close to the left and right edges the plots get
31 const double kLRMargin = .05;
32 const double kTopMargin = .05;
33 const double kBottomMargin = .1;
34 // Fraction of pad width
35 const double kMatrixWidth = .15;
36 // Twice this padding is added between all the panels
37 const double kXPadding = .003;
38 const double kYPadding = .005;
39 const Color_t kPred = kBlue;
40 //const Color_t kPred = kMagenta - 9;
41 //kViolet - 4;
42 //kMagenta - 9;
43 //const Color_t kMC = kAzure + 5;
44 const Color_t kMC = kRed;
45 
46 void no_titles(std::vector<TH1*> hs)
47 {
48  for(TH1* h: hs){
49  h->SetTitle("");
50  h->GetXaxis()->SetTitleSize(0);
51  h->GetYaxis()->SetTitleSize(0);
52  }
53 }
54 
55 TPad* pad(std::string name, double x0, double y0, double x1, double y1)
56 {
57  TPad* p = new TPad(name.c_str(), "", 0, 0, 1, 1);
58  p->SetLeftMargin(x0+kXPadding);
59  p->SetRightMargin(1-x1+kXPadding);
60  p->SetTopMargin(1-y1+kYPadding);
61  p->SetBottomMargin(y0+kYPadding);
62  p->SetFillStyle(0);
63  p->Draw();
64  return p;
65 }
66 
67 void latex(double x, double y, std::string txt, double ang = 0, int align = 22, double size = 0.042)
68 {
69  TLatex* ltx = new TLatex(x, y, txt.c_str());
70  ltx->SetNDC();
71  ltx->SetTextSize(size);
72  ltx->SetTextAlign(align);
73  ltx->SetTextAngle(ang);
74  ltx->Draw();
75 }
76 
77 TGraph* rotate_curve(TH1* h)
78 {
79  TGraph* g = new TGraph;
80  for(int i = 0; i < h->GetNbinsX()+2; ++i){
81  g->SetPoint(i, h->GetBinContent(i), h->GetXaxis()->GetBinCenter(i));
82  }
83  g->SetLineWidth(2);
84  g->SetLineColor(h->GetLineColor());
85  return g;
86 }
87 
88 TGraph* rotate_hist(TH1* h, bool neg = false)
89 {
90  int sign = neg ? -1 : +1;
91 
92  TGraph* g = new TGraph;
93  for(int i = 0; i < h->GetNbinsX()+2; ++i){
94  g->SetPoint(2*i, sign*h->GetBinContent(i), h->GetXaxis()->GetBinLowEdge(i));
95  g->SetPoint(2*i+1, sign*h->GetBinContent(i), h->GetXaxis()->GetBinUpEdge(i));
96  }
97  g->SetLineWidth(2);
98  g->SetLineColor(h->GetLineColor());
99 
100  return g;
101 }
102 
103 void reset_label(TH2* h, std::vector<int> bin)
104 {
105  for(int i:bin) h->GetXaxis()->ChangeLabel(i, -1, .0, -1, -1, -1, "");
106 /*
107  int nbin = h->GetNbinsX();
108  if(nbin<=2) return;
109  for(int i = 2; i< nbin - 1; i++) h->GetXaxis()->ChangeLabel(i, -1, .0, -1, -1, -1, "");
110 */
111 }
112 
114 {
116  calc->SetdCP(0.186 *(TMath::Pi()));
117  calc->SetTh23( asin(sqrt(0.586)));
118  calc->SetDmsq32( 2.502e-3);
119 
120  bool isNue = true;
121  if(ana == "numu") isNue = false;
122 
123  TFile* f = isNue?new TFile("extrap_hists_nue.root", "read"):new TFile("extrap_hists_numu.root", "read");
124 
125  gStyle->SetLabelSize(.0248, "xy");
126  gStyle->SetLabelOffset(0.006, "xy");
127 
128  std::string Extrap = isNue?"MEextrap":"MMextrap";
129 
130  TH1* h_reco_dat_nd = (TH1*)f->Get((TString)Extrap+"/NDDatReco");
131  TH1* h_reco_mc_nd = (TH1*)f->Get((TString)Extrap+"/NDMCReco");
132  h_reco_dat_nd->Scale(1, "width");
133  h_reco_mc_nd->Scale(1, "width");
134 
135  TH1* h_true_mc_nd = (TH1*)f->Get((TString)Extrap+"/NDMCTrue");
136  TH1* h_true_rwt_nd = (TH1*)f->Get((TString)Extrap+"/NDRwtTrue");
137 
138  TH1* h_true_mc_fd = (TH1*)f->Get((TString)Extrap+"/FDMCTrue");
139  TH1* h_true_ext_fd = (TH1*)f->Get((TString)Extrap+"/FDExtTrue");
140 
141  TH1* h_reco_ext_fd = (TH1*)f->Get((TString)Extrap+"/FDExtReco");
142  TH1* h_reco_mc_fd = (TH1*)f->Get((TString)Extrap+"/FDMCReco");
143 
144  TH1* h_nd_matrix = (TH1*)f->Get((TString)Extrap+"/NDRecoToTrue");
145 
146  TH1* h_fd_matrix = (TH1*)f->Get((TString)Extrap+"/FDTrueToReco");
147 //HACK
148  OscillatableSpectrum* trueToRecoMC = LoadFromFile<OscillatableSpectrum>(f, Extrap+"/trueToRecoMC").release();
149  OscillatableSpectrum* trueToRecoExt = LoadFromFile<OscillatableSpectrum>(f, Extrap+"/trueToRecoExt").release();
150 
151  int to = (ana=="nue")?12:14;
152  auto spec_matrixmc = trueToRecoMC->Oscillated(calc, 14, to);
153  auto spec_matrixext = trueToRecoExt->Oscillated(calc, 14, to);
154 
155  h_reco_ext_fd = spec_matrixmc.ToTH1(9e20);
156  h_reco_mc_fd = spec_matrixext.ToTH1(9e20);
157 
158  if(!isNue){
159  h_reco_ext_fd->Scale(1, "width");
160  h_reco_mc_fd->Scale(1, "width");
161  }
162 
163  TH1* h_osc = (TH1*)f->Get("osc");
164 
165  TH1* h_fn = (TH1*)h_true_mc_fd->Clone();
166  h_fn->Divide(h_true_mc_nd);
167 
168  h_true_mc_fd->Multiply(h_osc);
169  h_true_ext_fd->Multiply(h_osc);
170 
171  h_reco_dat_nd->Scale(1e-5);
172  h_reco_mc_nd->Scale(1e-5);
173 
174  h_true_mc_nd->Scale(1e-5);
175  h_true_rwt_nd->Scale(1e-5);
176 
177  no_titles({h_reco_dat_nd, h_reco_mc_nd, h_true_mc_nd, h_true_rwt_nd,
178  h_true_mc_fd, h_true_ext_fd, h_reco_ext_fd, h_reco_mc_fd,
179  h_nd_matrix, h_fd_matrix,
180  h_osc, h_fn});
181 
182 
183  TCanvas* c = new TCanvas("a", "b", 1800, 600);
184 
185  TPad* pNDReco = pad("NDReco", kLRMargin, kVertSplit, kLRMargin+kMatrixWidth, 1-kTopMargin);
186  TPad* pFDReco = pad("FDReco", 1-kLRMargin-kMatrixWidth, kVertSplit, 1-kLRMargin, 1-kTopMargin);
187 
188  TPad* pNDMatrix = pad("NDMatrix", kLRMargin, kBottomMargin, kLRMargin+kMatrixWidth, kVertSplit);
189  TPad* pFDMatrix = pad("FDMatrix", 1-kLRMargin-kMatrixWidth, kBottomMargin, 1-kLRMargin, kVertSplit);
190 
191  const double kNDTrueEdge = (.5+(kLRMargin+kMatrixWidth))/2.;
192  const double kCenterPlotWidth = .5-kNDTrueEdge;
193 
194  TPad* pNDTrue = pad("NDTrue", kLRMargin+kMatrixWidth, kBottomMargin, kNDTrueEdge, kVertSplit);
195  TPad* pFDTrue = pad("FDTrue", 1-kNDTrueEdge, kBottomMargin, 1-kLRMargin-kMatrixWidth, kVertSplit);
196 
197  TPad* pFN = pad("FN", kNDTrueEdge, kBottomMargin, .5, kVertSplit);
198  TPad* pOsc = pad("Osc", .5, kBottomMargin, 1-kNDTrueEdge, kVertSplit);
199 
200 
201  Int_t myPalette[1000];
202 
203  Double_t r[] = {1,.70, .2};
204  Double_t g[] = {1,.85, .15};
205  Double_t b[] = {1,.90, .35};
206  Double_t stop[] = {0., .15, 1.0};
207  Int_t FI = TColor::CreateGradientColorTable(3, stop, r, g, b, 1000);
208  for (int i=0;i<1000;i++) myPalette[i] = FI+i;
209  gStyle->SetPalette(1000, myPalette);
210 
211  c->cd();
212 
213  pNDReco->cd();
214  h_reco_mc_nd->SetLineColor(kMC);
215  h_reco_mc_nd->GetXaxis()->SetLabelSize(0);
216  h_reco_dat_nd->Draw("ep");
217  h_reco_dat_nd->GetYaxis()->SetRangeUser(0,1.3*h_reco_dat_nd->GetMaximum());
218  h_reco_mc_nd->Draw("hist same");
219 
220  pFDReco->cd();
221  h_reco_mc_fd->SetLineColor(kMC);
222  h_reco_mc_fd->GetXaxis()->SetLabelSize(0);
223  h_reco_mc_fd->Draw("hist Y+");
224  h_reco_mc_fd->GetYaxis()->SetRangeUser(0,1.3*h_reco_mc_fd->GetMaximum());
225  h_reco_ext_fd->SetLineColor(kPred);
226  h_reco_ext_fd->Draw("hist same");
227 
228 
229  double biny = isNue? 5: 4.75;
230 
231  pNDTrue->cd();
232  h_true_mc_nd->SetLineColor(kMC);
233  h_true_rwt_nd->SetLineColor(kPred);
234  TH2* axes = new TH2F("", "", 100, 0, 1, 100, 0.5, biny);
235  axes->GetYaxis()->SetLabelSize(0);
236  axes->Draw();
237  rotate_hist(h_true_mc_nd)->Draw("l same");
238  rotate_hist(h_true_rwt_nd)->Draw("l same");
239  axes->GetXaxis()->ChangeLabel(1, -1, .0248, 1);
240  reset_label(axes,{2,3,4,5});
241 
242 
243  pFDTrue->cd();
244  h_true_mc_fd->SetLineColor(kMC);
245  h_true_ext_fd->SetLineColor(kPred);
246  axes = new TH2F("", "", 100,-12, 0, 100, 0.5, biny);
247  axes->GetYaxis()->SetLabelSize(0);
248  axes->Draw();
249  rotate_hist(h_true_mc_fd, true)->Draw("l same");
250  rotate_hist(h_true_ext_fd, true)->Draw("l same");
251  axes->GetXaxis()->ChangeLabel(1, -1, .0248, 1, -1, -1, "12");
252  reset_label(axes,{2,3,4,5,6});
253 
254 
255  Int_t myPalette1[1000];
256 
257  Double_t r1[] = {1,1, .5};
258  Double_t g1[] = {1,0, .0};
259  Double_t b1[] = {1,.0, .0};
260  Double_t stop1[] = {0., .25, 1.0};
261  Int_t FI1 = TColor::CreateGradientColorTable(3, stop1, r1, g1, b1, 1000);
262  for (int i=0;i<1000;i++) myPalette1[i] = FI1+i;
263  gStyle->SetPalette(1000, myPalette1);
264 
265  pNDMatrix->cd();
266  axes = new TH2F("", "", 100, 0, 5, 100, 0.5, 4.75);
267  h_nd_matrix->GetXaxis()->SetRangeUser(0,5);
268  h_nd_matrix->GetYaxis()->SetRangeUser(0.5,4.75);
269  h_nd_matrix->Draw("col same");
270  axes->GetXaxis()->ChangeLabel(1, -1, .0248, 1);
271  axes->Draw("same");
272 
273  double recoE = isNue?27:4.75;
274 
275  pFDMatrix->cd();
276 
277  axes = new TH2F("", "", 100, 0, 5, 100, 0.5, recoE);
278  if(!isNue) h_fd_matrix->GetXaxis()->SetRangeUser(0,5);
279  h_fd_matrix->GetYaxis()->SetRangeUser(0.5,4.75);
280  h_fd_matrix->Draw("col Y+");
281  axes->GetXaxis()->ChangeLabel(1, -1, .0248, 1);
282  axes->Draw("Y+ same");
283 
284  double ratio = 2;
285  pFN->cd();
286  h_fn->SetLineColor(kMC);
287  h_fn->Scale(1000);
288  axes = new TH2F("", "", 100, 0, ratio, 100, 0.5, biny);
289  axes->GetYaxis()->SetLabelSize(0);
290  axes->Draw();
291  rotate_curve(h_fn)->Draw("l same");
292  axes->GetXaxis()->ChangeLabel(1, -1, .0248, 1);
293  axes->GetXaxis()->ChangeLabel(-1, -1, .0248, -1);
294  reset_label(axes,{2,3,4});
295 
296  pOsc->cd();
297  h_osc->SetLineColor(kMC);
298  axes = new TH2F("", "", 100, 0, isNue?0.102:1, 100, 0.5, biny);
299  axes->GetYaxis()->SetLabelSize(0);
300  axes->Draw();
301  rotate_curve(h_osc)->Draw("l same");
302  axes->GetXaxis()->ChangeLabel(1, -1, .0248, 1);
303  axes->GetXaxis()->ChangeLabel(-1, -1, .0248, -1);
304  reset_label(axes,{2,3,4,5});
305 
306  c->cd();
307 
308  latex(.4*kLRMargin, .5*kVertSplit + .5*(1-kTopMargin), "10^{5} ND Events/1 GeV", 90);
309 
310  latex(.4*kLRMargin, .5*kVertSplit + .5*(kBottomMargin), "True Energy (GeV)", 90);
311  latex(1-.50*kLRMargin, .5*kVertSplit + .5*(kBottomMargin), "True Energy (GeV)", 90);
312 
313  latex(kLRMargin+kMatrixWidth/2, .5*kBottomMargin, "ND Reco Energy (GeV)");
314  if(ana == "nue"){
315  latex(1-kLRMargin-kMatrixWidth/2, .5*kBottomMargin, "FD Analysis Bin");
316  latex(1-.50*kLRMargin, .5*kVertSplit + .5*(1-kTopMargin), "FD Events", 90);
317  }
318  else{
319  latex(1-kLRMargin-kMatrixWidth/2, .5*kBottomMargin, "FD Reco Energy (GeV)");
320  latex(1-.50*kLRMargin, .5*kVertSplit + .5*(1-kTopMargin), "FD Events/1 GeV", 90);
321  }
322  latex(kNDTrueEdge-.5*kCenterPlotWidth, .5*kBottomMargin, "10^{5} ND Events");
323  latex(1-kNDTrueEdge+.5*kCenterPlotWidth, .5*kBottomMargin, "FD Events");
324 
325  latex(.5-kCenterPlotWidth/2, .5*kBottomMargin, " 10^{-3} F/N Ratio");
326 
327  latex(.5+kCenterPlotWidth/2, .5*kBottomMargin, isNue?"P(#nu_{#mu}#rightarrow#nu_{e})":"P(#nu_{#mu}#rightarrow#nu_{#mu})");
328 
329 
330  // Draw legend by hand
331  TLine* lin = new TLine(.37, .85, .415, .85);
332  lin->SetLineWidth(2);
333  lin->Draw();
334  lin = new TLine(.37, .75, .415, .75);
335  lin->SetLineWidth(2);
336  lin->SetLineColor(kMC);
337  lin->Draw();
338  lin = new TLine(.37, .65, .415, .65);
339  lin->SetLineWidth(2);
340  lin->SetLineColor(kPred);
341  lin->Draw();
342 
343  latex(.42, .85, "ND data", 0, 12, 0.065);
344  latex(.42, .75, "Base Simulation", 0, 12, 0.065);
345  latex(.42, .65, "Data-Driven Prediction", 0, 12, 0.065);
346 
347  c->Print((TString)ana+"_extrap.pdf");
348 }
const XML_Char * name
Definition: expat.h:151
void latex(double x, double y, std::string txt, double ang=0, int align=22, double size=0.042)
Simple analyzer module to print to text file information required for c++ based block alignment study...
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const double kTopMargin
Float_t y1[n_points_granero]
Definition: compare.C:5
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
Float_t x1[n_points_granero]
Definition: compare.C:5
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH1 * ratio(TH1 *h1, TH1 *h2)
const Color_t kMC
const Color_t kPred
TGraph * rotate_hist(TH1 *h, bool neg=false)
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
virtual void SetDmsq32(const T &dmsq32)=0
void no_titles(std::vector< TH1 * > hs)
const double kXPadding
const double kYPadding
TPad * pad(std::string name, double x0, double y0, double x1, double y1)
const double kLRMargin
const double kVertSplit
float bin[41]
Definition: plottest35.C:14
TGraph * rotate_curve(TH1 *h)
void reset_label(TH2 *h, std::vector< int > bin)
Spectrum Oscillated(osc::IOscCalc *calc, int from, int to) const
virtual void SetTh23(const T &th23)=0
const hit & b
Definition: hits.cxx:21
const double kBottomMargin
TRandom3 r(0)
Float_t e
Definition: plot.C:35
Spectrum with true energy information, allowing it to be oscillated
void make_extrap_figure(std::string ana)
enum BeamMode kBlue
const double kMatrixWidth
def sign(x)
Definition: canMan.py:197
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string