plot_fd_datamc.C
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TH1.h"
3 #include "TH2.h"
4 #include "TFile.h"
5 #include "TLegend.h"
6 #include "TLatex.h"
7 #include "TStyle.h"
8 #include "TAttMarker.h"
9 #include "TArrow.h"
10 #include "TMarker.h"
11 
14 #include "CAFAna/Core/Cut.h"
16 #include "CAFAna/Analysis/Calcs.h"
17 #include "CAFAna/Analysis/Plots.h"
18 #include "CAFAna/Analysis/Style.h"
20 #include "CAFAna/Core/Spectrum.h"
21 #include "CAFAna/Cuts/SpillCuts.h"
22 #include "CAFAna/Cuts/TimingCuts.h"
23 #include "CAFAna/Vars/Vars.h"
26 
27 using namespace ana;
28 
29 // Put a "NOvA Preliminary" tag in the corner
31 {
32  TLatex* prelim = new TLatex(.9, .95, "NO#nuA Preliminary");
33  prelim->SetTextColor(kBlue);
34  prelim->SetNDC();
35  prelim->SetTextSize(2/30.);
36  prelim->SetTextAlign(32);
37  prelim->Draw();
38 }
39 
40 void POTStamp()
41 {
42  TLatex* prelim = new TLatex(.125, .85, "2.74#times 10^{20} POT-equiv.");
43  //prelim->SetTextColor(k);
44  prelim->SetNDC();
45  prelim->SetTextSize(0.04);
46  prelim->SetTextFont(42);
47  prelim->SetTextAlign(12);
48  prelim->Draw();
49 }
50 
52 {
53  TLatex* tag = new TLatex(.1, .95, (pid).c_str());
54  tag->SetTextColor(kGray+1);
55  tag->SetNDC();
56  tag->SetTextSize(0.04);
57  tag->SetTextFont(42);
58  tag->SetTextAlign(12);
59  tag->Draw();
60 }
61 
62 void CenterTitles(TH1* histo)
63 {
64  histo->GetXaxis()->CenterTitle();
65  histo->GetYaxis()->CenterTitle();
66  histo->GetZaxis()->CenterTitle();
67 }
68 
69 void Legend()
70 {
71  TLegend* leg = new TLegend(0.6,0.6,0.85,0.85,NULL,"brNDC");
72  leg->SetFillStyle(0);
73  TH1* dummy = new TH1F("", "", 1, 0, 1);
74  dummy->SetMarkerStyle(kFullCircle);
75 
76  leg->SetTextSize(0.040);
77  leg->AddEntry(dummy->Clone(), "FD data", "lep");
78  dummy->SetLineColor(kTotalMCColor);
79  leg->AddEntry(dummy->Clone(), "Signal prediction", "l");
80  dummy->SetLineColor(kNCBackgroundColor);
81  leg->AddEntry(dummy->Clone(), "Background", "l");
82  leg->Draw();
83 }
84 
86 {
87  TLegend* leg = new TLegend(0.61,0.625,0.86,0.875,NULL,"brNDC");
88  leg->SetFillStyle(0);
89  TH1* dummy = new TH1F("", "", 1, 0, 1);
90  TArrow* arr = new TArrow(2.355,1.920,2.355,1.755,0.025,"|>");
91  arr->SetLineWidth(3);
92  arr->SetFillStyle(1001);
93  arr->SetFillColor(1);
94 
95  leg->SetTextSize(0.040);
96  leg->AddEntry("NULL", "FD data", "");
97  dummy->SetLineColor(kRed);
98  leg->AddEntry(dummy->Clone(), "Signal prediction", "l");
99  dummy->SetLineColor(kNCBackgroundColor);
100  leg->AddEntry(dummy->Clone(), "Background", "l");
101  leg->Draw();
102  arr->Draw();
103 }
104 
105 
106 
107 void plot_fd_datamc(bool makeall = true)
108 {
109  std::cout << "Let's go...." << std::endl;
110 
111  std::string fname = "nue_fa_fd_data_mc_comparisons.root";
112  std::string fname2 = "/nova/ana/nu_e_ana/FirstAna/nue_first_ana_sens_newposc3.root";
113  const double kPOT = 3.45e20;
114 
115  // Use Default Oscillations for prediction
117 
118  // Prediction and cosmic spectra from extrapolation file
119  std::unique_ptr<IPrediction> predLEM = LoadFromFile<IPrediction>(fname2, "predictionLEM");
120  std::unique_ptr<IPrediction> predLID = LoadFromFile<IPrediction>(fname2, "predictionLID");
121 
122  std::unique_ptr<Spectrum> cosmicLEM = LoadFromFile<Spectrum>(fname2, "fdCosmicLEMScale");
123  std::unique_ptr<Spectrum> cosmicLID = LoadFromFile<Spectrum>(fname2, "fdCosmicLIDScale");
124 
125  // Hard-coded arrow values for calE for pictorial representation only.
126  // So that values that overlap can be seen on the plot.
127  TArrow* arrLID[6];
128  TArrow* arrLEM[5];
129  TArrow* arrLEM2[5];
130  TArrow* arrLEM3[5];
131 
132  double caleLEM[5] = {1.67,2.56,1.38,1.49,2.01};
133  double caleLID[6] = {2.210,2.05,2.10,1.98,1.79,2.230};
134 
135  for(int e = 0; e < 5; ++e){
136  arrLEM2[e] = new TArrow(caleLEM[e],0.005,caleLEM[e],0.06,0.025,"<|");
137  arrLEM[e] = new TArrow(caleLEM[e],0.005,caleLEM[e],0.3,0.025,"<|");
138  arrLEM3[e] = new TArrow(caleLEM[e],0.005,caleLEM[e],0.4,0.025,"<|");
139 
140  // for drawing LEM arrows twice
141  // One set just dashed lines no arrow
142  // Second set just arrow-head with black, no-fill arow-head
143  arrLEM[e]->SetFillColor(1);
144  arrLEM[e]->SetFillStyle(0);
145  arrLEM[e]->SetLineStyle(2);
146  arrLEM[e]->SetLineWidth(3);
147 
148  arrLEM2[e]->SetFillColor(1);
149  arrLEM2[e]->SetFillStyle(0);
150  arrLEM2[e]->SetLineWidth(3);
151 
152  // For drawing LEM arrows only, solid black lines, filled arrow-heads
153  // For LEM-only, must also draw arrLID as they are included in LEM
154  // sample
155  arrLEM3[e]->SetFillColor(1);
156  arrLEM3[e]->SetFillStyle(1001);
157  arrLEM3[e]->SetLineWidth(3);
158 
159  }
160 
161  for(int e = 0; e < 6; ++e){
162  arrLID[e] = new TArrow(caleLID[e],0.005,caleLID[e],0.4,0.025,"<|");
163 
164  arrLID[e]->SetFillColor(1);
165  arrLID[e]->SetFillStyle(1001);
166  arrLID[e]->SetLineWidth(3);
167  }
168 
169 
170  if(makeall){
171  struct PlotDef
172  {
173  std::string suffix;
174  std::string xlabel;
175  std::string ylabel;
176  };
177 
178  const int kNumPlots = 32; ///33;
179 
181  {
182  {"nhits","Number of Hits in Slice","Events / bin"},
183  {"nplanes","Number of Planes in Slice","Events / bin"},
184  //{"calE","Calorimetric Energy (GeV)","Events / 0.25 GeV"},
185 
186  {"vtxX","Vertex X Position (cm)","Events / bin"},
187  {"vtxY","Vertex Y Position (cm)","Events / bin"},
188  {"vtxZ","Vertex Z Position (cm)","Events / bin"},
189  {"lid","LID","Events / bin"},
190  {"ptp","p_{T}/p","Events / bin"},
191  {"costheta","cos#theta_{beam}","Events / bin"},
192  {"hits_per_plane","Number of Hits per Plane","Events / bin"},
193  {"pi0mass","#pi^{0} mass (GeV)", "Events / 0.1 GeV"},
194  {"mulll","#mu Longitudinal Log-Likelihood","Events / bin"},
195  {"elll","e Longitudinal Log-Likelihood","Events / bin"},
196  {"emulll","e/#mu Longitudinal Log-Likelihood","Events / bin"},
197  {"eglll","e/#gamma Longitudinal Log-Likelihood","Events / bin"},
198  {"eplll","e/p Longitudinal Log-Likelihood","Events / bin"},
199  {"enlll","e/n Longitudinal Log-Likelihood","Events / bin"},
200  {"epilll","e/#pi Longitudinal Log-Likelihood","Events / bin"},
201  {"epi0lll","e/#pi^{0} Longitudinal Log-Likelihood","Events / bin"},
202  {"mullt","#mu Transverse Log-Likelihood","Events / bin"},
203  {"ellt","e Transverse Log-Likelihood","Events / bin"},
204  {"emullt","e/#mu Transverse Log-Likelihood","Events / bin"},
205  {"egllt","e/#gamma Transverse Log-Likelihood","Events / bin"},
206  {"epllt","e/p Transverse Log-Likelihood","Events / bin"},
207  {"enllt","e/n Transverse Log-Likelihood","Events / bin"},
208  {"epillt","e/#pi Transverse Log-Likelihood","Events / bin"},
209  {"epi0llt","e/#pi^{0} Transverse Log-Likelihood","Events / bin"},
210  {"lem","LEM","Events / bin"},
211  {"lem_pidexp","Fraction of Signal Matches (LEM)","Events / bin"},
212  {"lem_meany","Mean Hadronic y (LEM)","Events / bin"},
213  {"lem_qfrac","Mean Fraction of Charge Matched (LEM)","Events / bin"},
214  {"lem_ediff","Signal - Background Match Potential (LEM)","Events / bin"},
215  {"lem_enrich","Fraction of Enriched Matches (LEM)","Events / bin"}
216  };
217 
218  const int kNumSels = 2;
219  const std::string selNames[kNumSels] = {"LID", "LEM"};
220 
221  // Spectra, data values and predictions for many variables
222  // not just calE that one gets from Erika's extrapolation file
223  Spectrum* dataSpects[kNumSels][kNumPlots];
224  Spectrum* fdCosmic[kNumSels][kNumPlots];
225 
226  IPrediction* predict[kNumSels][kNumPlots];
227 
228  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
229  const char* selName = selNames[selIdx].c_str();
230  for(int varIdx = 0; varIdx < kNumPlots; ++varIdx){
231  const char* name = plots[varIdx].suffix.c_str();
232  dataSpects[selIdx][varIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/data_spect_%s", selName, name).Data()).release();
233  fdCosmic[selIdx][varIdx] = LoadFromFile<Spectrum>(fname, TString::Format("fdCosmic%sScale/cosmic_spect_%s", selName, name).Data()).release();
234  predict[selIdx][varIdx] = LoadFromFile<IPrediction>(fname, TString::Format("prediction%s/prediction_%s", selName, name).Data()).release();
235  } //varIdx
236  }// selIdx
237 
238  // POT of the selected FD Data CAFs for the 11 events
239  const double kDPOT= 3.1844e16;
240 
241  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
242  for(int varIdx = 0; varIdx < kNumPlots; ++varIdx){
243  const std::string name = plots[varIdx].suffix;
244  const std::string pid = selNames[selIdx];
245 
246  Spectrum* data = dataSpects[selIdx][varIdx];
247  Spectrum* cosmic = fdCosmic[selIdx][varIdx];
248  IPrediction* pred = predict[selIdx][varIdx];
249 
250  TH1* hBkgd = cosmic->ToTH1(kPOT);
256  hNC->Add(hNueBeam); // All background (but cosmics) drawn together
257  hNC->Add(hNumu); // All background (but cosmics) drawn together
258  hNC->Add(hNutau);
259 
260 
261  hMC->SetLineColor(kTotalMCColor);
262  hBkgd->SetLineColor(kNCBackgroundColor);
263  // Add all backgrounds to cosmic spectra for total background
264  hBkgd->Add(hNC);
265 
266  TCanvas *c = new TCanvas("c",(name).c_str());
267 
268  TH1* hData = data->ToTH1(kDPOT);
269  hData->SetMarkerStyle(kFullCircle);
270  hData->Draw("ep");
271  hMC->Draw("hist same");
272  hBkgd->Draw("hist same");
273  hData->GetXaxis()->SetTitle((plots[varIdx].xlabel).c_str());
274  hData->GetYaxis()->SetTitle((plots[varIdx].ylabel).c_str());
275  CenterTitles(hData);
276 
277  POTStamp();
278 
279  Legend();
280 
281  std::string out_name = "data_mc_" + pid + "_" + name;
282  // Save plots
283  c->SaveAs(("plots/"+out_name+".pdf").c_str());
284  c->SaveAs(("plots/"+out_name+".eps").c_str());
285  c->SaveAs(("plots/"+out_name+".png").c_str());
286  c->Close();
287  }
288  }
289 
290  } // makeall
291 
292  // LEM/LID CalE on same plot
293 
294 
295  TH1* hBkgdLEM = cosmicLEM->ToTH1(kPOT);
296  TH1* hBkgdLID = cosmicLID->ToTH1(kPOT);
297 
298  TH1* hMCLEM = predLEM->PredictComponent(calc, ana::Flavors::kNuMuToNuE, ana::Current::kCC, ana::Sign::kBoth).ToTH1(kPOT);
299  TH1* hNutauLEM = predLEM->PredictComponent(calc, ana::Flavors::kAllNuTau, ana::Current::kCC, ana::Sign::kBoth).ToTH1(kPOT);
300  TH1* hNueBeamLEM = predLEM->PredictComponent(calc, ana::Flavors::kNuEToNuE, ana::Current::kCC, ana::Sign::kBoth).ToTH1(kPOT);
301  TH1* hNumuLEM = predLEM->PredictComponent(calc, ana::Flavors::kAllNuMu, ana::Current::kCC, ana::Sign::kBoth).ToTH1(kPOT);
302  TH1* hNCLEM = predLEM->PredictComponent(calc, ana::Flavors::kAll, ana::Current::kNC, ana::Sign::kBoth).ToTH1(kPOT);
303  hNCLEM->Add(hNueBeamLEM); // All background (but cosmics) drawn together
304  hNCLEM->Add(hNumuLEM); // All background (but cosmics) drawn together
305  hNCLEM->Add(hNutauLEM);
306 
307 
308  TH1* hMCLID = predLID->PredictComponent(calc, ana::Flavors::kNuMuToNuE, ana::Current::kCC, ana::Sign::kBoth).ToTH1(kPOT);
309  TH1* hNutauLID = predLID->PredictComponent(calc, ana::Flavors::kAllNuTau, ana::Current::kCC, ana::Sign::kBoth).ToTH1(kPOT);
310  TH1* hNueBeamLID = predLID->PredictComponent(calc, ana::Flavors::kNuEToNuE, ana::Current::kCC, ana::Sign::kBoth).ToTH1(kPOT);
311  TH1* hNumuLID = predLID->PredictComponent(calc, ana::Flavors::kAllNuMu, ana::Current::kCC, ana::Sign::kBoth).ToTH1(kPOT);
312  TH1* hNCLID = predLID->PredictComponent(calc, ana::Flavors::kAll, ana::Current::kNC, ana::Sign::kBoth).ToTH1(kPOT);
313  hNCLID->Add(hNueBeamLID); // All background (but cosmics) drawn together
314  hNCLID->Add(hNumuLID); // All background (but cosmics) drawn together
315  hNCLID->Add(hNutauLID);
316 
317  TCanvas *cBoth = new TCanvas("c","both");
318 
319  hMCLEM->SetTitle("");
320  hMCLEM->SetLineColor(kTotalMCColor);
321  hMCLEM->SetLineStyle(2);
322  hMCLEM->GetYaxis()->SetRangeUser(0.,2.);
323  hMCLEM->GetXaxis()->SetRangeUser(1.,3.);
324  hMCLEM->Draw("hist");
325 
326  hBkgdLEM->SetTitle("");
327  hBkgdLEM->Add(hNCLEM);
328  hBkgdLEM->SetLineColor(kNCBackgroundColor);
329  hBkgdLEM->SetLineStyle(2);
330  hBkgdLEM->Draw("hist same");
331 
332  hMCLID->SetTitle("");
333  hMCLID->SetLineColor(kTotalMCColor);
334  hMCLID->Draw("hist same");
335 
336  hBkgdLID->SetTitle("");
337  hBkgdLID->Add(hNCLID);
338  hBkgdLID->SetLineColor(kNCBackgroundColor);
339  hBkgdLID->Draw("hist same");
340 
341  for(int e = 0; e < 5; ++e){
342  arrLEM[e]->Draw();
343  arrLEM2[e]->Draw();
344  }
345 
346  for(int e = 0; e < 6; ++e){
347  arrLID[e]->Draw();
348  }
349 
350  CenterTitles(hMCLEM);
351 
352  hMCLEM->GetYaxis()->SetTitle("Events / 0.25 GeV");
353  hMCLEM->GetXaxis()->SetTitle("Calorimetric energy (GeV)");
354 
355  POTStamp();
356  //Preliminary();
357  LegendSpecial();
358 
359  cBoth->SaveAs("plots/data_mc_both_calE.png");
360  cBoth->SaveAs("plots/data_mc_both_calE.eps");
361  cBoth->SaveAs("plots/data_mc_both_calE.pdf");
362  cBoth->Close();
363 
364  // LID-only
365  TCanvas *cLID = new TCanvas("clid","lid");
366  hMCLID->GetYaxis()->SetRangeUser(0.,2.);
367  hMCLID->GetXaxis()->SetRangeUser(1.,3.);
368 
369  hMCLID->Draw();
370  hBkgdLID->Draw("hist same");
371  for(int e = 0; e < 6; ++e){
372  arrLID[e]->Draw();
373  }
374  CenterTitles(hMCLID);
375  hMCLID->GetYaxis()->SetTitle("Events / 0.25 GeV");
376  hMCLID->GetXaxis()->SetTitle("Calorimetric energy (GeV)");
377  POTStamp();
378  //Preliminary();
379  LegendSpecial();
380  PIDTag("LID");
381  cLID->SaveAs("plots/data_mc_LID_calE.png");
382  cLID->SaveAs("plots/data_mc_LID_calE.eps");
383  cLID->SaveAs("plots/data_mc_LID_calE.pdf");
384  cLID->Close();
385 
386  // LEM-only
387  TCanvas *cLEM = new TCanvas("clem","lem");
388  hMCLEM->SetLineStyle(1);
389  hBkgdLEM->SetLineStyle(1);
390  hMCLEM->Draw();
391  hBkgdLEM->Draw("hist same");
392  for(int e = 0; e < 5; ++e){
393  arrLEM3[e]->Draw();
394  }
395  for(int e = 0; e < 6; ++e){
396  arrLID[e]->Draw();
397  }
398 
399  CenterTitles(hMCLEM);
400  hMCLEM->GetYaxis()->SetTitle("Events / 0.25 GeV");
401  hMCLEM->GetXaxis()->SetTitle("Calorimetric energy (GeV)");
402  POTStamp();
403  //Preliminary();
404  LegendSpecial();
405  PIDTag("LEM");
406  cLEM->SaveAs("plots/data_mc_LEM_calE.png");
407  cLEM->SaveAs("plots/data_mc_LEM_calE.eps");
408  cLEM->SaveAs("plots/data_mc_LEM_calE.pdf");
409  cLEM->Close();
410 
411 }
void Legend()
const XML_Char * name
Definition: expat.h:151
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
(&#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
void POTStamp()
(&#39;beam &#39;)
Definition: IPrediction.h:15
void PIDTag(std::string pid)
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1128
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
const XML_Char const XML_Char * data
Definition: expat.h:268
Charged-current interactions.
Definition: IPrediction.h:39
General interface to any calculator that lets you set the parameters.
Sum up livetimes from individual cosmic triggers.
const int kNumPlots
Definition: GetSpectra.h:1
TH2D * histo
void plot_fd_datamc(bool makeall=true)
const int kNumSels
Definition: vars.h:43
TLatex * prelim
Definition: Xsec_final.C:133
OStream cout
Definition: OStream.cxx:6
const std::vector< PlotDef > plots
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Neutral-current interactions.
Definition: IPrediction.h:40
std::string suffix
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
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
All neutrinos, any flavor.
Definition: IPrediction.h:26
const std::string selNames[kNumSels]
Definition: vars.h:46
void Preliminary()
Float_t e
Definition: plot.C:35
const Color_t kNCBackgroundColor
Definition: Style.h:22
void LegendSpecial()