genie_contours.C
Go to the documentation of this file.
1 #include "CAFAna/Core/Registry.h"
3 #include "CAFAna/Fit/Fit.h"
6 #include "CAFAna/Core/Spectrum.h"
7 #include "CAFAna/FC/FCSurface.h"
13 #include "CAFAna/Systs/Systs.h"
19 #include "CAFAna/Systs/BeamSysts.h"
21 #include "CAFAna/Vars/FitVars.h"
22 #include "OscLib/IOscCalc.h"
23 #include "TFile.h"
24 #include "TCanvas.h"
25 #include "TH2.h"
26 #include "TLegend.h"
27 #include "TLatex.h"
28 #include "TRandom.h"
29 #include "TRandom3.h"
30 
31 using namespace ana;
32 /*void make_contours(MultiExperiment* expt, osc::IOscCalcAdjustable* calc, std::vector<const ISyst*> pc_syst,
33  FrequentistSurface surf_numu, TH1* prof_numu, int npc)
34 {
35  std::vector<const ISyst*> numuSysts = pc_syst;
36 
37  // ResetOscCalcToDefault(calc);
38  FrequentistSurface surf_numu_syst(expt, calc, &kFitSinSqTheta23, 50, 0, 1,
39  &kFitDmSq32Scaled, 50, 2, 3.5,
40  {&kFitDeltaInPiUnits, &kFitSinSq2Theta13},numuSysts,
41  {{&kFitDeltaInPiUnits,{0.,0.5,1.,1.5}}});
42 
43  TCanvas* cnumu = new TCanvas("numu", "numu", 200, 10, 1600, 600);
44  cnumu->Divide(2,1);
45  cnumu->cd(1);
46  surf_numu.DrawContour(Gaussian90Percent2D(surf_numu), kSolid, kBlue);
47  surf_numu.DrawBestFit(kBlue);
48  surf_numu_syst.DrawContour(Gaussian90Percent2D(surf_numu_syst), kDashed, kRed);
49  surf_numu_syst.DrawBestFit(kRed);
50 
51  TH1* prof_numu_syst = Profile(expt, calc, &kFitSinSqTheta23, 50, 0.3, 0.7, -1,
52  {&kFitDmSq32Scaled, &kFitSinSq2Theta13}, numuSysts);
53  cnumu->cd(2);
54  prof_numu->SetLineColor(kBlue);
55  prof_numu_syst->SetLineColor(kRed);
56  prof_numu_syst->SetLineStyle(kDashed);
57  TLegend* leg_numu = new TLegend(0.5, 0.4, 0.7, 0.8);
58  std::string cv_leg_numu = "CV : "+std::to_string(sqrt(prof_numu->Interpolate(0.5)));
59  std::string pc_leg_numu = "PC"+std::to_string(npc)+" : "+std::to_string(sqrt(prof_numu_syst->Interpolate(0.5)));
60  leg_numu->AddEntry(prof_numu, cv_leg_numu.c_str(), "l");
61  leg_numu->AddEntry(prof_numu_syst, pc_leg_numu.c_str(), "l");
62  prof_numu->Draw("hist");
63  prof_numu->SetTitle("#theta_{23} Sensitivity");
64  prof_numu_syst->Draw("hist same");
65  leg_numu->Draw();
66  std::string numufile = "contours/numu_pc"+std::to_string(npc)+".root";
67  std::string numufile2 = "contours/numu_pc"+std::to_string(npc)+".png";
68  cnumu->Print(numufile.c_str());
69  cnumu->Print(numufile2.c_str());
70 
71 
72  // ResetOscCalcToDefault(calc);
73  // FrequentistSurface surf_nue_syst(expt, calc,
74  // &kFitDeltaInPiUnits,100,0,2,
75  // &kFitSinSqTheta23, 50, 0, 1,
76  // {&kFitDmSq32Scaled, &kFitSinSq2Theta13}, nueSysts);
77  //
78  // TCanvas* cnue = new TCanvas("nue", "nue", 200, 10, 1600, 600);
79  // cnue->Divide(2,1);
80  // cnue->cd(1);
81  // surf_nue.DrawContour(Gaussian90Percent2D(surf_nue), kSolid, kBlue);
82  // surf_nue.DrawBestFit(kBlue);
83  // surf_nue_syst.DrawContour(Gaussian90Percent2D(surf_nue_syst), kDashed, kRed);
84  // surf_nue_syst.DrawBestFit(kRed);
85  //
86  //
87  // TH1* prof_nue_syst = Profile(expt, calc, &kFitDeltaInPiUnits, 100, 0, 2, -1,
88  // {&kFitDmSq32Scaled, &kFitSinSq2Theta13, &kFitSinSqTheta23}, numuSysts);
89  // cnue->cd(2);
90  // prof_nue->SetLineColor(kBlue);
91  // prof_nue_syst->SetLineColor(kRed);
92  // prof_nue_syst->SetLineStyle(kDashed);
93  // TLegend* leg_nue = new TLegend(0.3, 0.4, 0.5, 0.8);
94  // std::string cv_leg_nue = "CV : "+std::to_string(sqrt(prof_nue->Interpolate(0)));
95  // std::string pc_leg_nue = "PC"+std::to_string(npc)+" : "+std::to_string(sqrt(prof_nue_syst->Interpolate(0)));
96  // leg_nue->AddEntry(prof_nue, cv_leg_nue.c_str(), "l");
97  // leg_nue->AddEntry(prof_nue_syst, pc_leg_nue.c_str(), "l");
98  // prof_nue->Draw("hist");
99  // prof_nue->SetTitle("CP Sensitivity");
100  // prof_nue_syst->Draw("hist same");
101  // leg_nue->Draw();
102  // std::string nuefile = "contours/nue_pc"+std::to_string(npc)+".root";
103  // std::string nuefile2 = "contours/nue_pc"+std::to_string(npc)+".png";
104  // cnue->Print(nuefile.c_str());
105  // cnue->Print(nuefile2.c_str());
106 
107 }
108 
109 void make_nom_expt(TFile* predFile, std::string nue_folder, std::string numu_folder, const Spectrum& nue_data, const Spectrum& numu_data)
110 {
111  auto nuePred = ana::LoadFrom<IPrediction>(predFile, nue_folder.c_str());
112  auto numuPred = ana::LoadFrom<IPrediction>(predFile, numu_folder.c_str());
113 
114  osc::IOscCalcAdjustable* calc = DefaultOscCalc();
115  const double potFD = 9E20;
116 
117  ResetOscCalcToDefault(calc);
118  calc->SetTh23(asin(sqrt(0.56)));
119  calc->SetDmsq32(2.48e-3);
120  calc->SetdCP(.0*M_PI);
121  auto nueExpt = new SingleSampleExperiment(nuePred.get(), nue_data);
122  auto numuExpt = new SingleSampleExperiment(numuPred.get(), numu_data);
123  auto reactor = WorldReactorConstraint2019();
124  auto expt = new MultiExperiment({nueExpt, numuExpt, reactor});
125 
126  // ResetOscCalcToDefault(calc);
127  FrequentistSurface surf_numu(expt, calc, &kFitSinSqTheta23, 50, 0, 1,
128  &kFitDmSq32Scaled, 50, 2, 3.5,
129  {&kFitDeltaInPiUnits, &kFitSinSq2Theta13},{},
130  {{&kFitDeltaInPiUnits, {0.,0.5,1.,1.5}}});
131 
132  surf_numu.DrawContour(Gaussian90Percent2D(surf_numu), kSolid, kGray);
133  surf_numu.DrawBestFit(kGray);
134 
135  gPad->Update();
136 
137 }
138 
139 TH1* make_expt_profile(TFile* predFile, std::string nue_folder, std::string numu_folder, const Spectrum& nue_data, const Spectrum& numu_data)
140 {
141  auto nuePred = ana::LoadFrom<IPrediction>(predFile, nue_folder.c_str());
142  auto numuPred = ana::LoadFrom<IPrediction>(predFile, numu_folder.c_str());
143 
144  osc::IOscCalcAdjustable* calc = DefaultOscCalc();
145  const double potFD = 9E20;
146 
147  ResetOscCalcToDefault(calc);
148  calc->SetTh23(asin(sqrt(0.56)));
149  calc->SetDmsq32(2.48e-3);
150  calc->SetdCP(.0*M_PI);
151  auto nueExpt = new SingleSampleExperiment(nuePred.get(), nue_data);
152  auto numuExpt = new SingleSampleExperiment(numuPred.get(), numu_data);
153  auto reactor = WorldReactorConstraint2019();
154  auto expt = new MultiExperiment({nueExpt, numuExpt, reactor});
155 
156  // ResetOscCalcToDefault(calc);
157  TH1* prof_numu = Profile(expt, calc, &kFitSinSqTheta23, 50, 0.3, 0.7, -1,
158  {&kFitDmSq32Scaled, &kFitSinSq2Theta13}, {});
159  prof_numu->SetLineColor(kGray+1);
160  return prof_numu;
161 }*/
162 
163 double get_univ_chisq(TFile* predFileNue, TFile* predFileNumu, std::string nue_fhc_folder, std::string numu_fhc_folder, std::string nue_rhc_folder, std::string numu_rhc_folder, const Spectrum& nue_fhc_data, const Spectrum& numu_fhc_data, const Spectrum& nue_rhc_data, const Spectrum& numu_rhc_data, const SystShifts shift, osc::IOscCalcAdjustable* calc)
164 {
165  // auto nuePred_fhc = ana::LoadFrom<PredictionInterp>(predFile, nue_fhc_folder.c_str());
166  // auto numuPred_fhc = ana::LoadFrom<PredictionInterp>(predFile, numu_fhc_folder.c_str());
167  // auto nuePred_rhc = ana::LoadFrom<PredictionInterp>(predFile, nue_rhc_folder.c_str());
168  // auto numuPred_rhc = ana::LoadFrom<PredictionInterp>(predFile, numu_rhc_folder.c_str());
169  auto nuePred_fhc = PredictionInterp::LoadFrom(predFileNue , nue_fhc_folder.c_str()).release();
170  auto numuPred_fhc = PredictionInterp::LoadFrom(predFileNumu, numu_fhc_folder.c_str()).release();
171  auto nuePred_rhc = PredictionInterp::LoadFrom(predFileNue , nue_rhc_folder.c_str()).release();
172  auto numuPred_rhc = PredictionInterp::LoadFrom(predFileNumu, numu_rhc_folder.c_str()).release();
173 
174  // osc::IOscCalcAdjustable* calc = DefaultOscCalc();
175  const double potFD_fhc = 14E20;
176  const double potFD_rhc = 12.33E20;
177 
178  auto nueExpt_fhc = new SingleSampleExperiment(nuePred_fhc, nue_fhc_data);
179  auto numuExpt_fhc = new SingleSampleExperiment(numuPred_fhc, numu_fhc_data);
180  auto nueExpt_rhc = new SingleSampleExperiment(nuePred_rhc, nue_rhc_data);
181  auto numuExpt_rhc = new SingleSampleExperiment(numuPred_rhc, numu_rhc_data);
182  auto reactor = WorldReactorConstraint2019();
183  auto expt = new MultiExperiment({nueExpt_fhc, numuExpt_fhc, nueExpt_rhc, numuExpt_rhc, reactor});
184 
185  // ResetOscCalcToDefault(calc);
186  return sqrt(expt->ChiSq(calc, shift));
187 }
188 
190 {
191 
192  // osc::IOscCalcAdjustable* calc = DefaultOscCalc();
193  // ResetOscCalcToDefault(calc);
194  return sqrt(expt->ChiSq(calc, shift));
195 }
196 
197 void genie_contours(int npc = 5)
198 {
199 
200  //TString filename1 = "genie_preds_000_499.root";
201  //TString filename2 = "genie_preds_500_999.root";
202  TString filenue = "genie_syst_pred_nue_PCs.root";
203  TString filenue100 = "genie_syst_pred_nue_univs_0_100.root";
204  TString filenue200 = "genie_syst_pred_nue_univs_100_200.root";
205  TString filenue300 = "genie_syst_pred_nue_univs_200_300.root";
206  TString filenue400 = "genie_syst_pred_nue_univs_300_400.root";
207  TString filenue500 = "genie_syst_pred_nue_univs_400_500.root";
208  TString filenue600 = "genie_syst_pred_nue_univs_500_600.root";
209  TString filenue700 = "genie_syst_pred_nue_univs_600_700.root";
210  TString filenue800 = "genie_syst_pred_nue_univs_700_800.root";
211  TString filenue900 = "genie_syst_pred_nue_univs_800_900.root";
212  TString filenue1000 = "genie_syst_pred_nue_univs_900_1000.root";
213  TString filenumu = "genie_syst_pred_numu.root";
214  std::vector<const ISyst*> pc_no_sup_shifts;
215  std::vector<const ISyst*> pc_fhc_sup_shifts;
216  std::vector<const ISyst*> pc_fhc_ws_sup_shifts;
217  // for(int i = 0; i < 10; i++)
218  // pc_fhc_sup_shifts.push_back(GetPPFXPrincipalsFHCSup(i));
219  // for(int i = 0; i < 10; i++)
220  // pc_fhc_ws_sup_shifts.push_back(GetPPFXPrincipalsFHCWSSup(i)); for(int i = 0; i < 10; i++)
221  for(int i = 0; i < 20; i++)
222  pc_no_sup_shifts.push_back(GetGeniePrincipals2020Small(i));
223  //TFile* predFile = new TFile(filename1, "read");
224  //TFile* predFile2 = new TFile(filename2, "read");
225  TFile* predFileNue = new TFile(filenue, "read");
226  TFile* predFileNue100 = new TFile(filenue100, "read");
227  TFile* predFileNue200 = new TFile(filenue200, "read");
228  TFile* predFileNue300 = new TFile(filenue300, "read");
229  TFile* predFileNue400 = new TFile(filenue400, "read");
230  TFile* predFileNue500 = new TFile(filenue500, "read");
231  TFile* predFileNue600 = new TFile(filenue600, "read");
232  TFile* predFileNue700 = new TFile(filenue700, "read");
233  TFile* predFileNue800 = new TFile(filenue800, "read");
234  TFile* predFileNue900 = new TFile(filenue900, "read");
235  TFile* predFileNue1000 = new TFile(filenue1000, "read");
236  TFile* predFileNumu = new TFile(filenumu, "read");
237 
238  std::cout << FindCAFAnaDir()+"/data/xs/" << std::endl;
239  std::cout << "point 1" << std::endl;
240 
241  auto nuePred_fhc = PredictionInterp::LoadFrom(predFileNue , "nue_prop_fhc").release();
242  auto numuPred_fhc = PredictionInterp::LoadFrom(predFileNumu, "numu_fhc").release();
243  auto nuePred_rhc = PredictionInterp::LoadFrom(predFileNue , "nue_prop_rhc").release();
244  auto numuPred_rhc = PredictionInterp::LoadFrom(predFileNumu, "numu_rhc").release();
245 
246  std::cout << "point 2" << std::endl;
247 
249  const double potFD_fhc = 14E20;
250  const double potFD_rhc = 12.33E20;
251 
252  ResetOscCalcToDefault(calc);
253  calc->SetTh23(asin(sqrt(0.56)));
254  calc->SetDmsq32(2.48e-3);
255  calc->SetdCP(.0*M_PI);
256  auto nueFakeData_fhc = nuePred_fhc ->Predict(calc).FakeData(potFD_fhc);
257  auto numuFakeData_fhc = numuPred_fhc->Predict(calc).FakeData(potFD_fhc);
258  auto nueExpt_fhc = new SingleSampleExperiment(nuePred_fhc, nueFakeData_fhc);
259  auto numuExpt_fhc = new SingleSampleExperiment(numuPred_fhc, numuFakeData_fhc);
260  auto nueFakeData_rhc = nuePred_rhc ->Predict(calc).FakeData(potFD_rhc);
261  auto numuFakeData_rhc = numuPred_rhc->Predict(calc).FakeData(potFD_rhc);
262  auto nueExpt_rhc = new SingleSampleExperiment(nuePred_rhc, nueFakeData_rhc);
263  auto numuExpt_rhc = new SingleSampleExperiment(numuPred_rhc, numuFakeData_rhc);
264  auto reactor = WorldReactorConstraint2019();
265  auto expt = new MultiExperiment({nueExpt_fhc, numuExpt_fhc, nueExpt_rhc, numuExpt_rhc, reactor});
266 
267  std::cout << "point 3" << std::endl;
268 
269  /*
270  std::string title = "NH, 2.48e-3";
271  std::string label = "nh";
272  */
273  //calc->SetdCP(0);
274 
275  std::string title = "IH, -2.48e-3";
276  std::string label = "ih";
277  calc->SetDmsq32(-2.48e-3);
278 
279  /*
280  std::string title = "Max-Mixing Rejection";
281  std::string label = "mmix";
282  calc->SetTh23(asin(sqrt(0.5)));
283  */
284  /*
285  std::string title = "IH Rejection";
286  std::string label = "ih";
287  calc->SetDmsq32(-2.48e-3);
288  */
289  // ResetOscCalcToDefault(calc);
290  // for(int npc_tot = 6; npc_tot > 0; npc_tot--){
291 
292  int npc_tot = npc;
293  double rlow = 4.0;
294  double rhigh = 5.0;
295  TH1D* sens_syst_no_sup = new TH1D("sens_syst_no_sup", "sens_syst_no_sup", 200, rlow, rhigh);
296  //gRandom->SetSeed(0);
297 
298  std::cout << "point 4" << std::endl;
299  int tempseed = 0;
300  for(int i = 0; i < 10000; i++){
301  SystShifts shift_no_sup;
302  for(int npc = 0; npc < npc_tot; npc++){
303  //double shift = 2.0*gRandom->Gaus();
304  TRandom3 rand (tempseed++);
305  //shift_no_sup.SetShift(pc_no_sup_shifts[npc], shift);
306  shift_no_sup.SetShift(pc_no_sup_shifts[npc], (double) rand.Gaus());
307  }
308  //std::cout << "-----" << std::endl;
309  double chisq_no_sup = get_chisq(expt, shift_no_sup, calc);
310  sens_syst_no_sup->Fill(chisq_no_sup);
311  std::cout << i << " " << chisq_no_sup << std::endl;
312  }
313  std::cout << "point 5" << std::endl;
314  TH1D* sensitivity = new TH1D("sensitivity", "sensitivity", 20, rlow, rhigh);
315  std::string fnameUnivs = "genie_univs_hist_"+label+".root";
316  TFile fUnivs(fnameUnivs.c_str(),"read");
317  if(fUnivs.IsZombie())
318  {
319  for(int i = 0; i < 1000; i++){
320  //std::cout << "univ " << i << std::endl;
321  double chisq;
322  if (i<100)
323  chisq = get_univ_chisq(predFileNue100, predFileNumu, "nue_prop_univ_fhc_"+std::to_string(i), "numu_univ_fhc_"+std::to_string(i), "nue_prop_univ_rhc_"+std::to_string(i), "numu_univ_rhc_"+std::to_string(i), nueFakeData_fhc, numuFakeData_fhc, nueFakeData_rhc, numuFakeData_rhc, kNoShift, calc);
324  else if (i<200)
325  chisq = get_univ_chisq(predFileNue200, predFileNumu, "nue_prop_univ_fhc_"+std::to_string(i), "numu_univ_fhc_"+std::to_string(i), "nue_prop_univ_rhc_"+std::to_string(i), "numu_univ_rhc_"+std::to_string(i), nueFakeData_fhc, numuFakeData_fhc, nueFakeData_rhc, numuFakeData_rhc, kNoShift, calc);
326  else if (i<300)
327  chisq = get_univ_chisq(predFileNue300, predFileNumu, "nue_prop_univ_fhc_"+std::to_string(i), "numu_univ_fhc_"+std::to_string(i), "nue_prop_univ_rhc_"+std::to_string(i), "numu_univ_rhc_"+std::to_string(i), nueFakeData_fhc, numuFakeData_fhc, nueFakeData_rhc, numuFakeData_rhc, kNoShift, calc);
328  else if (i<400)
329  chisq = get_univ_chisq(predFileNue400, predFileNumu, "nue_prop_univ_fhc_"+std::to_string(i), "numu_univ_fhc_"+std::to_string(i), "nue_prop_univ_rhc_"+std::to_string(i), "numu_univ_rhc_"+std::to_string(i), nueFakeData_fhc, numuFakeData_fhc, nueFakeData_rhc, numuFakeData_rhc, kNoShift, calc);
330  else if (i<500)
331  chisq = get_univ_chisq(predFileNue500, predFileNumu, "nue_prop_univ_fhc_"+std::to_string(i), "numu_univ_fhc_"+std::to_string(i), "nue_prop_univ_rhc_"+std::to_string(i), "numu_univ_rhc_"+std::to_string(i), nueFakeData_fhc, numuFakeData_fhc, nueFakeData_rhc, numuFakeData_rhc, kNoShift, calc);
332  else if (i<600)
333  chisq = get_univ_chisq(predFileNue600, predFileNumu, "nue_prop_univ_fhc_"+std::to_string(i), "numu_univ_fhc_"+std::to_string(i), "nue_prop_univ_rhc_"+std::to_string(i), "numu_univ_rhc_"+std::to_string(i), nueFakeData_fhc, numuFakeData_fhc, nueFakeData_rhc, numuFakeData_rhc, kNoShift, calc);
334  else if (i<700)
335  chisq = get_univ_chisq(predFileNue700, predFileNumu, "nue_prop_univ_fhc_"+std::to_string(i), "numu_univ_fhc_"+std::to_string(i), "nue_prop_univ_rhc_"+std::to_string(i), "numu_univ_rhc_"+std::to_string(i), nueFakeData_fhc, numuFakeData_fhc, nueFakeData_rhc, numuFakeData_rhc, kNoShift, calc);
336  else if (i<800)
337  chisq = get_univ_chisq(predFileNue800, predFileNumu, "nue_prop_univ_fhc_"+std::to_string(i), "numu_univ_fhc_"+std::to_string(i), "nue_prop_univ_rhc_"+std::to_string(i), "numu_univ_rhc_"+std::to_string(i), nueFakeData_fhc, numuFakeData_fhc, nueFakeData_rhc, numuFakeData_rhc, kNoShift, calc);
338  else if (i<900)
339  chisq = get_univ_chisq(predFileNue900, predFileNumu, "nue_prop_univ_fhc_"+std::to_string(i), "numu_univ_fhc_"+std::to_string(i), "nue_prop_univ_rhc_"+std::to_string(i), "numu_univ_rhc_"+std::to_string(i), nueFakeData_fhc, numuFakeData_fhc, nueFakeData_rhc, numuFakeData_rhc, kNoShift, calc);
340  else if (i<1000)
341  chisq = get_univ_chisq(predFileNue1000, predFileNumu, "nue_prop_univ_fhc_"+std::to_string(i), "numu_univ_fhc_"+std::to_string(i), "nue_prop_univ_rhc_"+std::to_string(i), "numu_univ_rhc_"+std::to_string(i), nueFakeData_fhc, numuFakeData_fhc, nueFakeData_rhc, numuFakeData_rhc, kNoShift, calc);
342  //if( i==499 ) continue;
343  //if( i==999 ) continue;
344  //double chisq = i < 500 ?
345  // get_univ_chisq(predFile, "nue_prop_univ_fhc_"+std::to_string(i), "numu_univ_fhc_"+std::to_string(i), "nue_prop_univ_rhc_"+std::to_string(i), "numu_univ_rhc_"+std::to_string(i), nueFakeData_fhc, numuFakeData_fhc, nueFakeData_rhc, numuFakeData_rhc, kNoShift, calc) :
346  // get_univ_chisq(predFile2, "nue_prop_univ_fhc_"+std::to_string(i), "numu_univ_fhc_"+std::to_string(i), "nue_prop_univ_rhc_"+std::to_string(i), "numu_univ_rhc_"+std::to_string(i), nueFakeData_fhc, numuFakeData_fhc, nueFakeData_rhc, numuFakeData_rhc, kNoShift, calc) ;
347  //double chisq = get_univ_chisq(predFileNue, predFileNumu, "nue_prop_univ_fhc_"+std::to_string(i), "numu_univ_fhc_"+std::to_string(i), "nue_prop_univ_rhc_"+std::to_string(i), "numu_univ_rhc_"+std::to_string(i), nueFakeData_fhc, numuFakeData_fhc, nueFakeData_rhc, numuFakeData_rhc, kNoShift, calc);
348  sensitivity->Fill(chisq);
349  std::cout << i << " " << chisq << std::endl;
350  }
351  sensitivity->Scale(1./sensitivity->Integral());
352  sensitivity->SetMaximum(2.0*sensitivity->GetMaximum());
353  } else {
354  sensitivity = (TH1D*) fUnivs.Get("sensitivity")->Clone();
355  }
356 
357  std::cout << "point 6" << std::endl;
358 
359  TCanvas *c_no_sup = new TCanvas("c_no_sup", "c_no_sup", 200, 10, 800, 600);
360  c_no_sup->cd();
361  //sensitivity->Scale(1./sensitivity->Integral());
362  //sensitivity->SetMaximum(2.0*sensitivity->GetMaximum());
363  sensitivity->Draw("hist same");
364  sens_syst_no_sup->Scale(sens_syst_no_sup->GetNbinsX()/(sensitivity->GetNbinsX()*sens_syst_no_sup->Integral()));
365  sens_syst_no_sup->Draw("hist same");
366  sensitivity->SetTitle(Form("%s w/ %d PC", title.c_str(), npc_tot));
367  sensitivity->GetXaxis()->SetTitle("#sigma");
368  std::cout << "PC Mean : " << sens_syst_no_sup->GetMean() << std::endl;
369  std::cout << "PC RMS : " << sens_syst_no_sup->GetRMS() << std::endl;
370  std::cout << "Genie Univ Mean : " << sensitivity->GetMean() << std::endl;
371  std::cout << "Genie Univ RMS : " << sensitivity->GetRMS() << std::endl;
372  sens_syst_no_sup->SetLineColor(kRed);
373  c_no_sup->Update();
374  c_no_sup->Print(Form("plots/%s_joint_no_sup_%dpc.pdf", label.c_str(), npc_tot));
375 
376  /*
377  TFile* f = new TFile("genie_univs_hist_cp.root", "recreate");
378  f->cd();
379  sensitivity->Write();
380  f->Close();
381  */
382  // TCanvas *c_fhc_sup = new TCanvas("c_fhc_sup", "c_fhc_sup", 200, 10, 800, 600);
383  // c_fhc_sup->cd();
384  // sensitivity->Scale(1./sensitivity->Integral());
385  // sensitivity->Draw("hist same");
386  // sens_syst_fhc_sup->Scale(sens_syst_fhc_sup->GetNbinsX()/(sensitivity->GetNbinsX()*sens_syst_fhc_sup->Integral()));
387  // sens_syst_fhc_sup->Draw("hist same");
388  // sensitivity->SetTitle(Form("CP violation w/ %d PC", npc_tot));
389  // sensitivity->GetXaxis()->SetTitle("#sigma");
390  // std::cout << "PC Mean : " << sens_syst_fhc_sup->GetMean() << std::endl;
391  // std::cout << "PC RMS : " << sens_syst_fhc_sup->GetRMS() << std::endl;
392  // std::cout << "PPFX Univ Mean : " << sensitivity->GetMean() << std::endl;
393  // std::cout << "PPFX Univ RMS : " << sensitivity->GetRMS() << std::endl;
394  // sens_syst_fhc_sup->SetLineColor(kRed);
395  // c_fhc_sup->Update();
396  // c_fhc_sup->Print(Form("deltacp_fhc_only_fhc_sup_%dpc.pdf", npc_tot));
397  //
398  // TCanvas *c_fhc_ws_sup = new TCanvas("c_fhc_ws_sup", "c_fhc_ws_sup", 200, 10, 800, 600);
399  // c_fhc_ws_sup->cd();
400  // sensitivity->Scale(1./sensitivity->Integral());
401  // sensitivity->Draw("hist same");
402  // sens_syst_fhc_ws_sup->Scale(sens_syst_fhc_ws_sup->GetNbinsX()/(sensitivity->GetNbinsX()*sens_syst_fhc_ws_sup->Integral()));
403  // sens_syst_fhc_ws_sup->Draw("hist same");
404  // sensitivity->SetTitle(Form("CP violation w/ %d PC", npc_tot));
405  // sensitivity->GetXaxis()->SetTitle("#sigma");
406  // std::cout << "PC Mean : " << sens_syst_fhc_ws_sup->GetMean() << std::endl;
407  // std::cout << "PC RMS : " << sens_syst_fhc_ws_sup->GetRMS() << std::endl;
408  // std::cout << "PPFX Univ Mean : " << sensitivity->GetMean() << std::endl;
409  // std::cout << "PPFX Univ RMS : " << sensitivity->GetRMS() << std::endl;
410  // sens_syst_fhc_ws_sup->SetLineColor(kRed);
411  // c_fhc_ws_sup->Update();
412  // c_fhc_ws_sup->Print(Form("deltacp_fhc_only_fhc_ws_sup_%dpc.pdf", npc_tot));
413  // }
414 
415 
416 }
virtual double ChiSq(osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const
Definition: IExperiment.h:18
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
T sqrt(T number)
Definition: d0nt_math.hpp:156
double get_univ_chisq(TFile *predFileNue, TFile *predFileNumu, std::string nue_fhc_folder, std::string numu_fhc_folder, std::string nue_rhc_folder, std::string numu_rhc_folder, const Spectrum &nue_fhc_data, const Spectrum &numu_fhc_data, const Spectrum &nue_rhc_data, const Spectrum &numu_rhc_data, const SystShifts shift, osc::IOscCalcAdjustable *calc)
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
std::string FindCAFAnaDir()
Definition: Utilities.cxx:204
virtual void SetDmsq32(const T &dmsq32)=0
#define M_PI
Definition: SbMath.h:34
void genie_contours(int npc=5)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
expt
Definition: demo5.py:34
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
double get_chisq(IExperiment *expt, const SystShifts shift, osc::IOscCalcAdjustable *calc)
Combine multiple component experiments.
static std::unique_ptr< PredictionInterp > LoadFrom(TDirectory *dir, const std::string &name)
virtual void SetTh23(const T &th23)=0
Base class defining interface for experiments.
Definition: IExperiment.h:14
const ReactorExperiment * WorldReactorConstraint2019()
Reactor constraint from PDG 2019.
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Float_t e
Definition: plot.C:35
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:81
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
GeniePCASyst * GetGeniePrincipals2020Small(int PCIdx)
enum BeamMode string