plotDataPred.C
Go to the documentation of this file.
1 #include "includes.h"
2 #include "tools.h"
3 #include "varsandcuts.h"
4 
5 using namespace ana;
6 
7 
8 void plotDataPred(bool newpred = true, bool plotunosc = true, bool plotchi = true, bool plotshifts = true){
9 
10  bool preliminary = true;
11 
12 
16  //osc::IOscCalcAdjustable* calcNoOsc = DefaultOscCalc();
17  //osc::NoOscillations* calcNoOsc;
19  auto calcNoOsc = osc.Copy();
20 
21  std::string CalcName = "3ABestFit";
22  // 3a best fit
23  calc->SetTh23(asin(sqrt(0.559455)));
24  calc->SetDmsq32(0.00244);
25  // sa best fit
26  calcSA->SetDmsq32(2.6746e-3);
27  calcSA->SetTh23(0.68696);
28  // maximal mixing scenario
29  calcMaxMix->SetTh23(asin(sqrt(0.5))); // sin2theta23 = 0.5
30  calcMaxMix->SetDmsq32(0.00244); // 3a best fit
31 
32 
33  // sample
34  std::string SampleName = "9e20";
35  double pot = kAna2017POT;
36  double livetime = kAna2017Livetime;
37 
38 
39  //inputs
40  TString outDir = "";
41  std::string inDirData = "/nova/ana/nu_mu_ana/Ana2017/Data";
42  std::string inDir = "/nova/ana/nu_mu_ana/Ana2017/Predictions";
43  std::string inName = "Prediction_StatsOnly_3ACut_3AEnergy_OptBinning_";
44 
45  // data
46  std::vector<Spectrum*> sData;
47  std::vector<TH1*> hData;
48  std::vector<TH1*> hDataUnscaled;
49  TFile* fDataQ0 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant0.root").c_str() );
50  TFile* fDataQ1 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant1.root").c_str() );
51  TFile* fDataQ2 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant2.root").c_str() );
52  TFile* fDataQ3 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant3.root").c_str() );
53  TFile* fDataQ4 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant4.root").c_str() );
54 
55  sData.push_back(LoadFrom< Spectrum >(fDataQ0, "fd_data").release());
56  sData.push_back(LoadFrom< Spectrum >(fDataQ1, "fd_data").release());
57  sData.push_back(LoadFrom< Spectrum >(fDataQ2, "fd_data").release());
58  sData.push_back(LoadFrom< Spectrum >(fDataQ3, "fd_data").release());
59  sData.push_back(LoadFrom< Spectrum >(fDataQ4, "fd_data").release());
60 
61  for(int quantId = 0; quantId<nQuantPlus; quantId++){
62  hData.push_back(sData[quantId]->ToTH1(pot));
63  hDataUnscaled.push_back(sData[quantId]->ToTH1(pot));
64  }
65 
66  fDataQ1 -> Close();
67  fDataQ2 -> Close();
68  fDataQ3 -> Close();
69  fDataQ4 -> Close();
70  delete fDataQ1;
71  delete fDataQ2;
72  delete fDataQ3;
73  delete fDataQ4;
74 
75 
76 
77 
78  //cosmics
79  std::string CosmicsSpecFile;
80  std::string CosmicsHistName;
81  std::string cosmicsName = "NoCosmics";
82  std::cout << "\nCosmics" << std::endl;
83  CosmicsSpecFile = "/nova/ana/nu_mu_ana/Ana2017/Cosmics/cosmicsScaled_NewBinning_" + SampleName + ".root";
84  CosmicsHistName = "cosmics3A";
85 
86  std::cout << "\nLoading cosmic distributions" << std::endl;
87  TFile inCosmicsFile(CosmicsSpecFile.c_str());
88 
89  std::vector<Spectrum> sCosmics;
90  std::vector<TH1*> hCosmics;
91  hCosmics.push_back((TH1*)inCosmicsFile.Get("q1all"));
92  hCosmics.push_back((TH1*)inCosmicsFile.Get("q1all"));
93  hCosmics.push_back((TH1*)inCosmicsFile.Get("q2all"));
94  hCosmics.push_back((TH1*)inCosmicsFile.Get("q3all"));
95  hCosmics.push_back((TH1*)inCosmicsFile.Get("q4all"));
96  for(int quantId = 2; quantId < nQuantPlus; quantId++){
97  hCosmics[0]->Add(hCosmics[quantId]);
98  }
99  for(int quantId = 0; quantId<nQuantPlus; quantId++){
100  Spectrum cosmics(hCosmics[quantId], pot, livetime);
101  sCosmics.push_back(cosmics);
102  }
103 
104 
105 
106 
107  std::cout << "\nDefining predictions and systematics" << std::endl;
108  std::vector<PredictionInterp*> predictionVec;
110  ENumu2017ExtrapType extrap = kNuMu;
111 
112  //if(newpred){
113  PredictionSystNumu2017 predQ1(extrap, calc, 1);
114  PredictionSystNumu2017 predQ2(extrap, calc, 2);
115  PredictionSystNumu2017 predQ3(extrap, calc, 3);
116  PredictionSystNumu2017 predQ4(extrap, calc, 4);
117  predictionVec.push_back(&predQ1);
118  predictionVec.push_back(&predQ2);
119  predictionVec.push_back(&predQ3);
120  predictionVec.push_back(&predQ4);
121  //}
122  /*
123  else{
124  systs = getAllDirectNumuSysts2017();
125  for(int quantId = 0; quantId < nQuant; quantId++){
126  std::cout << "Quantile " << quantId + 1 << std::endl;
127  std::string input = inDir + inName + Form("Quant%d_", quantId+1) + "full_" + SampleName + ".root";
128  TFile* file = new TFile(input.c_str());
129  std::cout << "\nPrediction vectors" << std::endl;
130  predictionVec.push_back(LoadFrom<PredictionInterp>(file, "prediction")).release();
131  file->Close();
132  std::cout << "Read and closed input files" << std::endl;
133  }//loop quant
134  }
135  */
136  std::cout << "\nList of systematics:" << std::endl;
137  for (auto & sys:systs) std::cout <<sys->ShortName() << std::endl;
138 
139 
140  TFile* fout = new TFile(outDir + "plotDataPred.root","RECREATE");
141 
142  std::vector<TH1*> hNC;
143  std::vector<TH1*> hCC;
144  std::vector<TH1*> hNuMu;
145  std::vector<TH1*> hNoOsc;
146  std::vector<TH1*> hPred;
147  std::vector<TH1*> hBkg;
148  std::vector<TH1*> hAllBkg;
149 
150  hNC.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
151  hCC.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
152  hNuMu.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
153  hBkg.push_back(predictionVec[0]->Predict(calc).ToTH1(kAna2017POT, kGray));
154  hAllBkg.push_back(predictionVec[0]->Predict(calc).ToTH1(kAna2017POT, kGray));
155  hPred.push_back(predictionVec[0]->Predict(calc).ToTH1(kAna2017POT, kGray));
156  hNoOsc.push_back(predictionVec[0]->Predict(calcNoOsc).ToTH1(kAna2017POT, kGray));
157 
158  for(int quantId = 0; quantId < nQuant; quantId++){
159  hNC.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
160  hCC.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
161  hNuMu.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
162  hBkg.push_back(predictionVec[quantId]->Predict(calc).ToTH1(kAna2017POT, kGray));
163  hAllBkg.push_back(predictionVec[quantId]->Predict(calc).ToTH1(kAna2017POT, kGray));
164  hPred.push_back(predictionVec[quantId]->Predict(calc).ToTH1(kAna2017POT, kGray));
165  hNoOsc.push_back(predictionVec[quantId]->Predict(calcNoOsc).ToTH1(kAna2017POT, kGray));
166  if(quantId > 0){ // already added the first quantile
167  hNC[0] ->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
168  hCC[0] ->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
169  hNuMu[0] ->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
170  hBkg[0] ->Add(predictionVec[quantId]->Predict(calc).ToTH1(kAna2017POT, kGray));
171  hAllBkg[0]->Add(predictionVec[quantId]->Predict(calc).ToTH1(kAna2017POT, kGray));
172  hPred[0] ->Add(predictionVec[quantId]->Predict(calc).ToTH1(kAna2017POT, kGray));
173  hNoOsc[0] ->Add(predictionVec[quantId]->Predict(calcNoOsc).ToTH1(kAna2017POT, kGray));
174  }
175  }
176  for(int quantId = 0; quantId < nQuantPlus; quantId++){
177  hBkg[quantId]->Add(hNuMu[quantId], -1);
178  hAllBkg[quantId]->Add(hNuMu[quantId], -1);
179  hAllBkg[quantId]->Add(hCosmics[quantId]);
180  }
181 
182 
183 
184  // scale everything, everything!
185  std::cout << "\nScaling histograms to 0.1 GeV" << std::endl;
186  for(int quantId = 0; quantId < nQuantPlus; quantId++){
187  hData[quantId]->Scale(0.1, "width");
188  hCosmics[quantId]->Scale(0.1, "width");
189  hNC[quantId]->Scale(0.1, "width");
190  hCC[quantId]->Scale(0.1, "width");
191  hBkg[quantId]->Scale(0.1, "width");
192  hAllBkg[quantId]->Scale(0.1, "width");
193  hNuMu[quantId]->Scale(0.1, "width");
194  hPred[quantId]->Scale(0.1, "width");
195  hNoOsc[quantId]->Scale(0.1, "width");
196  }
197 
198 
199  // data graphs
200  std::cout << "\nDefine data graphs\n" << std::endl;
201  std::vector<TGraph*> grData;
202  for(int quantId = 0; quantId < nQuantPlus; quantId++){
203  grData.push_back(graphAsymmError(hData[quantId], hDataUnscaled[quantId]));
204  }
205 
206  /*
207  // add in the cosmic syst by hand
208  TH1 * htempup = predSACC.Predict(calc).ToTH1(pot);
209  htempup->Add(hCosmics[0],1 + (0.33/2.81));
210  Spectrum specup(htempup, pot, livetime);
211  ups.push_back(specup);
212  TH1 * htempdn = predSACC.Predict(calc).ToTH1(pot);
213  htempdn->Add(hCosmics[0],1 - (0.33/2.81));
214  Spectrum specdn(htempdn, pot, livetime);
215  dns.push_back(specdn);
216  */
217 
218 
219 
220  // predictions 3a best fit
221  std::cout << "3A best fit prediction" << std::endl;
222  std::vector<TH1*> hPrediction;
223  hPrediction.push_back(predictionVec[0]->Predict(calc).ToTH1(pot));
224  for(int quantId = 0; quantId < nQuant; quantId++){
225  hPrediction.push_back(predictionVec[quantId]->Predict(calc).ToTH1(pot));
226  }
227  hPrediction[0]->Add(hPrediction[2]);
228  hPrediction[0]->Add(hPrediction[3]);
229  hPrediction[0]->Add(hPrediction[4]);
230  for(int quantId = 0; quantId < nQuantPlus; quantId++){
231  hPrediction[quantId]->Scale(0.1, "width");
232  hPrediction[quantId]->Add(hCosmics[quantId]); // add in nominal cosmics
233  }
234  std::vector<Spectrum> sPrediction;
235  for(int quantId = 0; quantId < nQuantPlus; quantId++){
236  Spectrum prediction(hPrediction[quantId], pot, livetime);
237  sPrediction.push_back(prediction);
238  }
239 
240 
241 
242  // predictions maximal mixing
243  std::vector<TH1*> hPredictionMaxMix;
244  hPredictionMaxMix.push_back(predictionVec[0]->Predict(calcMaxMix).ToTH1(pot));
245  for(int quantId = 0; quantId < nQuant; quantId++){
246  hPredictionMaxMix.push_back(predictionVec[quantId]->Predict(calcMaxMix).ToTH1(pot));
247  }
248  hPredictionMaxMix[0]->Add(hPredictionMaxMix[2]);
249  hPredictionMaxMix[0]->Add(hPredictionMaxMix[3]);
250  hPredictionMaxMix[0]->Add(hPredictionMaxMix[4]);
251  for(int quantId = 0; quantId < nQuantPlus; quantId++){
252  hPredictionMaxMix[quantId]->Scale(0.1, "width");
253  hPredictionMaxMix[quantId]->Add(hCosmics[quantId]); // add in nominal cosmics
254  }
255  std::vector<Spectrum> sPredictionMaxMix;
256  for(int quantId = 0; quantId < nQuantPlus; quantId++){
257  Spectrum prediction(hPredictionMaxMix[quantId], pot, livetime);
258  sPredictionMaxMix.push_back(prediction);
259  }
260 
261 
262  // predictions sa best fit
263  std::cout << "SA best fit prediction" << std::endl;
264  std::vector<TH1*> hPredictionSA;
265  hPredictionSA.push_back(predictionVec[0]->Predict(calcSA).ToTH1(pot));
266  for(int quantId = 0; quantId < nQuant; quantId++){
267  hPredictionSA.push_back(predictionVec[quantId]->Predict(calcSA).ToTH1(pot));
268  }
269  hPredictionSA[0]->Add(hPredictionSA[2]);
270  hPredictionSA[0]->Add(hPredictionSA[3]);
271  hPredictionSA[0]->Add(hPredictionSA[4]);
272  for(int quantId = 0; quantId < nQuantPlus; quantId++){
273  hPredictionSA[quantId]->Scale(0.1, "width");
274  hPredictionSA[quantId]->Add(hCosmics[quantId]); // add in nominal cosmics
275  }
276  std::vector<Spectrum> sPredictionSA;
277  for(int quantId = 0; quantId < nQuantPlus; quantId++){
278  Spectrum prediction(hPredictionSA[quantId], pot, livetime);
279  sPredictionSA.push_back(prediction);
280  }
281 
282 
283  // predictions no oscillations
284  std::cout << "No oscillations prediction" << std::endl;
285  std::vector<TH1*> hPredictionNoOsc;
286  hPredictionNoOsc.push_back(predictionVec[0]->Predict(calcNoOsc).ToTH1(pot));
287  for(int quantId = 0; quantId < nQuant; quantId++){
288  hPredictionNoOsc.push_back(predictionVec[quantId]->Predict(calcNoOsc).ToTH1(pot));
289  }
290  hPredictionNoOsc[0]->Add(hPredictionNoOsc[2]);
291  hPredictionNoOsc[0]->Add(hPredictionNoOsc[3]);
292  hPredictionNoOsc[0]->Add(hPredictionNoOsc[4]);
293  for(int quantId = 0; quantId < nQuantPlus; quantId++){
294  hPredictionNoOsc[quantId]->Scale(0.1, "width");
295  hPredictionNoOsc[quantId]->Add(hCosmics[quantId]); // add in nominal cosmics
296  }
297 
298 
299 
300  // edited plotWithSystErrorBand
301  // eats the top of user range
302  // and sets yaxis to events/0.1 gev
303  std::vector<float> maxGraph;
304  maxGraph.push_back(13.0);
305  maxGraph.push_back(5.0);
306  maxGraph.push_back(6.0);
307  maxGraph.push_back(4.5);
308  maxGraph.push_back(6.0);
309  std::vector<float> maxChi;
310  maxChi.push_back(5.0);
311  maxChi.push_back(1.5);
312  maxChi.push_back(5.0);
313  maxChi.push_back(2.5);
314  maxChi.push_back(2.5);
315 
316 
317 
318  std::cout << "\nDefining hOne histogram" << std::endl;
319  TH1D* hOne = new TH1D("hOne","hOne",20,0,5);
320  for(unsigned int k=0; k<20; k++){
321  hOne->SetBinContent(k+1,1);
322  }
323  hOne->SetLineStyle(7);
324 
325 
326 
327 
328  std::cout << "\nAbout to plot\n" << std::endl;
329  for(int quantId = 0; quantId < nQuantPlus; quantId++){
330 
331  // pimp in one go
332  hPredictionNoOsc[quantId]->GetXaxis()->SetTitle("Reconstructed Neutrino Energy (GeV)");
333  hPredictionNoOsc[quantId]->GetYaxis()->SetTitle("Events / 0.1 GeV");
334  hPredictionNoOsc[quantId]->GetYaxis()->CenterTitle();
335  hPredictionNoOsc[quantId]->GetXaxis()->CenterTitle();
336  hPredictionNoOsc[quantId]->SetTitle("");
337 
338  hPredictionNoOsc[quantId]->SetLineColor(kGreen+2);
339  hPredictionNoOsc[quantId]->SetLineStyle(kDashed);
340  hPredictionNoOsc[quantId]->SetLineWidth(3);
341  hPredictionSA[quantId]->SetLineColor(kGreen+2);
342  hPredictionSA[quantId]->SetLineStyle(kDashed);
343  hPredictionSA[quantId]->SetLineWidth(3);
344  hPredictionMaxMix[quantId]->SetLineColor(kRed);
345  hPredictionMaxMix[quantId]->SetLineWidth(3);
346  hPrediction[quantId]->SetLineColor(kRed);
347  hPrediction[quantId]->SetLineWidth(3);
348  hCosmics[quantId]->SetLineColor(kGray);
349  hCosmics[quantId]->SetLineWidth(2);
350  hBkg[quantId]->SetLineColor(kBlue);
351  hBkg[quantId]->SetLineStyle(kDashed);
352  hBkg[quantId]->SetLineWidth(2);
353  hAllBkg[quantId]->SetLineColor(kBlue);
354  hAllBkg[quantId]->SetLineStyle(kDashed);
355  hAllBkg[quantId]->SetLineWidth(2);
356 
357 
358  // spectrums for syst band best fit
359  std::vector<Spectrum> ups, dns;
360  for(const ISyst* syst: systs){
361  SystShifts shifts;
362  //if(syst->ShortName() == "cosmicScale") continue;
363  shifts.SetShift(syst, +1);
364  if(quantId==0){
365  TH1 * htempup = predictionVec[0]->PredictSyst(calc, shifts).ToTH1(pot);
366  htempup->Add(predictionVec[1]->PredictSyst(calc, shifts).ToTH1(pot));
367  htempup->Add(predictionVec[2]->PredictSyst(calc, shifts).ToTH1(pot));
368  htempup->Add(predictionVec[3]->PredictSyst(calc, shifts).ToTH1(pot));
369  htempup->Scale(0.1, "width");
370  htempup->Add(hCosmics[0]);
371  Spectrum specup(htempup, pot, livetime);
372  ups.push_back(specup);
373  shifts.SetShift(syst, -1);
374  TH1 * htempdn = predictionVec[0]->PredictSyst(calc, shifts).ToTH1(pot);
375  htempdn->Add(predictionVec[1]->PredictSyst(calc, shifts).ToTH1(pot));
376  htempdn->Add(predictionVec[2]->PredictSyst(calc, shifts).ToTH1(pot));
377  htempdn->Add(predictionVec[3]->PredictSyst(calc, shifts).ToTH1(pot));
378  htempdn->Scale(0.1, "width");
379  htempdn->Add(hCosmics[0]);
380  Spectrum specdn(htempdn, pot, livetime);
381  dns.push_back(specdn);
382  }
383  else{
384  TH1 * htempup = predictionVec[quantId-1]->PredictSyst(calc, shifts).ToTH1(pot);
385  htempup->Scale(0.1, "width");
386  htempup->Add(hCosmics[quantId]);
387  Spectrum specup(htempup, pot, livetime);
388  ups.push_back(specup);
389  shifts.SetShift(syst, -1);
390  TH1 * htempdn = predictionVec[quantId-1]->PredictSyst(calc, shifts).ToTH1(pot);
391  htempdn->Scale(0.1, "width");
392  htempdn->Add(hCosmics[quantId]);
393  Spectrum specdn(htempdn, pot, livetime);
394  dns.push_back(specdn);
395  }
396  }
397 
398 
399 
400  // spectrums for syst band maximal mixing
401  std::vector<Spectrum> upsMaxMix, dnsMaxMix;
402  for(const ISyst* syst: systs){
403  SystShifts shifts;
404  //if(syst->ShortName() == "cosmicScale") continue;
405  shifts.SetShift(syst, +1);
406  if(quantId==0){
407  TH1 * htempup = predictionVec[0]->PredictSyst(calcMaxMix, shifts).ToTH1(pot);
408  htempup->Add(predictionVec[1]->PredictSyst(calcMaxMix, shifts).ToTH1(pot));
409  htempup->Add(predictionVec[2]->PredictSyst(calcMaxMix, shifts).ToTH1(pot));
410  htempup->Add(predictionVec[3]->PredictSyst(calcMaxMix, shifts).ToTH1(pot));
411  htempup->Scale(0.1, "width");
412  htempup->Add(hCosmics[0]);
413  Spectrum specup(htempup, pot, livetime);
414  upsMaxMix.push_back(specup);
415  shifts.SetShift(syst, -1);
416  TH1 * htempdn = predictionVec[0]->PredictSyst(calcMaxMix, shifts).ToTH1(pot);
417  htempdn->Add(predictionVec[1]->PredictSyst(calcMaxMix, shifts).ToTH1(pot));
418  htempdn->Add(predictionVec[2]->PredictSyst(calcMaxMix, shifts).ToTH1(pot));
419  htempdn->Add(predictionVec[3]->PredictSyst(calcMaxMix, shifts).ToTH1(pot));
420  htempdn->Scale(0.1, "width");
421  htempdn->Add(hCosmics[0]);
422  Spectrum specdn(htempdn, pot, livetime);
423  dnsMaxMix.push_back(specdn);
424  }
425  else{
426  TH1 * htempup = predictionVec[quantId-1]->PredictSyst(calcMaxMix, shifts).ToTH1(pot);
427  htempup->Scale(0.1, "width");
428  htempup->Add(hCosmics[quantId]);
429  Spectrum specup(htempup, pot, livetime);
430  upsMaxMix.push_back(specup);
431  shifts.SetShift(syst, -1);
432  TH1 * htempdn = predictionVec[quantId-1]->PredictSyst(calcMaxMix, shifts).ToTH1(pot);
433  htempdn->Scale(0.1, "width");
434  htempdn->Add(hCosmics[quantId]);
435  Spectrum specdn(htempdn, pot, livetime);
436  dnsMaxMix.push_back(specdn);
437  }
438  }
439 
440 
441 
442  // spectrums for syst band sa best fit
443  std::vector<Spectrum> upsSA, dnsSA;
444  for(const ISyst* syst: systs){
445  SystShifts shifts;
446  //if(syst->ShortName() == "cosmicScale") continue;
447  shifts.SetShift(syst, +1);
448  if(quantId==0){
449  TH1 * htempup = predictionVec[0]->PredictSyst(calcSA, shifts).ToTH1(pot);
450  htempup->Add(predictionVec[1]->PredictSyst(calcSA, shifts).ToTH1(pot));
451  htempup->Add(predictionVec[2]->PredictSyst(calcSA, shifts).ToTH1(pot));
452  htempup->Add(predictionVec[3]->PredictSyst(calcSA, shifts).ToTH1(pot));
453  htempup->Scale(0.1, "width");
454  htempup->Add(hCosmics[0]);
455  Spectrum specup(htempup, pot, livetime);
456  upsSA.push_back(specup);
457  shifts.SetShift(syst, -1);
458  TH1 * htempdn = predictionVec[0]->PredictSyst(calcSA, shifts).ToTH1(pot);
459  htempdn->Add(predictionVec[1]->PredictSyst(calcSA, shifts).ToTH1(pot));
460  htempdn->Add(predictionVec[2]->PredictSyst(calcSA, shifts).ToTH1(pot));
461  htempdn->Add(predictionVec[3]->PredictSyst(calcSA, shifts).ToTH1(pot));
462  htempdn->Scale(0.1, "width");
463  htempdn->Add(hCosmics[0]);
464  Spectrum specdn(htempdn, pot, livetime);
465  dnsSA.push_back(specdn);
466  }
467  else{
468  TH1 * htempup = predictionVec[quantId-1]->PredictSyst(calcSA, shifts).ToTH1(pot);
469  htempup->Scale(0.1, "width");
470  htempup->Add(hCosmics[quantId]);
471  Spectrum specup(htempup, pot, livetime);
472  upsSA.push_back(specup);
473  shifts.SetShift(syst, -1);
474  TH1 * htempdn = predictionVec[quantId-1]->PredictSyst(calcSA, shifts).ToTH1(pot);
475  htempdn->Scale(0.1, "width");
476  htempdn->Add(hCosmics[quantId]);
477  Spectrum specdn(htempdn, pot, livetime);
478  dnsSA.push_back(specdn);
479  }
480  }
481 
482 
483 
484  // canvas size
485  int cx = 600; int cy = 500;
486  int cxunosc = 600; int cyunosc = 800;
487  int cxunoscr = 600; int cyunoscr = 400;
488  int cxchi = 600; int cychi = 800;
489 
490 
491  // data + best fit + sa pred + background
492  TString outName = outDir + Form("datapred_all_quant%d",quantId);
493  ofstream file;
494  TCanvas *canvas = new TCanvas("canvas","canvas", cx, cy);
495  canvas->cd();
496  PlotWithSystErrorBand(sPrediction[quantId], ups, dns, pot, -1, -1, maxGraph[quantId], true);
497  hPredictionSA[quantId]->Draw("hist same");
498  hAllBkg[quantId]->Draw("hist same");
499  grData[quantId]->Draw("P0");
500  grData[quantId]->Draw("P0");
501  if(preliminary) Preliminary();
502  //legend
503  TLegend *legend = new TLegend(0.60,0.55,0.85,0.85);
504  legend->AddEntry(grData[quantId], "Data","lep");
505  legend->AddEntry(hPrediction[quantId], "Prediction","l");
506  TLegendEntry *entry=legend->AddEntry("error","1-#sigma syst. range","f");
507  entry->SetFillStyle(1001);
508  entry->SetFillColor(kRed-10);
509  legend->AddEntry(hPredictionSA[quantId], "Previous best fit","l");
510  legend->AddEntry(hAllBkg[quantId], "Background","l");
511  legend->SetTextSize(0.04); //no border for legend
512  legend->SetBorderSize(0); //no border for legend
513  legend->SetFillColor(0); //fill colour is white
514  legend->Draw();
515  drawLabel(quantId);
516  gPad->Modified();
517  canvas->SaveAs(outName + ".eps");
518  canvas->SaveAs(outName + ".png");
519  canvas->SaveAs(outName + ".pdf");
520  file.open (outName + ".txt");
521  if(quantId==0) file << "Third analysis spectrum showing data, prediction at best fit with systematics band about stats-only best fit, prediction at previous best fit, and backgrounds.";
522  if(quantId>0) file<< "Third analysis spectrum showing data, prediction at best fit with systematics band about stats-only best fit, prediction at previous best fit, and backg\
523 rounds in quantile " << quantId << ".";
524  file.close();
525 
526 
527  // data + best fit + beam background + cosmic background
528  outName = outDir + Form("datapred_bestfit_quant%d",quantId);
529  ofstream file3;
530  TCanvas *canvas3 = new TCanvas("canvas3","canvas3", cx, cy);
531  canvas3->cd();
532  PlotWithSystErrorBand(sPrediction[quantId], ups, dns, pot, -1, -1, maxGraph[quantId], true);
533  hBkg[quantId]->Draw("hist same");
534  hCosmics[quantId]->Draw("hist same");
535  grData[quantId]->Draw("P0");
536  grData[quantId]->Draw("P0");
537  if(preliminary) Preliminary();
538  //legend
539  TLegend *legend3 = new TLegend(0.60,0.60,0.85,0.85);
540  legend3->AddEntry(grData[quantId], "Data","lep");
541  legend3->AddEntry(hPrediction[quantId], "Prediction","l");
542  TLegendEntry *entry3=legend3->AddEntry("error","1-#sigma syst. range","f");
543  entry3->SetFillStyle(1001);
544  entry3->SetFillColor(kRed-10);
545  legend3->AddEntry(hBkg[quantId], "Beam bkg.","l");
546  legend3->AddEntry(hCosmics[quantId], "Cosmic bkg.","l");
547  legend3->SetTextSize(0.04); //no border for legend
548  legend3->SetBorderSize(0); //no border for legend
549  legend3->SetFillColor(0); //fill colour is white
550  legend3->Draw();
551  drawLabel(quantId);
552  canvas3->SaveAs(outName + ".eps");
553  canvas3->SaveAs(outName + ".png");
554  canvas3->SaveAs(outName + ".pdf");
555  file3.open (outName + ".txt");
556  if(quantId==0) file3 << "Third analysis spectrum showing data, best fit prediction, systematics band about stats-only best fit, and beam and cosmic backgrounds.";
557  if(quantId>0) file3 << "Third analysis spectrum showing data, best fit prediction, systematics band about stats-only best fit, and beam and cosmic backgrounds in quantile " << quantId << ".";
558  file3.close();
559 
560 
561  // data + max mix predition + sa prediction
562  outName = outDir + Form("datapred_samaxmix_quant%d",quantId);
563  ofstream file2;
564  TCanvas *canvas2 = new TCanvas("canvas2","canvas2", cx, cy);
565  canvas2->cd();
566  PlotWithSystErrorBand(sPredictionMaxMix[quantId], upsMaxMix, dnsMaxMix, pot, -1, -1, maxGraph[quantId], true);
567  hPredictionSA[quantId]->Draw("hist same");
568  grData[quantId]->Draw("P0");
569  grData[quantId]->Draw("P0");
570  if(preliminary) Preliminary();
571  //legend
572  TLegend *legend2 = new TLegend(0.60,0.65,0.85,0.85);
573  legend2->AddEntry(grData[quantId], "Data","lep");
574  legend2->AddEntry(hPredictionMaxMix[quantId], "Max. mix. pred.","l");
575  TLegendEntry *entry2=legend2->AddEntry("error","1-#sigma syst. range","f");
576  entry2->SetFillStyle(1001);
577  entry2->SetFillColor(kRed-10);
578  legend2->AddEntry(hPredictionSA[quantId], "Previous best fit","l");
579  legend2->SetTextSize(0.04); //no border for legend
580  legend2->SetBorderSize(0); //no border for legend
581  legend2->SetFillColor(0); //fill colour is white
582  legend2->Draw();
583  drawLabel(quantId);
584  canvas2->SaveAs(outName + ".eps");
585  canvas2->SaveAs(outName + ".png");
586  canvas2->SaveAs(outName + ".pdf");
587  file2.open (outName + ".txt");
588  if(quantId==0) file2 << "Third analysis spectrum showing data, maximal mixing prediction with systematics band and prediction at previous best fit.";
589  if(quantId>0) file2 << "Third analysis spectrum showing data, maximal mixing prediction with systematics band and prediction at previous best fit in quantile " << quantId << ".";
590  file2.close();
591 
592 
593  if(plotunosc){
594 
595  // prediction ratios
596  TH1 *nooscdenom = hNoOsc[quantId];
597  nooscdenom->Add(hBkg[quantId],-1);
598  TH1 *predratio = hPred[quantId];
599  predratio->Add(hBkg[quantId],-1);
600  predratio->Divide(nooscdenom);
601  predratio->GetYaxis()->SetTitle("Reconstructed Neutrino Energy (GeV)");
602  predratio->GetYaxis()->SetTitle("Ratio to no oscillation");
603  predratio->GetYaxis()->CenterTitle();
604  predratio->GetXaxis()->CenterTitle();
605  predratio->SetLineColor(kRed);
606  predratio->SetLineWidth(3);
607  predratio->SetMaximum(1.5);
608  predratio->SetMinimum(-0.001);
609 
610 
611  // data ratio
612  TFeldmanCousins fc(0.6827);
613  TGraphAsymmErrors *dataratio = GraphWithPoissonErrors(hData[quantId]);
614  dataratio->SetMarkerStyle(kFullCircle);
615  dataratio->SetLineWidth(2);
616  for(int bin = 0; bin < hData[quantId]->GetNbinsX(); bin++) {
617  double oldx, oldy;
618  dataratio->GetPoint(bin,oldx,oldy);
619  oldy = oldy - hCosmics[quantId]->GetBinContent(bin+1) - hBkg[quantId]->GetBinContent(bin+1);
620  dataratio->SetPoint(bin,oldx,oldy/nooscdenom->GetBinContent(bin+1));
621  dataratio->SetPointEYhigh(bin,(fc.CalculateUpperLimit(oldy,0)-oldy)/nooscdenom->GetBinContent(bin+1));
622  dataratio->SetPointEYlow(bin,(oldy-fc.CalculateLowerLimit(oldy,0))/nooscdenom->GetBinContent(bin+1));
623  if(bin==1||bin==0) {
624  dataratio->SetPoint(bin,oldx,-50);
625  dataratio->SetPointEYhigh(bin,0);
626  dataratio->SetPointEYlow(bin,0);
627  }
628  }
629 
630 
631 
632  TString outNameUnosc = outDir + Form("spect_unosc_quant%d",quantId);
633  ofstream fileUnosc;
634  TCanvas *canunosc = new TCanvas("canunosc","canunosc", cxunosc, cyunosc);
635  canunosc->cd();
636  hPredictionNoOsc[quantId]->Draw("hist same");
637  hPrediction[quantId]->Draw("hist same");
638  hAllBkg[quantId]->Draw("hist same");
639  grData[quantId]->Draw("p0 same");
640  if(preliminary) Preliminary();
641  TLegend *legunosc = new TLegend(0.52,0.38,0.85,0.53);
642  legunosc->SetFillColor(10);
643  legunosc->SetFillStyle(0);
644  legunosc->SetBorderSize(0);
645  legunosc->SetTextSize(0.04);
646  legunosc->AddEntry(grData[quantId],"Data","lep");
647  legunosc->AddEntry(hAllBkg[quantId],"Background","l");
648  legunosc->AddEntry(hPrediction[quantId],"Prediction","l");
649  legunosc->AddEntry(hPredictionNoOsc[quantId],"Unoscilated pred.","l");
650  legunosc->Draw();
651  drawLabelUnosc(quantId);
652  canunosc->SaveAs(outNameUnosc + ".eps");
653  canunosc->SaveAs(outNameUnosc + ".pdf");
654  canunosc->SaveAs(outNameUnosc + ".png");
655  fileUnosc.open (outNameUnosc + ".txt");
656  if(quantId==0)fileUnosc << "Spectrum plot showing the unoscillated prediction, data, and the best fit prediction (no systematics).";
657  if(quantId>0)fileUnosc << "Spectrum plot showing the unoscillated prediction, data, and the best fit prediction (no systematics) in quantile " << quantId << ".";
658  fileUnosc.close();
659 
660 
661  outNameUnosc = outDir + Form("ratio_noosc_quant%d",quantId);
662  ofstream fileUnoscRatio;
663  TCanvas *canratio = new TCanvas("canratio","canratio", cxunoscr, cyunoscr);
664  canratio->cd();
665  predratio->Draw("hist same");
666  dataratio->Draw("pe same");
667  hOne->Draw("hist same");
668  if(preliminary) Preliminary();
669  TLegend *legratio = new TLegend(0.30,0.45,0.50,0.60);
670  legratio->SetFillColor(10);
671  legratio->SetFillStyle(0);
672  legratio->SetBorderSize(0);
673  legratio->SetTextSize(0.05);
674  legratio->AddEntry(dataratio,"Data","lep");
675  legratio->AddEntry(predratio,"Prediction","l");
676  legratio->Draw();
677  drawLabelRatio(quantId);
678  canratio->SaveAs(outNameUnosc + ".eps");
679  canratio->SaveAs(outNameUnosc + ".pdf");
680  canratio->SaveAs(outNameUnosc + ".png");
681  fileUnoscRatio.open (outNameUnosc + ".txt");
682  if(quantId==0) fileUnoscRatio << "Spectrum ratio of data over unoscillated prediction and best fit prediction over unoscillated prediction. All spectra have the predicted backgrounds subtracted out before division";
683  if(quantId>0) fileUnoscRatio << "Spectrum ratio of data over unoscillated prediction and best fit prediction over unoscillated prediction for quantile " << quantId << ". All spectra have the predicted backgrounds subtracted out before division";
684  fileUnoscRatio.close();
685 
686  } // if plotunoscilated
687 
688 
689 
690  if(plotchi){
691 
692  double binEdges[20] = { 0, 0.75, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.25, 2.5, 2.75, 3, 3.5, 4.0, 5.0};
693  TH1F *hChisq = new TH1F("hChisq","hChisq",19, binEdges);
694  float sumchi = 0;
695  for(int i = 1; i < 20; i++) {
696  float e = hPrediction[quantId]->GetBinContent(i);
697  float o = hData[quantId]->GetBinContent(i);
698  float tempchi = LogLikelihood(e, o);
699  hChisq->SetBinContent(i,tempchi);
700  sumchi+=tempchi;
701  }
702 
703  hChisq->SetMaximum(maxChi[quantId]);
704  hChisq->SetLineColor(kBlack);
705  hChisq->SetMarkerStyle(kFullCircle);
706  hChisq->SetMarkerSize(1);
707  hChisq->SetLineWidth(2);
708  hChisq->SetTitle("");
709  hChisq->GetXaxis()->SetTitle("Reconstructed Neutrino Energy (GeV)");
710  hChisq->GetYaxis()->SetTitle("#chi^{2} contribution");
711  hChisq->GetYaxis()->CenterTitle();
712  hChisq->GetXaxis()->CenterTitle();
713 
714 
715  TString outNameChisq = outDir + Form("chisq_quant%d",quantId);
716  ofstream fileChisq;
717  TCanvas *canchisq = new TCanvas("canchisq","canchisq", cxchi , cychi);
718  canchisq->cd();
719  TPad *pad1 = new TPad("spectrum","spectrum",0.0,0.4,1.0,1.0);
720  TPad *pad2 = new TPad("chisq","chisq",0.0,0.0,1.0,0.4);
721  pad1->SetBottomMargin(0.0);
722  pad1->SetLeftMargin(0.12);
723  pad1->SetRightMargin(0.05);
724  pad2->SetTopMargin(0.0);
725  pad2->SetBottomMargin(0.20);
726  pad2->SetLeftMargin(0.12);
727  pad2->SetRightMargin(0.05);
728  pad1->Draw();
729  pad2->Draw();
730  pad1->cd();
731  gPad->SetTickx(1);
732  gPad->SetTicky(1);
733  PlotWithSystErrorBand(sPrediction[quantId], ups, dns, pot, -1, -1, maxGraph[quantId], true);
734  hAllBkg[quantId]->Draw("hist same");
735  grData[quantId]->Draw("P0");
736  grData[quantId]->Draw("P0");
737  if(preliminary) Preliminary();
738  TLegend *legchisq = new TLegend(0.63,0.63,0.88,0.85);
739  legchisq->SetFillColor(10);
740  legchisq->SetFillStyle(0);
741  legchisq->SetBorderSize(0);
742  legchisq->SetTextSize(0.04);
743  legchisq->AddEntry(grData[quantId],"Data","lep");
744  legchisq->AddEntry(hPrediction[quantId],"Prediction","l");
745  TLegendEntry *enchisq=legchisq->AddEntry("error","1-#sigma syst. range","f");
746  enchisq->SetFillStyle(1001);
747  enchisq->SetFillColor(kRed-10);
748  legchisq->AddEntry(hAllBkg[quantId],"Background","l");
749  legchisq->Draw();
750  drawLabel(quantId);
751  pad2->cd();
752  gPad->SetTickx(1);
753  gPad->SetTicky(1);
754  hChisq->Draw("hist");
755  canchisq->SaveAs(outNameChisq + ".eps");
756  canchisq->SaveAs(outNameChisq + ".pdf");
757  canchisq->SaveAs(outNameChisq + ".png");
758  fileChisq.open (outNameChisq + ".txt");
759  if(quantId==0) fileChisq << "The top panel shows the third analysis spectrum including data, best fit prediction, systematics band about stats-only best fit and backgrounds. The bottom panel shows the chi-sq contribution to the final fit of each energy bin.";
760  if(quantId>0) fileChisq << "The top panel shows the third analysis spectrum including data, best fit prediction, systematics band about stats-only best fit and backgrounds. \
761 The bottom panel shows the chi-sq contribution to the final fit of each energy bin in quantile " << quantId << ".";
762  fileChisq.close();
763 
764  } // if plotchi
765 
766 
767 
768 
769  // clear shifts at the end as
770  // might be used for the unosc
771  ups.clear();
772  dns.clear();
773  upsSA.clear();
774  dnsSA.clear();
775  upsMaxMix.clear();
776  dnsMaxMix.clear();
777 
778 
779 
780  } // for quantiles
781 
782 
783 
784 
785 } // End of function
void drawLabel(int quantId=0)
Definition: tools.h:132
Pass neutrinos through unchanged.
Definition: IOscCalc.h:40
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
TGraphAsymmErrors * GraphWithPoissonErrors(const TH1 *h, bool noErrorsXaxis, bool drawEmptyBins)
Calculate statistical errors appropriate for small Poisson numbers.
Definition: Plots.cxx:910
void drawLabelRatio(int quantId=0)
Definition: tools.h:223
osc::OscCalcDumb calc
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
Definition: EigenUtils.cxx:36
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:33
TGraphAsymmErrors * PlotWithSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float headroom, bool newaxis, EBinType bintype)
Plot prediction with +/-1sigma error band.
Definition: Plots.cxx:304
TString inDir
void plotDataPred(bool newpred=true, bool plotunosc=true, bool plotchi=true, bool plotshifts=true)
Definition: plotDataPred.C:8
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:44
#define pot
const int nQuantPlus
Definition: varsandcuts.h:5
Int_t preliminary
Definition: SimpleIterate.C:63
float bin[41]
Definition: plottest35.C:14
c1 Close()
Oscillation probability calculators.
Definition: Calcs.h:5
void Preliminary()
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double livetime
Definition: saveFDMCHists.C:21
const int nQuant
Definition: varsandcuts.h:4
virtual void SetTh23(const T &th23)=0
void drawLabelUnosc(int quantId=0)
Definition: tools.h:180
Neutral-current interactions.
Definition: IPrediction.h:40
TFile * file
Definition: cellShifts.C:17
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
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
Float_t e
Definition: plot.C:35
def entry(str)
Definition: HTMLTools.py:26
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:78
TPad * pad1
Definition: analysis.C:13
T asin(T number)
Definition: d0nt_math.hpp:60