fd_plot_bless.C
Go to the documentation of this file.
7 #include "CAFAna/Core/Ratio.h"
9 #include "CAFAna/Core/Spectrum.h"
12 
14 #include "CAFAna/Vars/Vars.h"
15 #include "CAFAna/Vars/NueVars.h"
16 #include "CAFAna/Analysis/Style.h"
18 using namespace ana;
19 
21 
22 #include "TCanvas.h"
23 #include "TFile.h"
24 #include "TGraphAsymmErrors.h"
25 #include "TLatex.h"
26 #include "TLine.h"
27 #include "TLegend.h"
28 #include "TPad.h"
29 
30 //make a pretty legend
31 void Legend(double x0, double y0, double x1, double y1)
32 {
33  // Doesn't handle log-y well
34  // TLegend* leg = AutoPlaceLegend(x1-x0, y1-y0);
35  TLegend* leg = new TLegend(x0, y0, x1, y1);
36  leg->SetFillStyle(0);
37 
38  TH1* dummy = new TH1F("", "", 1, 0, 1);
39  TH1* dummyFill = new TH1F("", "", 1, 0, 1);
40  dummyFill->SetLineColor(kTotalMCColor);
41  dummyFill->SetFillColor(kTotalMCErrorBandColor);
42 
43  dummy->SetMarkerStyle(kFullCircle);
44 
45  leg->AddEntry(dummy->Clone(), "FD data", "lep");
46  dummy->SetLineColor(kGray+2);
47  dummy->SetMarkerColor(kGray+2);
48  // leg->AddEntry(dummy->Clone(), "ND data stagger", "lep");
49  dummy->SetLineColor(kTotalMCColor);
50  leg->AddEntry(dummy->Clone(), "Best Fit Pred.", "l");
51  // dummy->SetLineColor(kNueSignalColor);
52  // leg->AddEntry(dummy->Clone(), "Signal", "l");
53  // dummy->SetLineColor(kNCBackgroundColor);
54  // leg->AddEntry(dummy->Clone(), "NC", "l");
55  // dummy->SetLineColor(kBeamNueBackgroundColor);
56  // leg->AddEntry(dummy->Clone(), "Beam #nu_{e} CC", "l");
57  // dummy->SetLineColor(kNumuBackgroundColor);
58  // leg->AddEntry(dummy->Clone(), "#nu_{#mu} CC", "l");
59  // dummy->SetLineColor(kGray+1);
60  // leg->AddEntry(dummy->Clone(), "Cosmic", "l");
61  //dummy->SetFillColor(kNCBackgroundColor-9);
62  dummy->SetLineColor(kNCBackgroundColor);
63  FillWithDimColor(dummy);
64  leg->AddEntry(dummy->Clone(), "Total Background", "bf");
65 
66 
67  leg->Draw();
68 }
69 
70 //create a nova prelimary tag
72 {
73  TLatex* prelim = new TLatex(.935, .95, "NOvA Preliminary");
74  prelim->SetTextColor(kBlue);
75  prelim->SetNDC();
76  prelim->SetTextSize(2/30.);
77  prelim->SetTextAlign(32);
78  prelim->Draw();
79 }
80 
82 {
83  struct HistDef
84  {
86  };
87 
88  //which variables do you actually want to look at?
89  const int kNumVars =15;
90  const HistDef defs[kNumVars] = {
91  "nhits",
92  "inelast",
93  "vtxx",
94  "vtxy",
95  "vtxz",
96  "shw_len",
97  "shw_width",
98  "recoE",
99  "nshw",
100  "shw_e",
101  "had_e",
102  "ptp",
103  "cvne",
104  "cvne_zoom",
105  "costheta"
106  };
107 
108  //which cut tiers do you want to look at?
109  const int kNumSels = 2;
110  const std::string selNames[kNumSels] = {"presel","cvneloose"};
111 
113  IPrediction* preds[kNumSels][kNumVars];
114  Spectrum* cosms[kNumSels][kNumVars];
115 
116  //set oscillation parameters:
119 
120  calc->SetDmsq32(2.44e-3);
121  calc->SetTh13(asin(sqrt(.085))/2);
122  calc->SetTh23(asin(sqrt(0.45)));
123  calc->SetdCP(1.59*TMath::Pi());
124 
125  //pick up spectra/prediction objects from the input root file
126  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
127  const char* selName = selNames[selIdx].c_str();
128  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
129  const char* varName = defs[varIdx].name.c_str();
130 
131  std::cout<<TString::Format("%s/spect_%s", selName, varName).Data()<<std::endl;
132  spects[selIdx][varIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/spect_%s", selName, varName).Data()).release();
133  preds[selIdx][varIdx] = LoadFromFile<IPrediction>(fname, TString::Format("%s/pred_%s", selName, varName).Data()).release();
134  cosms[selIdx][varIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/cosm_%s", selName, varName).Data()).release();
135 
136 
137  if(selIdx == 0 && varIdx == 0)
138  std::cout << "data, MC POT: "
139  << spects[selIdx][varIdx]->POT() << " "
140  << preds[selIdx][varIdx]->Predict(&noosc).POT() << std::endl;
141 
142  } // end for varIdx
143  } // end for selIdx
144 
145  //loop over the different selector/cut levels to plot each one
146  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
147  const std::string selName = selNames[selIdx];
148  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
149 
150  const std::string varName = defs[varIdx].name;
151 
152  //play with the binning choice for given variables:
153  int rebinFactor=5;
154  if(varName=="nshw"){
155  rebinFactor=1;
156  }
157 
158  if(varName=="recoE"){
159  rebinFactor=4;
160  }
161 
162  if(varName=="shw_width"){
163  rebinFactor=1;
164  }
165 
166  if(varName=="shw_e"){
167  rebinFactor=2;
168  }
169 
170  if(varName=="had_e"){
171  rebinFactor=2;
172  }
173 
174  if(varName=="vtxx"||varName=="vtxy"||varName=="vtxz"){
175  rebinFactor=5;
176  }
177 
178  //create a canvas
179  TCanvas *c2 =new TCanvas("c2","c2",1000,800);
180  c2->SetLeftMargin(0.15); //default is 10%
181  c2->SetRightMargin(0.05);
182 
183  gPad->SetFillStyle(0);
184 
185  double legendxlow=0.6;
186  double legendxhigh=0.9;
187 
188  //Pull down the variable we want
189  Spectrum* spect = spects[selIdx][varIdx];
190  Spectrum* cosm = cosms[selIdx][varIdx];
191 
192  IPrediction* pred = preds[selIdx][varIdx];
193 
194  //Get components of the prediction and store them in histograms
195  TH1* hMC = pred->Predict(calc).ToTH1(spect->POT());
196  hMC->SetLineColor(kTotalMCColor);
197  hMC->Rebin(rebinFactor);
198 
199  TH1* hNC = pred->PredictComponent(calc,
201  Current::kNC,
202  Sign::kBoth).ToTH1(spect->POT());
203  hNC->SetLineColor(kNCBackgroundColor);
204  hNC->Rebin(rebinFactor);
205 
206  TH1* hCC = pred->PredictComponent(calc,
208  Current::kCC,
209  Sign::kBoth).ToTH1(spect->POT());
210  hCC->SetLineColor(kNumuBackgroundColor);
211  hCC->Rebin(rebinFactor);
212 
213  TH1* hTau = pred->PredictComponent(calc,
215  Current::kCC,
216  Sign::kBoth).ToTH1(spect->POT());
217  hTau->SetLineColor(kNumuBackgroundColor);
218  hTau->Rebin(rebinFactor);
219 
220  TH1* hBeam = pred->PredictComponent(calc,
222  Current::kCC,
223  Sign::kBoth).ToTH1(spect->POT());
224  hBeam->SetLineColor(kBeamNueBackgroundColor);
225  hBeam->Rebin(rebinFactor);
226 
227  TH1* hSig = pred->PredictComponent(calc,
229  Current::kCC,
230  Sign::kBoth).ToTH1(spect->POT());
231  hSig->SetLineColor(kNueSignalColor);
232  hSig->Rebin(rebinFactor);
233 
234  TH1* hData = spect->ToTH1(spect->POT());
235  hData->SetMarkerStyle(kFullCircle);
236  hData->Rebin(rebinFactor);
237 
238  TH1* hCos = cosm->ToTH1(spect->Livetime(), kGray+1, kSolid,
240  hCos->SetLineColor(kGray+2);
241  hCos->Rebin(rebinFactor);
242 
243  hSig->Scale(21.94/22.78);
244  hCC->Scale(0.71/0.62);
245  hBeam->Scale(3.11/3.02);
246  hNC->Scale(3.71/3.36);
247  hTau->Scale(1./1.);
248 
249  TH1 *hnew = (TH1*)hSig->Clone("hnew");
250  hMC=hnew;
251 
252  TH1 *hBack = (TH1*)hCC->Clone("hBack");
253 
254  hMC->SetLineColor(kRed);
255  hMC->Add(hCC);
256  hMC->Add(hNC);
257  hMC->Add(hBeam);
258  hMC->Add(hTau);
259  hMC->Add(hCos);
260 
261  hBack->Add(hNC);
262  hBack->Add(hBeam);
263  hBack->Add(hTau);
264  hBack->Add(hCos);
265  hBack->SetLineColor(kNCBackgroundColor);//+3);
266  FillWithDimColor(hBack);
267  //hBack->SetFillColor(kNCBackgroundColor-9);//+1);
268 
269  std::cout<<hBack->Integral()<<std::endl;
270  std::cout<<hMC->Integral()<<std::endl;
271 
272  hMC->GetXaxis()->CenterTitle();
273  hMC->GetYaxis()->CenterTitle();
274  //hMC->GetXaxis()->SetDecimals();
275  //hMC->GetYaxis()->SetDecimals();
276  hMC->GetYaxis()->SetTitleOffset(1.15);
277  //hMC->GetYaxis()->SetTitle(TString::Format("Events / %.02lf #times 10^{20} POT", /1E20));
278  hMC->GetYaxis()->SetTitle("Events / 6.05 #times 10^{20} POT-equiv");
279 
280  if(varName=="had_e"){
281  hMC->GetXaxis()->SetTitle("Hadronic calorimetric energy (GeV)");
282  }
283  if(varName=="shw_e"){
284  hMC->GetXaxis()->SetTitle("Primary shower calorimetric energy (GeV)");
285  }
286  if(varName=="recoE"){
287  hMC->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
288  }
289  if(varName=="vtxx"){
290  hMC->GetXaxis()->SetTitle("Vertex X (cm)");
291  }
292  if(varName=="vtxy"){
293  hMC->GetXaxis()->SetTitle("Vertex Y (cm)");
294  }
295  if(varName=="vtxz"){
296  hMC->GetXaxis()->SetTitle("Vertex Z (cm)");
297  }
298  if(varName=="costheta"){
299  hMC->GetXaxis()->SetTitle("cos#theta");
300  }
301  if(varName=="cvne"){
302  hMC->GetXaxis()->SetTitle("CVN #nu_{e} classifier");
303  }
304  if(varName=="cvne_zoom"){
305  hMC->GetXaxis()->SetTitle("CVN #nu_{e} classifier");
306  }
307 
308  if(varName=="inelast"){
309  hMC->GetXaxis()->SetTitle("Reconstructed inelasticity");
310  }
311 
312  //hMC->GetXaxis()->SetNdivisions(405,kFALSE);
313  hMC->SetMaximum(1.1*(hData->GetMaximum()+sqrt(hData->GetMaximum())));
314  //Draw the error band, MC, Data, MC components, and then Data one last time
315  if(varName=="vtxx"||varName=="vtxy"||varName=="vtxz"){
316  hMC->SetMaximum(1.9*hMC->GetMaximum());
317  }
318 
319  if(varName=="cvne_zoom"||varName=="costheta"){
320  legendxlow=0.2;
321  legendxhigh=0.60;
322  }
323  hMC ->SetMinimum(.001);
324 
325  if(varName=="shw_width"){
326  hMC->GetXaxis()->SetRangeUser(0,15);
327  }
328 
329  if(varName=="recoE"){
330  hMC->GetXaxis()->SetRangeUser(0,4);
331  }
332 
333  if(varName=="had_e"){
334  hMC->GetXaxis()->SetRangeUser(0,2);
335  }
336 
337  if(varName=="shw_e"){
338  hMC->GetXaxis()->SetRangeUser(0,4);
339  }
340  gPad->SetFrameLineColor(kGreen);
341  gPad->SetFrameLineWidth(1);
342 
343 
344  hMC ->Draw("hist");
345  //hBack ->Draw("hist");
346 
347  TGraphAsymmErrors* gr = GraphWithPoissonErrors(hData, hMC, false, true);
348  gr->SetMarkerStyle(kFullCircle);
349  gr->SetLineWidth(2);
350 
351  //gr->Draw("ep same");
352 
353  hBack->Draw("hist same");
354 
355  // hNC->Draw("hist same");
356  // hCC->Draw("hist same");
357  // hBeam->Draw("hist same");
358  // hCos->Draw("hist same");
359  // hSig->Draw("hist same");
360  hMC->Draw("hist same");
361  gr->Draw("ep same");
362 
363  hMC->Draw("same axis");
364  //Fix up the axis
365  c2->RedrawAxis();
366  c2->Update();
367 
368  //Add the legend and preliminary
369  //Legend(legendxlow, .5, legendxhigh, .85);
370  Legend(legendxlow, .7, legendxhigh, .85);
371  Preliminary();
372 
373  c2->Print(("plots/"+selName+"_"+varName+"_bless.pdf").c_str());
374  c2->Print(("plots/"+selName+"_"+varName+"_bless.eps").c_str());
375  c2->Print(("plots/"+selName+"_"+varName+"_bless.png").c_str());
376  c2->Print(("plots/"+selName+"_"+varName+"_bless.gif").c_str());
377  c2->Print(("plots/"+selName+"_"+varName+"_bless.C").c_str());
378 
379  c2->SetLogy();
380  hMC ->SetMinimum(0.1);
381  hMC ->SetMaximum(10*hMC ->GetMaximum());
382  hMC ->DrawCopy("hist");
383 
384  gr->Draw("ep same");
385  hBack->Draw("hist same");
386  // hNC->Draw("hist same");
387  // hCC->Draw("hist same");
388  // hBeam->Draw("hist same");
389  // hCos->Draw("hist same");
390  // hSig->Draw("hist same");
391  hMC->Draw("hist same");
392  gr->Draw("ep same");
393 
394  //Fix up the axis
395  c2->RedrawAxis();
396  c2->Update();
397 
398  //Add the legend and preliminary
399  //Legend(legendxlow, .6, legendxhigh, .85);
400  Legend(legendxlow, .75, legendxhigh, .85);
401  Preliminary();
402 
403 
404  c2->Print(("plots/"+selName+"_"+varName+"_bless_log.pdf").c_str());
405  c2->Print(("plots/"+selName+"_"+varName+"_bless_log.eps").c_str());
406  c2->Print(("plots/"+selName+"_"+varName+"_bless_log.png").c_str());
407  c2->Print(("plots/"+selName+"_"+varName+"_bless_log.gif").c_str());
408  c2->Print(("plots/"+selName+"_"+varName+"_bless_log.C").c_str());
409 
410  } // end for varIdx
411  } // end for selIdx
412 
413 }
Pass neutrinos through unchanged.
const int kNumVars
Definition: vars.h:14
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
osc::IOscCalculatorAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
Float_t y1[n_points_granero]
Definition: compare.C:5
(&#39; appearance&#39;)
Definition: IPrediction.h:18
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:553
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
Float_t x1[n_points_granero]
Definition: compare.C:5
(&#39;beam &#39;)
Definition: IPrediction.h:15
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual void SetdCP(const T &dCP)=0
TGraphAsymmErrors * GraphWithPoissonErrors(const TH1 *h, bool noErrorsXaxis, bool drawEmptyBins)
Calculate statistical errors appropriate for small Poisson numbers.
Definition: Plots.cxx:842
const Color_t kTotalMCErrorBandColor
Definition: Style.h:18
void FillWithDimColor(TH1 *h, bool usealpha)
const Color_t kNumuBackgroundColor
Definition: Style.h:31
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
const Color_t kNueSignalColor
Definition: Style.h:19
c2
Definition: demo5.py:33
Charged-current interactions.
Definition: IPrediction.h:39
General interface to any calculator that lets you set the parameters.
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
virtual void SetTh23(const T &th23)=0
const HistDef defs[kNumVars]
Definition: vars.h:15
virtual void SetDmsq32(const T &dmsq32)=0
const int kNumSels
Definition: vars.h:43
osc::NoOscillations noosc
Definition: NueSystFuncs.h:44
TLatex * prelim
Definition: Xsec_final.C:133
double POT() const
Definition: Spectrum.h:263
std::string name
Definition: NuePlotLists.h:12
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void fd_plot_bless(std::string fname)
Definition: fd_plot_bless.C:81
Neutral-current interactions.
Definition: IPrediction.h:40
void Preliminary()
Definition: fd_plot_bless.C:71
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
virtual void SetTh13(const T &th13)=0
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const Color_t kTotalMCColor
Definition: Style.h:17
void spects(osc::IOscCalculatorAdjustable *calc, IPrediction *predFHC, IPrediction *predRHC)
Definition: sensitivity.C:100
All neutrinos, any flavor.
Definition: IPrediction.h:26
const std::string selNames[kNumSels]
Definition: vars.h:46
Float_t e
Definition: plot.C:35
const Color_t kNCBackgroundColor
Definition: Style.h:22
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:266
T asin(T number)
Definition: d0nt_math.hpp:60
void Legend(double x0, double y0, double x1, double y1)
Definition: fd_plot_bless.C:31