plotDataPred_SplitCanvas.C
Go to the documentation of this file.
1 #include "includes.h"
2 #include "numu_tools.h"
3 #include "varsandcuts.h"
4 
5 using namespace ana;
6 
7 
9 
10 
11  bool preliminary = false;
12 
15  auto calcNoOsc = osc.Copy();
16 
17  std::string CalcName = "3ABestFit";
18  // 3a best fit
19  calc->SetTh23(asin(sqrt(0.559455)));
20  calc->SetDmsq32(0.00244);
21 
22  // sample
23  std::string SampleName = "9e20";
24  double pot = kAna2017POT;
25  double livetime = kAna2017Livetime;
26 
27  //in-outputs
28  TString outDir = "";
29  std::string inDirData = "/nova/ana/nu_mu_ana/Ana2017/Data";
30  std::string inDir = "/nova/ana/users/dmendez/Numu/Sensitivities2017/WithSysts/SABestFit/";
31  std::string inName = "Prediction_SABestFit_WithSysts_Extrap_OptimisedAna__";
32 
33  // data
34  std::vector<Spectrum*> sData;
35  std::vector<TH1*> hData;
36  std::vector<TH1*> hDataUnscaled;
37  std::cout << "\nLoading data" << std::endl;
38  TFile* fDataQ0 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant0.root").c_str() );
39  TFile* fDataQ1 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant1.root").c_str() );
40  TFile* fDataQ2 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant2.root").c_str() );
41  TFile* fDataQ3 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant3.root").c_str() );
42  TFile* fDataQ4 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant4.root").c_str() );
43  sData.push_back(LoadFrom< Spectrum >(fDataQ0, "fd_data").release());
44  sData.push_back(LoadFrom< Spectrum >(fDataQ1, "fd_data").release());
45  sData.push_back(LoadFrom< Spectrum >(fDataQ2, "fd_data").release());
46  sData.push_back(LoadFrom< Spectrum >(fDataQ3, "fd_data").release());
47  sData.push_back(LoadFrom< Spectrum >(fDataQ4, "fd_data").release());
48 
49  for(int quantId = 0; quantId<nQuantPlus; quantId++){
50  hData.push_back(sData[quantId]->ToTH1(pot));
51  hDataUnscaled.push_back(sData[quantId]->ToTH1(pot));
52  }
53 
54  fDataQ1 -> Close();
55  fDataQ2 -> Close();
56  fDataQ3 -> Close();
57  fDataQ4 -> Close();
58  delete fDataQ1;
59  delete fDataQ2;
60  delete fDataQ3;
61  delete fDataQ4;
62 
63 
64 
65  // edited plotWithSystErrorBand
66  // eats the top of user range
67  // and sets yaxis to events/0.1 gev
68  std::vector<float> maxGraph;
69  maxGraph.push_back(13.0);
70  maxGraph.push_back(6.0);
71  maxGraph.push_back(6.0);
72  maxGraph.push_back(6.0);
73  maxGraph.push_back(6.0);
74 
75 
76 
77  //cosmics
78  std::string CosmicsSpecFile;
79  std::string CosmicsHistName;
80  std::string cosmicsName = "NoCosmics";
81  std::cout << "\nCosmics" << std::endl;
82  CosmicsSpecFile = "/nova/ana/nu_mu_ana/Ana2017/Cosmics/cosmicsScaled_NewBinning_" + SampleName + ".root";
83  CosmicsHistName = "cosmics3A";
84 
85  std::cout << "\nLoading cosmic distributions" << std::endl;
86  TFile inCosmicsFile(CosmicsSpecFile.c_str());
87  std::vector<Spectrum> sCosmics;
88  std::vector<TH1*> hCosmics;
89  hCosmics.push_back((TH1*)inCosmicsFile.Get("q1all"));
90  hCosmics.push_back((TH1*)inCosmicsFile.Get("q1all"));
91  hCosmics.push_back((TH1*)inCosmicsFile.Get("q2all"));
92  hCosmics.push_back((TH1*)inCosmicsFile.Get("q3all"));
93  hCosmics.push_back((TH1*)inCosmicsFile.Get("q4all"));
94  for(int quantId = 2; quantId < nQuantPlus; quantId++){
95  hCosmics[0]->Add(hCosmics[quantId]);
96  }
97  for(int quantId = 0; quantId<nQuantPlus; quantId++){
98  Spectrum cosmics(hCosmics[quantId], pot, livetime);
99  sCosmics.push_back(cosmics);
100  }
101 
102 
103 
104  std::cout << "\nDefining predictions and systematics" << std::endl;
105  std::vector<PredictionInterp*> predictionVec;
107  ENumu2017ExtrapType extrap = kNuMu;
108  PredictionSystNumu2017 predQ1(extrap, calc, 1);
109  PredictionSystNumu2017 predQ2(extrap, calc, 2);
110  PredictionSystNumu2017 predQ3(extrap, calc, 3);
111  PredictionSystNumu2017 predQ4(extrap, calc, 4);
112  predictionVec.push_back(&predQ1);
113  predictionVec.push_back(&predQ2);
114  predictionVec.push_back(&predQ3);
115  predictionVec.push_back(&predQ4);
116 
117 
118  /*
119  //predictions
120  std::vector <string> epoch = {"full"};
121  for(int quantId = 0; quantId < 4; quantId++){
122  std::cout << "Quantile " << quantId + 1 << std::endl;
123  std::string input = inDir + inName + Form("Quant%d_", quantId+1) + "full.root";
124  TFile* file = new TFile(input.c_str());
125  std::cout << "\nPrediction vectors" << std::endl;
126  predictionVec.push_back(LoadFrom<PredictionInterp>(file, "prediction")).release();
127  file->Close();
128  std::cout << "Read and closed input files" << std::endl;
129  }//loop quant
130  */
131 
132 
133  std::vector<TH1*> hNC;
134  std::vector<TH1*> hCC;
135  std::vector<TH1*> hNuMu;
136  std::vector<TH1*> hNoOsc;
137  std::vector<TH1*> hPred;
138  std::vector<TH1*> hBkg;
139  std::vector<TH1*> hAllBkg;
140 
141  hNC.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
142  hCC.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
143  hNuMu.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
144  hBkg.push_back(predictionVec[0]->Predict(calc).ToTH1(kAna2017POT, kGray));
145  hAllBkg.push_back(predictionVec[0]->Predict(calc).ToTH1(kAna2017POT, kGray));
146  hPred.push_back(predictionVec[0]->Predict(calc).ToTH1(kAna2017POT, kGray));
147  hNoOsc.push_back(predictionVec[0]->Predict(calcNoOsc).ToTH1(kAna2017POT, kGray));
148 
149  for(int quantId = 0; quantId < nQuant; quantId++){
150  hNC.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
151  hCC.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
152  hNuMu.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
153  hBkg.push_back(predictionVec[quantId]->Predict(calc).ToTH1(kAna2017POT, kGray));
154  hAllBkg.push_back(predictionVec[quantId]->Predict(calc).ToTH1(kAna2017POT, kGray));
155  hPred.push_back(predictionVec[quantId]->Predict(calc).ToTH1(kAna2017POT, kGray));
156  hNoOsc.push_back(predictionVec[quantId]->Predict(calcNoOsc).ToTH1(kAna2017POT, kGray));
157  if(quantId > 0){ // already added the first quantile
158  hNC[0] ->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
159  hCC[0] ->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
160  hNuMu[0] ->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
161  hBkg[0] ->Add(predictionVec[quantId]->Predict(calc).ToTH1(kAna2017POT, kGray));
162  hAllBkg[0]->Add(predictionVec[quantId]->Predict(calc).ToTH1(kAna2017POT, kGray));
163  hPred[0] ->Add(predictionVec[quantId]->Predict(calc).ToTH1(kAna2017POT, kGray));
164  hNoOsc[0] ->Add(predictionVec[quantId]->Predict(calcNoOsc).ToTH1(kAna2017POT, kGray));
165  }
166  }
167  for(int quantId = 0; quantId < nQuantPlus; quantId++){
168  hBkg[quantId]->Add(hNuMu[quantId], -1);
169  hAllBkg[quantId]->Add(hNuMu[quantId], -1);
170  hAllBkg[quantId]->Add(hCosmics[quantId]);
171  }
172 
173 
174 
175  // scale everything, everything!
176  std::cout << "\nScaling histograms to 0.1 GeV" << std::endl;
177  for(int quantId = 0; quantId < nQuantPlus; quantId++){
178  hData[quantId]->Scale(0.1, "width");
179  hCosmics[quantId]->Scale(0.1, "width");
180  hNC[quantId]->Scale(0.1, "width");
181  hCC[quantId]->Scale(0.1, "width");
182  hBkg[quantId]->Scale(0.1, "width");
183  hAllBkg[quantId]->Scale(0.1, "width");
184  hNuMu[quantId]->Scale(0.1, "width");
185  hPred[quantId]->Scale(0.1, "width");
186  hNoOsc[quantId]->Scale(0.1, "width");
187  }
188 
189 
190  // data graphs
191  std::cout << "\nDefine data graphs\n" << std::endl;
192  std::vector<TGraph*> grData;
193  for(int quantId = 0; quantId < nQuantPlus; quantId++){
194  grData.push_back(graphAsymmError(hData[quantId], hDataUnscaled[quantId]));
195  }
196 
197 
198  // predictions 3a best fit
199  std::cout << "3A best fit prediction" << std::endl;
200  std::vector<TH1*> hPrediction;
201  hPrediction.push_back(predictionVec[0]->Predict(calc).ToTH1(pot));
202  for(int quantId = 0; quantId < nQuant; quantId++){
203  hPrediction.push_back(predictionVec[quantId]->Predict(calc).ToTH1(pot));
204  }
205  hPrediction[0]->Add(hPrediction[2]);
206  hPrediction[0]->Add(hPrediction[3]);
207  hPrediction[0]->Add(hPrediction[4]);
208  for(int quantId = 0; quantId < nQuantPlus; quantId++){
209  hPrediction[quantId]->Scale(0.1, "width");
210  hPrediction[quantId]->Add(hCosmics[quantId]); // add in nominal cosmics
211  }
212  std::vector<Spectrum> sPrediction;
213  for(int quantId = 0; quantId < nQuantPlus; quantId++){
214  Spectrum prediction(hPrediction[quantId], pot, livetime);
215  sPrediction.push_back(prediction);
216  }
217 
218 
219 
220  TCanvas *canvas;
221  TPad *pad1, *pad2, *pad3, *pad4;
222  SplitCanvasQuant(canvas, pad1, pad2, pad3, pad4);
223 
224  std::cout << "\nAbout to plot\n" << std::endl;
225  TString outName = outDir + "datapred_bestfit_4pads";
226  for(int quantId = 3; quantId <=4; quantId++){
227 
228  std::cout << "quantile " << quantId << std::endl;
229  hPrediction[quantId]->SetLineColor(kRed);
230  hPrediction[quantId]->SetLineWidth(3);
231  hCosmics[quantId]->SetLineColor(kGray+2);
232  hCosmics[quantId]->SetLineWidth(2);
233  hBkg[quantId]->SetLineColor(kBlue);
234  hBkg[quantId]->SetLineStyle(kDashed);
235  hBkg[quantId]->SetLineWidth(2);
236  hAllBkg[quantId]->SetLineColor(kBlue);
237  hAllBkg[quantId]->SetLineStyle(kDashed);
238  hAllBkg[quantId]->SetLineWidth(2);
239 
240  // spectrums for syst band best fit
241  std::vector<Spectrum> ups, dns;
242  for(const ISyst* syst: systs){
243  SystShifts shifts;
244  shifts.SetShift(syst, +1);
245  TH1 * htempup = predictionVec[quantId-1]->PredictSyst(calc, shifts).ToTH1(pot);
246  htempup->Scale(0.1, "width");
247  htempup->Add(hCosmics[quantId]);
248  Spectrum specup(htempup, pot, livetime);
249  ups.push_back(specup);
250  shifts.SetShift(syst, -1);
251  TH1 * htempdn = predictionVec[quantId-1]->PredictSyst(calc, shifts).ToTH1(pot);
252  htempdn->Scale(0.1, "width");
253  htempdn->Add(hCosmics[quantId]);
254  Spectrum specdn(htempdn, pot, livetime);
255  dns.push_back(specdn);
256  }
257 
258  if(quantId==1) pad1->cd();
259  if(quantId==2) pad2->cd();
260  if(quantId==3) pad3->cd();
261  if(quantId==4) pad4->cd();
262  hCosmics[quantId]->GetXaxis()->SetLabelSize(0.10);
263  hCosmics[quantId]->GetXaxis()->SetTitleSize(0.10);
264  PlotWithSystErrorBand_Quant(quantId, sPrediction[quantId], ups, dns, pot, -1, -1, maxGraph[quantId], true);
265  hCosmics[quantId]->Draw("hist same");
266  hBkg[quantId]->Draw("hist same");
267  grData[quantId]->Draw("p0 same");
268  drawQuantLabel(quantId);
269  } // for quantiles
270 
271 
272 
273  for(int quantId = 1; quantId <=2; quantId++){
274 
275  std::cout << "quantile " << quantId << std::endl;
276  hPrediction[quantId]->SetLineColor(kRed);
277  hPrediction[quantId]->SetLineWidth(3);
278  hCosmics[quantId]->SetLineColor(kGray+2);
279  hCosmics[quantId]->SetLineWidth(2);
280  hBkg[quantId]->SetLineColor(kBlue);
281  hBkg[quantId]->SetLineStyle(kDashed);
282  hBkg[quantId]->SetLineWidth(2);
283  hAllBkg[quantId]->SetLineColor(kBlue);
284  hAllBkg[quantId]->SetLineStyle(kDashed);
285  hAllBkg[quantId]->SetLineWidth(2);
286 
287  // spectrums for syst band best fit
288  std::vector<Spectrum> ups, dns;
289  for(const ISyst* syst: systs){
290  SystShifts shifts;
291  shifts.SetShift(syst, +1);
292  if(quantId==0){
293  TH1 * htempup = predictionVec[0]->PredictSyst(calc, shifts).ToTH1(pot);
294  htempup->Add(predictionVec[1]->PredictSyst(calc, shifts).ToTH1(pot));
295  htempup->Add(predictionVec[2]->PredictSyst(calc, shifts).ToTH1(pot));
296  htempup->Add(predictionVec[3]->PredictSyst(calc, shifts).ToTH1(pot));
297  htempup->Scale(0.1, "width");
298  htempup->Add(hCosmics[0]);
299  Spectrum specup(htempup, pot, livetime);
300  ups.push_back(specup);
301  shifts.SetShift(syst, -1);
302  TH1 * htempdn = predictionVec[0]->PredictSyst(calc, shifts).ToTH1(pot);
303  htempdn->Add(predictionVec[1]->PredictSyst(calc, shifts).ToTH1(pot));
304  htempdn->Add(predictionVec[2]->PredictSyst(calc, shifts).ToTH1(pot));
305  htempdn->Add(predictionVec[3]->PredictSyst(calc, shifts).ToTH1(pot));
306  htempdn->Scale(0.1, "width");
307  htempdn->Add(hCosmics[0]);
308  Spectrum specdn(htempdn, pot, livetime);
309  dns.push_back(specdn);
310  }
311  else{
312  TH1 * htempup = predictionVec[quantId-1]->PredictSyst(calc, shifts).ToTH1(pot);
313  htempup->Scale(0.1, "width");
314  htempup->Add(hCosmics[quantId]);
315  Spectrum specup(htempup, pot, livetime);
316  ups.push_back(specup);
317  shifts.SetShift(syst, -1);
318  TH1 * htempdn = predictionVec[quantId-1]->PredictSyst(calc, shifts).ToTH1(pot);
319  htempdn->Scale(0.1, "width");
320  htempdn->Add(hCosmics[quantId]);
321  Spectrum specdn(htempdn, pot, livetime);
322  dns.push_back(specdn);
323  }
324  }
325 
326 
327  if(quantId==1) pad1->cd();
328  if(quantId==2) pad2->cd();
329  if(quantId==3) pad3->cd();
330  if(quantId==4) pad4->cd();
331  hCosmics[quantId]->GetXaxis()->SetLabelSize(0.10);
332  hCosmics[quantId]->GetXaxis()->SetTitleSize(0.10);
333  PlotWithSystErrorBand_Quant(quantId, sPrediction[quantId], ups, dns, pot, -1, -1, maxGraph[quantId], true);
334  hCosmics[quantId]->Draw("hist same");
335  hBkg[quantId]->Draw("hist same");
336  grData[quantId]->Draw("p0 same");
337  drawQuantLabel(quantId);
338  if(quantId==1){
339  pad1->cd();
340  TLegend *legend = new TLegend(0.09,0.73,0.27,0.87,"","NDC");
341  legend->AddEntry(grData[quantId], "Data","lep");
342  legend->AddEntry(hPrediction[quantId], "Prediction","l");
343  TLegendEntry *entry=legend->AddEntry("error","1-#sigma syst. range","f");
344  entry->SetFillStyle(1001);
345  entry->SetFillColor(kRed-10);
346  legend->AddEntry(hBkg[quantId], "Beam bkg.","l");
347  legend->AddEntry(hCosmics[quantId], "Cosmic bkg.","l");
348  legend->SetTextSize(0.025); //no border for legend
349  legend->SetBorderSize(0); //no border for legend
350  legend->SetFillStyle(0); //fill colour is white
351  legend->Draw();
352  gPad->Update();
353  }
354  } // for quantiles
355 
356 
357 
358  canvas->cd();
359  TLatex* xTitle = new TLatex(0.5, 0.03, "Reconstructed Energy (GeV)");
360  xTitle->SetNDC();
361  xTitle->SetTextAlign(22);
362  xTitle->Draw();
363  canvas->cd();
364  TLatex* yTitle = new TLatex(0.02, 0.5, "Events / 0.1 GeV");
365  yTitle->SetNDC();
366  yTitle->SetTextAlign(22);
367  yTitle->SetTextAngle(90.);
368  yTitle->Draw();
369 
370 
371  canvas->SaveAs(outName + "_testpads.png");
372  canvas->SaveAs(outName + "_testpads.eps");
373  canvas->SaveAs(outName + "_testpads.pdf");
374 
375 
376 
377 
378 } // End of function
Pass neutrinos through unchanged.
Definition: IOscCalc.h:45
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
TGraph * graphAsymmError(TH1 *hDataScaled, TH1 *hData)
Definition: tools.h:59
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double kAna2017POT
Definition: Exposures.h:174
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
virtual void SetDmsq32(const T &dmsq32)=0
std::string outDir
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
TString inDir
Charged-current interactions.
Definition: IPrediction.h:39
const double kAna2017Livetime
Definition: Exposures.h:200
Loads shifted spectra from files.
virtual _IOscCalc< T > * Copy() const override
Definition: IOscCalc.h:49
#define pot
const int nQuantPlus
Definition: varsandcuts.h:5
Int_t preliminary
Definition: SimpleIterate.C:63
c1 Close()
Oscillation probability calculators.
Definition: Calcs.h:5
void plotDataPred_SplitCanvas()
OStream cout
Definition: OStream.cxx:6
double livetime
Definition: saveFDMCHists.C:21
TGraphAsymmErrors * PlotWithSystErrorBand_Quant(const int quant, IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float maxy, bool newaxis)
Definition: Plots.cxx:541
const int nQuant
Definition: varsandcuts.h:4
TPaveText * drawQuantLabel(int quantId=0)
Definition: numu_tools.h:4
virtual void SetTh23(const T &th23)=0
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
TPad * pad3
Definition: analysis.C:13
void SplitCanvasQuant(TCanvas *&canvas, TPad *&pad1, TPad *&pad2, TPad *&pad3, TPad *&pad4)
Definition: Plots.cxx:1437
TPad * pad2
Definition: analysis.C:13
std::vector< const ISyst * > getAllAna2017Systs(const EAnaType2017 ana, const bool smallgenie)
All neutrinos, any flavor.
Definition: IPrediction.h:26
enum BeamMode kBlue
def entry(str)
Definition: HTMLTools.py:26
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
TPad * pad1
Definition: analysis.C:13
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string