demo_numuOnly.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 void demo_numuOnly(const TString options = "demo1")
3 {
4  std::cout << "Sorry, you must run in compiled mode" << std::endl;
5 }
6 #else
7 
8 
10 #include "CAFAna/Analysis/Fit.h"
11 #include "CAFAna/Analysis/Surface.h"
13 #include "CAFAna/Core/Spectrum.h"
14 #include "CAFAna/FC/FCSurface.h"
20 #include "CAFAna/Vars/FitVars.h"
22 #include "TFile.h"
23 #include "TCanvas.h"
24 #include "TH2.h"
25 
26 using namespace ana;
27 
28 
30  Spectrum&,
31  TString, bool);
32 
33 void demo_numuOnly(const TString options = "demo1")
34 {
35 
36  TString predFileName = "pred_tutorial_nue_numu.root";
37  TFile * predFile = new TFile(predFileName,"read");
38  auto predictionNumu = ana::LoadFrom<IPrediction>(predFile->GetDirectory("prediction_numu"));
40  const double potFD = 6E20;
41 
42  if(options.Contains("demo1")) {
43 
44  calc->SetTh23(asin(sqrt(0.4)));
45  auto fakeDataNumu = predictionNumu->Predict(calc).FakeData(potFD);
46 
48 
49  PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeDataNumu, "demo1", false);
50  }
51 
52  if(options.Contains("demo2")){
53 
54  for(int potfactor:{1,10,100}){
55  double pot = potFD * potfactor;
56  //Create different types of fake data
58  auto fakeData_max = predictionNumu->Predict(calc).FakeData(pot);
59  auto mockData_max = predictionNumu->Predict(calc).MockData(pot);
60  calc->SetTh23(asin(sqrt(0.4)));
61  auto fakeData_nonMax = predictionNumu->Predict(calc).FakeData(pot);
62 
63  //reduce copy paste
64  std::vector <TString> labels = {"fakeMax","mockMax","fakeNonMax"};
65  std::vector <Spectrum> data = {fakeData_max,mockData_max,fakeData_nonMax};
66 
67  for(int i=0;i<3;++i){
69  PlotNumuPredData(predictionNumu.get(), calc, pot, data[i], "demo2_"+labels[i]+"_potX"+std::to_string(potfactor), true);
70  auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), data[i]);
71 
72  std::cout << "Fit: prediction maximal mixing, data " << labels[i] << std::endl;
74  //We are not fitting any parameters
75  //Whatever is in the calculator, is used for the prediction
76  Fitter fitter (exptNumu,{},{});
77  fitter.Fit(calc);
78  }
79  }
80  }
81 
82  ////////////////////////////////////////////////////////////
83  if(options.Contains("demo3")){
84 
86  calc->SetTh23(asin(sqrt(0.4)));
87  auto fakeData = predictionNumu->Predict(calc).FakeData(potFD);
88  auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeData);
89 
91  PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeData, "demo3_before", true);
92 
93  Fitter fitter (exptNumu,{&kFitSinSqTheta23},{});
94  fitter.Fit(calc);
95  //Fit result is printed on screen and saved in calc :)
96  PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeData, "demo3_after_th23", true);
97 
98  //Crazy alternative: fit everything but th23
100  Fitter fitter2 (exptNumu,{&kFitSinSq2Theta13, &kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta12, &kFitDmSq21},{});
101  fitter2.Fit(calc);
102  PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeData, "demo3_after_other", true);
103  }
104 
105  ////////////////////////////////////////////////////////////
106  if(options.Contains("demo4")){
107 
108  ResetOscCalcToDefault(calc);
109  calc->SetTh23(asin(sqrt(0.4)));
110  auto fakeData = predictionNumu->Predict(calc).FakeData(potFD);
111  auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeData);
112 
113  auto reactor = WorldReactorConstraint2015();
114  auto solar = new SolarConstraints();
115 
116  auto exptNumuConstraint = new MultiExperiment({exptNumu,reactor,solar});
117 
118  ResetOscCalcToDefault(calc);
119  PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeData, "demo4_before", true);
120 
121  Fitter fitter1 (exptNumuConstraint,{&kFitSinSq2Theta13, &kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta12, &kFitDmSq21},{});
122  fitter1.Fit(calc);
123  PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeData, "demo4_constraints_noth23", true);
124 
125  ResetOscCalcToDefault(calc);
126  Fitter fitter2 (exptNumuConstraint,{&kFitSinSqTheta23, &kFitSinSq2Theta13, &kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta12, &kFitDmSq21},{});
127  fitter2.Fit(calc);
128  PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeData, "demo4_constraints_all", true);
129 
130  ResetOscCalcToDefault(calc);
131  Fitter fitter3 (exptNumu,{ &kFitSinSqTheta23, &kFitSinSq2Theta13, &kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta12, &kFitDmSq21},{});
132  fitter3.Fit(calc);
133  PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeData, "demo4_noconstraints_all", true);
134  }
135 
136  if(options.Contains("demo5")){
137 
138  ResetOscCalcToDefault(calc);
139  calc->SetTh23(asin(sqrt(0.4)));
140  auto fakeData = predictionNumu->Predict(calc).FakeData(potFD);
141  auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeData);
142 
143  ResetOscCalcToDefault(calc);
144 
145  auto prof_th23 = Profile (exptNumu, calc, &kFitSinSqTheta23, 30, 0, 1,
146  -1, {});
147 
148  prof_th23->SetMaximum(450);
149  prof_th23->Draw("hist");
150  gPad->Print("plot_demo5_profile_th23.pdf");
151 
152  ResetOscCalcToDefault(calc);
153  auto prof_th23_dmsq = Profile (exptNumu, calc, &kFitSinSqTheta23, 30, 0, 1,
154  -1, {&kFitDmSq32});
155 
156  prof_th23_dmsq->SetMaximum(450);
157  prof_th23_dmsq->Draw("hist");
158  gPad->Print("plot_demo5_profile_th23_dmsq.pdf");
159 
160  ResetOscCalcToDefault(calc);
161 
162  auto prof_dCP = Profile (exptNumu, calc, &kFitDeltaInPiUnits, 30, 0, 2,
163  -1, {});
164 
165  prof_dCP->SetMaximum(20);
166  prof_dCP->Draw("hist");
167  gPad->Print("plot_demo5_profile_dCP.pdf");
168 
169  ResetOscCalcToDefault(calc);
170  auto prof_dCP_th23_dmsq = Profile (exptNumu, calc, &kFitDeltaInPiUnits, 30, 0, 2,
171  -1, {&kFitDmSq32,&kFitSinSqTheta23});
172 
173  prof_dCP_th23_dmsq->SetMaximum(20);
174  prof_dCP_th23_dmsq->Draw("hist");
175  gPad->Print("plot_demo5_profile_dCP_th23_dmsq.pdf");
176  }
177 
178  if(options.Contains("demo6")){
179 
180  ResetOscCalcToDefault(calc);
181  calc->SetTh23(asin(sqrt(0.4)));
182  auto fakeData = predictionNumu->Predict(calc).FakeData(potFD);
183  auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeData);
184  auto reactor = WorldReactorConstraint2015();
185 
186  auto exptNumuConstraint = new MultiExperiment({exptNumu,reactor});
187  ResetOscCalcToDefault(calc);
188 
189  Surface surf (exptNumu, calc,
190  &kFitSinSqTheta23, 30, 0, 1,
191  &kFitDmSq32Scaled, 30, 2,3.5,
193  new TCanvas;
194  surf.Draw();
195  gPad->Print("plot_demo6_surface.pdf");
196  auto margs = surf.GetProfiledHists();
197  margs[0]->Draw("colz");
198  gPad->Print("plot_demo6_surface_prof_dCP.pdf");
199 
200  margs[1]->Draw("colz");
201  gPad->Print("plot_demo6_surface_prof_th13.pdf");
202 
203 
204  ResetOscCalcToDefault(calc);
205 
206  Surface surf2 (exptNumuConstraint, calc,
207  &kFitSinSqTheta23, 30, 0, 1,
208  &kFitDmSq32Scaled, 30, 2,3.5,
210  {{&kFitDeltaInPiUnits,{0.,0.5,1.,1.5}}});
211  new TCanvas;
212  surf2.Draw();
213  gPad->Print("plot_demo6_surface2.pdf");
214  auto margs2 = surf2.GetProfiledHists();
215  margs2[0]->Draw("colz");
216  gPad->Print("plot_demo6_surface2_prof_dCP.pdf");
217 
218  margs2[1]->Draw("colz");
219  gPad->Print("plot_demo6_surface2_prof_th13.pdf");
220  }
221  if(options.Contains("demo7")){
222 
223  ResetOscCalcToDefault(calc);
224  calc->SetTh23(asin(sqrt(0.4)));
225  auto fakeData = predictionNumu->Predict(calc).FakeData(potFD);
226  auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeData);
227  ResetOscCalcToDefault(calc);
228  Surface surf (exptNumu, calc,
229  &kFitSinSqTheta23, 30, 0, 1,
230  &kFitDmSq32Scaled, 30, 2,3.5);
231  new TCanvas;
232  TH2 * surf_1Sigma = Gaussian68Percent2D(surf);
233  TH2 * surf_2Sigma = Gaussian2Sigma2D(surf);
234  TH2 * surf_3Sigma = Gaussian3Sigma2D(surf);
235  surf.DrawContour(surf_1Sigma, kSolid, kRed);
236  surf.DrawContour(surf_2Sigma, kSolid, kViolet);
237  surf.DrawContour(surf_3Sigma, kSolid, kBlue);
238  gPad->Print("plot_demo7_contours.pdf");
239 
240 
241  }
242 
243 
244 
245 }
246 
247 
248 
250  double potFD,
251  Spectrum & fakeDataNumu,
252  TString suffix, bool sameCanvas=false){
253  new TCanvas;
254 
255  auto hpred = predictionNumu->Predict(calc).ToTH1(potFD,kRed);
256  auto hfake = fakeDataNumu.ToTH1(potFD);
257 
258  hpred->SetMaximum(1.2*max(hpred->GetMaximum(),hfake->GetMaximum()));
259  hfake->SetMaximum(1.2*max(hpred->GetMaximum(),hfake->GetMaximum()));
260 
261  hfake->SetMarkerStyle(kFullCircle);
262  hfake->SetTitle("Fake Data");
263 
264  if(sameCanvas){
265  hpred->Draw("hist");
266  hfake->Draw("same");
267  gPad->Print("plot_"+suffix+"_overlay.pdf");
268  }
269  else{
270  hpred->SetTitle("Prediction");
271  hpred->Draw("hist");
272  gPad->Print("plot_"+suffix+"_pred.pdf");
273  hfake->Draw("p");
274  gPad->Print("plot_"+suffix+"_data.pdf");
275  }
276 }
277 
278 
279 
280 
281 #endif
T max(const caf::Proxy< T > &a, T b)
const FitSinSq2Theta12 kFitSinSq2Theta12
Definition: FitVars.cxx:23
TGraph * Profile(const IChiSqExperiment *expt, osc::IOscCalculatorAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double input_minchi, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &profSystsMap, MinuitFitter::FitOpts opts)
scan in one variable, profiling over all others
Definition: Fit.cxx:48
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
Constraints on the parameters and from solar experiments.
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
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:16
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
const ReactorExperiment * WorldReactorConstraint2015()
Weighted average of all experiments as of first nue paper writing.
void ResetOscCalcToDefault(osc::IOscCalculatorAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
TH2 * Gaussian2Sigma2D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 2D in gaussian approximation.
void demo_numuOnly(const TString options="demo1")
Definition: demo_numuOnly.C:33
const FitDmSq21 kFitDmSq21
Definition: FitVars.cxx:24
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
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
General interface to any calculator that lets you set the parameters.
void PlotNumuPredData(IPrediction *, osc::IOscCalculatorAdjustable *, double, Spectrum &, TString, bool)
#define pot
TH2 * Gaussian3Sigma2D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 2D in gaussian approximation.
const FitDmSq32Scaled kFitDmSq32Scaled
Definition: FitVars.cxx:17
OStream cout
Definition: OStream.cxx:6
Combine multiple component experiments.
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir)
Definition: IPrediction.cxx:38
const FitSinSqTheta23 kFitSinSqTheta23
Definition: FitVars.cxx:15
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const FitDeltaInPiUnits kFitDeltaInPiUnits
Definition: FitVars.cxx:14
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
surf
Definition: demo5.py:35