demo_nueNumuSysts.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 void demo_nueNumuSysts(const TString options = "demoS1")
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/Systs/Systs.h"
21 #include "CAFAna/Systs/NumuSysts.h"
22 #include "CAFAna/Vars/FitVars.h"
24 #include "TFile.h"
25 #include "TCanvas.h"
26 #include "TH2.h"
27 
28 using namespace ana;
29 
30 
31 void demo_nueNumuSysts(const TString options = "demoS1")
32 {
33 
34  TString predFileName = "pred_tutorial_nue_numu.root";
35  TFile * predFile = new TFile(predFileName,"read");
36  auto predictionNumu = ana::LoadFrom<IPrediction>(predFile->GetDirectory("prediction_systs_numu"));
37  auto predictionNue = ana::LoadFrom<IPrediction>(predFile->GetDirectory("prediction_systs_nue"));
38 
40  const double potFD = 6E20;
41 
42  //copied from make_prediction
43  MECq0ShapeSyst2017 kMECq0ShapeSyst2017;
44  std::vector< const ISyst * > numuSysts = {&kMECq0ShapeSyst2017, &kNumuNCScaleSyst, &kBeamAll};
45  std::vector< const ISyst * > nueSysts = {&kMECq0ShapeSyst2017, &kBeamAllTransport, &kNA49};
46 
47  if(options.Contains("demoS1")) {
48  for(TString part:{"numu","nue"}){
49  IPrediction * pred = predictionNumu.get();
50  std::vector< const ISyst * > systs = numuSysts;
51  if(part=="nue"){
52  pred = predictionNue.get();
53  systs = nueSysts;
54  }
55  for(auto syst:systs){
56  auto hnom = pred->Predict(calc).ToTH1(potFD);
57  hnom->SetMaximum(1.5*hnom->GetMaximum());
58  hnom->Draw("hist");
59  for(double sigma:{-2.,-1.,1.,2.}){
60  SystShifts shift(syst,sigma);
61  int color = (sigma>0?kRed:kBlue);
62  auto hss = pred->PredictSyst(calc,shift).ToTH1(potFD,color,fabs(sigma));
63  hss->Draw("hist same");
64  }//end sigma
65  gPad->Print((TString)"plot_demoS1_"+part+"Syst_"+syst->ShortName()+".pdf");
66 
67  }//end syst
68  }//end particle
69 
70 
71 
72  }//end demos1
73 
74  if(options.Contains("demoS2")){
76  calc->SetTh23(asin(sqrt(0.4)));
77  auto fakeData = predictionNumu->Predict(calc).FakeData(potFD);
78  auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeData);
80  Surface surf (exptNumu, calc,
81  &kFitSinSqTheta23, 30, 0, 1,
82  &kFitDmSq32Scaled, 30, 2,3.5,{},numuSysts);
83 
84  TH2 * surf_1Sigma = Gaussian68Percent2D(surf);
85  TH2 * surf_2Sigma = Gaussian2Sigma2D(surf);
86  TH2 * surf_3Sigma = Gaussian3Sigma2D(surf);
87  surf.DrawContour(surf_1Sigma, kSolid, kRed);
88  surf.DrawContour(surf_2Sigma, kSolid, kViolet);
89  surf.DrawContour(surf_3Sigma, kSolid, kBlue);
90  gPad->Print("plot_demoS2_contours_syst.pdf");
91  }
92 
93 
94  if(options.Contains("demoS3")){
96  calc->SetTh23(asin(sqrt(0.4)));
97  auto fakeDataNumu = predictionNumu->Predict(calc).FakeData(potFD);
98  auto fakeDataNue = predictionNue->Predict(calc).FakeData(potFD);
99  auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeDataNumu);
100  auto exptNue = new SingleSampleExperiment(predictionNue.get(), fakeDataNue);
101  auto reactor = WorldReactorConstraint2015(); auto exptAll = new MultiExperiment({exptNumu,exptNue,reactor});
102 
103  ResetOscCalcToDefault(calc);
104  Surface surf (exptAll, calc,
105  &kFitSinSqTheta23, 30, 0, 1,
106  &kFitDmSq32Scaled, 30, 2,3.5,
108  {{&kFitDeltaInPiUnits,{0.,0.5,1.,1.5}}});
109 
110  TH2 * surf_1Sigma = Gaussian68Percent2D(surf);
111  TH2 * surf_2Sigma = Gaussian2Sigma2D(surf);
112  TH2 * surf_3Sigma = Gaussian3Sigma2D(surf);
113  surf.DrawContour(surf_1Sigma, kSolid, kRed);
114  surf.DrawContour(surf_2Sigma, kSolid, kViolet);
115  surf.DrawContour(surf_3Sigma, kSolid, kBlue);
116  gPad->Print("plot_demoS3_contours_joint.pdf");
117 
118 
119  ResetOscCalcToDefault(calc);
120  Surface surf2 (exptAll, calc,
121  &kFitDeltaInPiUnits,30,0,2,
122  &kFitSinSqTheta23, 30, 0, 1,
123  {&kFitDmSq32Scaled, &kFitSinSq2Theta13});
124 
125  TH2 * surf2_1Sigma = Gaussian68Percent2D(surf2);
126  TH2 * surf2_2Sigma = Gaussian2Sigma2D(surf2);
127  TH2 * surf2_3Sigma = Gaussian3Sigma2D(surf2);
128  surf2.DrawContour(surf2_1Sigma, kSolid, kRed);
129  surf2.DrawContour(surf2_2Sigma, kSolid, kViolet);
130  surf2.DrawContour(surf2_3Sigma, kSolid, kBlue);
131  gPad->Print("plot_demoS3_contours_joint_dcp.pdf");
132  }
133 
134  if(options.Contains("demoS4")){
135  ResetOscCalcToDefault(calc);
136  calc->SetTh23(asin(sqrt(0.4)));
137  auto fakeDataNumu = predictionNumu->Predict(calc).FakeData(potFD);
138  auto fakeDataNue = predictionNue->Predict(calc).FakeData(potFD);
139  auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeDataNumu);
140  auto exptNue = new SingleSampleExperiment(predictionNue.get(), fakeDataNue);
141  auto reactor = WorldReactorConstraint2015(); auto exptAll = new MultiExperiment({exptNumu,exptNue,reactor});
142 
143  std::vector<const ISyst*> allSysts= {&kMECq0ShapeSyst2017, &kNumuNCScaleSyst, &kNA49};
144  //numu expt correlations
145  exptAll->SetSystCorrelations(0,
146  {
147  {&kNA49, &kBeamAll}
148  });
149  //nue expt correlations
150  exptAll->SetSystCorrelations(1,
151  {
152  {&kNumuNCScaleSyst,NULL}
153  });
154 
155  ResetOscCalcToDefault(calc);
156  Surface surf (exptAll, calc,
157  &kFitSinSqTheta23, 20, 0, 1,
158  &kFitDmSq32Scaled, 20, 2,3.5,
160  {{&kFitDeltaInPiUnits,{0.,0.5,1.,1.5}}});
161 
162  TH2 * surf_1Sigma = Gaussian68Percent2D(surf);
163  TH2 * surf_2Sigma = Gaussian2Sigma2D(surf);
164  TH2 * surf_3Sigma = Gaussian3Sigma2D(surf);
165  surf.DrawContour(surf_1Sigma, kSolid, kRed);
166  surf.DrawContour(surf_2Sigma, kSolid, kViolet);
167  surf.DrawContour(surf_3Sigma, kSolid, kBlue);
168  gPad->Print("plot_demoS4_contours_joint_systs.pdf");
169  }
170 
171 
172 
173 /*
174  if(options.Contains("demo6")){
175 
176  ResetOscCalcToDefault(calc);
177  calc->SetTh23(asin(sqrt(0.4)));
178  auto fakeData = predictionNumu->Predict(calc).FakeData(potFD);
179  auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeData);
180  auto reactor = WorldReactorConstraint2015();
181 
182  auto exptNumuConstraint = new MultiExperiment({exptNumu,reactor});
183  ResetOscCalcToDefault(calc);
184 
185  Surface surf (exptNumu, calc,
186  &kFitSinSqTheta23, 30, 0, 1,
187  &kFitDmSq32Scaled, 30, 2,3.5,
188  {&kFitDeltaInPiUnits, &kFitSinSq2Theta13});
189  new TCanvas;
190  surf.Draw();
191  gPad->Print("plot_demo6_surface.pdf");
192  auto margs = surf.GetProfiledHists();
193  margs[0]->Draw("colz");
194  gPad->Print("plot_demo6_surface_prof_dCP.pdf");
195 
196  margs[1]->Draw("colz");
197  gPad->Print("plot_demo6_surface_prof_th13.pdf");
198 
199 
200  ResetOscCalcToDefault(calc);
201 
202  Surface surf2 (exptNumuConstraint, calc,
203  &kFitSinSqTheta23, 30, 0, 1,
204  &kFitDmSq32Scaled, 30, 2,3.5,
205  {&kFitDeltaInPiUnits, &kFitSinSq2Theta13},{},
206  {{&kFitDeltaInPiUnits,{0.,0.5,1.,1.5}}});
207  new TCanvas;
208  surf2.Draw();
209  gPad->Print("plot_demo6_surface2.pdf");
210  auto margs2 = surf2.GetProfiledHists();
211  margs2[0]->Draw("colz");
212  gPad->Print("plot_demo6_surface2_prof_dCP.pdf");
213 
214  margs2[1]->Draw("colz");
215  gPad->Print("plot_demo6_surface2_prof_th13.pdf");
216  }
217  if(options.Contains("demo7")){
218 
219  ResetOscCalcToDefault(calc);
220  calc->SetTh23(asin(sqrt(0.4)));
221  auto fakeData = predictionNumu->Predict(calc).FakeData(potFD);
222  auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeData);
223  ResetOscCalcToDefault(calc);
224  Surface surf (exptNumu, calc,
225  &kFitSinSqTheta23, 30, 0, 1,
226  &kFitDmSq32Scaled, 30, 2,3.5);
227  new TCanvas;
228  TH2 * surf_1Sigma = Gaussian68Percent2D(surf);
229  TH2 * surf_2Sigma = Gaussian2Sigma2D(surf);
230  TH2 * surf_3Sigma = Gaussian3Sigma2D(surf);
231  surf.DrawContour(surf_1Sigma, kSolid, kRed);
232  surf.DrawContour(surf_2Sigma, kSolid, kViolet);
233  surf.DrawContour(surf_3Sigma, kSolid, kBlue);
234  gPad->Print("plot_demo7_contours.pdf");
235 
236 
237  }
238 
239 */
240 
241 }
242 
243 
244 
245 
246 
247 #endif
const NumuNCScaleSyst kNumuNCScaleSyst
Definition: NumuSysts.cxx:36
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
std::vector< SystGroupDef > systs
Definition: syst_header.h:384
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
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
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
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.
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
virtual Spectrum PredictSyst(osc::IOscCalculator *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:98
void demo_nueNumuSysts(const TString options="demoS1")
TString part[npart]
Definition: Style.C:32
General interface to any calculator that lets you set the parameters.
const BeamSyst kBeamAllTransport((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"totErr","beamTransportComb","Combined Beam Transport Systematics")
All Beam Transport systematics combined in quadratures.
Definition: BeamSysts.h:122
TH2 * Gaussian3Sigma2D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 2D in gaussian approximation.
static const BeamSyst & kBeamAll
Definition: BeamSysts.h:124
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 NOvARwgtSyst kMECq0ShapeSyst2017("MECq0Shape","MEC q_{0} shape", novarwgt::kMECq0ShapeSyst2017)
Definition: MECSysts.h:27
const FitSinSqTheta23 kFitSinSqTheta23
Definition: FitVars.cxx:15
string syst
Definition: plotSysts.py:176
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
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