getSensitivity.C
Go to the documentation of this file.
1 
3 #include "CAFAna/Fit/Fit.h"
6 #include "CAFAna/Core/Loaders.h"
9 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Cuts/Cuts.h"
12 #include "CAFAna/Cuts/SpillCuts.h"
13 #include "CAFAna/Cuts/TimingCuts.h"
29 #include "CAFAna/Vars/FitVars.h"
33 #include "CAFAna/Vars/Vars.h"
34 
35 #include "OscLib/OscCalcPMNSOpt.h"
36 
37 #include "TCanvas.h"
38 #include "TH1.h"
39 #include "TLegend.h"
40 #include "TStyle.h"
41 
42 using namespace ana;
43 
44 
46 {
47 
48  // get predictions:
49  TString predDir = "/pnfs/nova/persistent/users/lvinton1/numu/pd/numu/3A/predsFakeND/";
50  TString q1Pred = "pred3A_Quant1_realNDData_full_nonMax_cutAllPeriods_tmp2017Syst_noq0shp_wDirAbsHadE_S17-09-24.root";
51  TString q2Pred = "pred3A_Quant2_realNDData_full_nonMax_cutAllPeriods_tmp2017Syst_noq0shp_wDirAbsHadE_S17-09-24.root";
52  TString q3Pred = "pred3A_Quant3_realNDData_full_nonMax_cutAllPeriods_tmp2017Syst_noq0shp_wDirAbsHadE_S17-09-24.root";
53  TString q4Pred = "pred3A_Quant4_realNDData_full_nonMax_cutAllPeriods_tmp2017Syst_noq0shp_wDirAbsHadE_S17-09-24.root";
54 
55  TString q1Pred_maxMix = "pred3A_Quant1_realNDData_full_nonMax_cutAllPeriods_tmp2017Syst_noq0shp_wDirAbsHadE_S17-09-24.root";
56  TString q2Pred_maxMix = "pred3A_Quant2_realNDData_full_nonMax_cutAllPeriods_tmp2017Syst_noq0shp_wDirAbsHadE_S17-09-24.root";
57  TString q3Pred_maxMix = "pred3A_Quant3_realNDData_full_nonMax_cutAllPeriods_tmp2017Syst_noq0shp_wDirAbsHadE_S17-09-24.root";
58  TString q4Pred_maxMix = "pred3A_Quant4_realNDData_full_nonMax_cutAllPeriods_tmp2017Syst_noq0shp_wDirAbsHadE_S17-09-24.root";
59 
60  DontAddDirectory guard;
61  TFile* fPred_Q1 = new TFile(predDir + q1Pred);
62  TFile* fPred_Q2 = new TFile(predDir + q2Pred);
63  TFile* fPred_Q3 = new TFile(predDir + q3Pred);
64  TFile* fPred_Q4 = new TFile(predDir + q4Pred);
65  TFile* fPred_maxMix_Q1 = new TFile(predDir + q1Pred_maxMix);
66  TFile* fPred_maxMix_Q2 = new TFile(predDir + q2Pred_maxMix);
67  TFile* fPred_maxMix_Q3 = new TFile(predDir + q3Pred_maxMix);
68  TFile* fPred_maxMix_Q4 = new TFile(predDir + q4Pred_maxMix);
69 
70  PredictionInterp* prediction_Q1 = LoadFrom<PredictionInterp>(fPred_Q1, "prediction").release();
71  PredictionInterp* prediction_Q2 = LoadFrom<PredictionInterp>(fPred_Q2, "prediction").release();
72  PredictionInterp* prediction_Q3 = LoadFrom<PredictionInterp>(fPred_Q3, "prediction").release();
73  PredictionInterp* prediction_Q4 = LoadFrom<PredictionInterp>(fPred_Q4, "prediction").release();
74 
75  PredictionInterp* prediction_maxMix_Q1 = LoadFrom<PredictionInterp>(fPred_maxMix_Q1, "prediction").release();
76  PredictionInterp* prediction_maxMix_Q2 = LoadFrom<PredictionInterp>(fPred_maxMix_Q2, "prediction").release();
77  PredictionInterp* prediction_maxMix_Q3 = LoadFrom<PredictionInterp>(fPred_maxMix_Q3, "prediction").release();
78  PredictionInterp* prediction_maxMix_Q4 = LoadFrom<PredictionInterp>(fPred_maxMix_Q4, "prediction").release();
79 
80  fPred_Q1 -> Close();
81  fPred_Q2 -> Close();
82  fPred_Q3 -> Close();
83  fPred_Q4 -> Close();
84  fPred_maxMix_Q1 -> Close();
85  fPred_maxMix_Q2 -> Close();
86  fPred_maxMix_Q3 -> Close();
87  fPred_maxMix_Q4 -> Close();
88  delete fPred_Q1;
89  delete fPred_Q2;
90  delete fPred_Q3;
91  delete fPred_Q4;
92  delete fPred_maxMix_Q1;
93  delete fPred_maxMix_Q2;
94  delete fPred_maxMix_Q3;
95  delete fPred_maxMix_Q4;
96 
97  // And we're ready to go
98 
99  // custom osc. parameters (non-max mixing) Get the numu result parameters!
100  osc::OscCalcPMNSOpt calc_nonMax;
101  ana::ResetOscCalcToDefault(&calc_nonMax);
102  calc_nonMax.SetDmsq32(2.67e-3);
103  calc_nonMax.SetTh23(asin(sqrt(0.402))); // PRL-published non-max mixing
104 
105  osc::OscCalcPMNSOpt calc_maxMix;
106  ana::ResetOscCalcToDefault(&calc_maxMix);
107  calc_maxMix.SetDmsq32(2.5e-3);
108  calc_maxMix.SetTh23(asin(sqrt(0.5)));
109 
110  std::vector<const ISyst*> systs = getAllDirectNumuSysts2017();
111  systs =
112  {
117  //&kBeamAllTransport,
118  //kMECq0ShapeSyst2017,
121  &kRPACCQEEnhSyst2017, // new RPA error (1-sided more enhancement)
122  &kRPACCQESuppSyst2017, // new RPA error (1-sided more suppression)
124  &kMuEScaleSyst2017, // correlated muon E scale applied to muon track length
125  &kRelMuEScaleSyst2017, // uncorrelated muon E scale applied to muon track length
126  &kDirectRelHadEScaleSyst2017, // directly applied 5% hadronic E uncorrelated syst
127  &kDirectHadEScaleSyst2017, // will be provided elsewhere soon.
128  kMAQEGenieReducedSyst2017, // replacing fReweightNormCCQE and fReweightMaCCQEshape from SA. Reduced M_A uncertainty (as compared to GENIE M_A knob)
129  };
130  systs.push_back(&kNumuSummedSmallGENIESyst); // horrible and needs to be re-looked at for TA
131  std::set<rwgt::ReweightLabel_t> largeGENIEParameters;
132  largeGENIEParameters.insert(rwgt::fReweightMaNCEL);
133  largeGENIEParameters.insert(rwgt::fReweightMaCCRES);
134  largeGENIEParameters.insert(rwgt::fReweightMvCCRES);
135  largeGENIEParameters.insert(rwgt::fReweightMaNCRES);
136  largeGENIEParameters.insert(rwgt::fReweightMvNCRES);
137  largeGENIEParameters.insert(rwgt::fReweightCCQEPauliSupViaKF);
138  for (const auto & parameter:largeGENIEParameters)
139  systs.push_back(GetGenieKnobSyst(parameter));
140 
141  systs.push_back(&kCosmicBkgScaleSyst);
142 
143  std::vector<const IFitVar*> oscpar = getNumuMarginalisedOscParam();
144  std::vector<const IFitVar*> margPars_ssth23 = getNumuMarginalisedOscParam();
145  margPars_ssth23.push_back(&kFitDmSq32Scaled);
146  std::vector<const IFitVar*> margPars_dmsq32 = getNumuMarginalisedOscParam();
147  margPars_dmsq32.push_back(&kFitSinSqTheta23);
148 
149 
150  Spectrum fakePred_Q1 = prediction_Q1->Predict(&calc_nonMax);
151  Spectrum fakeData_Q1 = fakePred_Q1.FakeData(kAna2017POT);
152  Spectrum fakePred_Q2 = prediction_Q2->Predict(&calc_nonMax);
153  Spectrum fakeData_Q2 = fakePred_Q2.FakeData(kAna2017POT);
154  Spectrum fakePred_Q3 = prediction_Q3->Predict(&calc_nonMax);
155  Spectrum fakeData_Q3 = fakePred_Q3.FakeData(kAna2017POT);
156  Spectrum fakePred_Q4 = prediction_Q4->Predict(&calc_nonMax);
157  Spectrum fakeData_Q4 = fakePred_Q4.FakeData(kAna2017POT);
158 
159  //calc_maxMix
160  Spectrum fakePred_maxMix_Q1 = prediction_maxMix_Q1->Predict(&calc_maxMix);
161  Spectrum fakeData_maxMix_Q1 = fakePred_maxMix_Q1.FakeData(kAna2017POT);
162  Spectrum fakePred_maxMix_Q2 = prediction_maxMix_Q2->Predict(&calc_maxMix);
163  Spectrum fakeData_maxMix_Q2 = fakePred_maxMix_Q2.FakeData(kAna2017POT);
164  Spectrum fakePred_maxMix_Q3 = prediction_maxMix_Q3->Predict(&calc_maxMix);
165  Spectrum fakeData_maxMix_Q3 = fakePred_maxMix_Q3.FakeData(kAna2017POT);
166  Spectrum fakePred_maxMix_Q4 = prediction_maxMix_Q4->Predict(&calc_maxMix);
167  Spectrum fakeData_maxMix_Q4 = fakePred_maxMix_Q4.FakeData(kAna2017POT);
168 
169  // Get the cosmic distributions:
170  TFile* fCosmic = new TFile("/nova/app/users/bays/dev_xchecks/work/done5/cos_bkg_hists_double.root");
171  TH1D* hCos_Q1 = (TH1D*) fCosmic -> Get("q1all");
172  TH1D* hCos_Q2 = (TH1D*) fCosmic -> Get("q2all");
173  TH1D* hCos_Q3 = (TH1D*) fCosmic -> Get("q3all");
174  TH1D* hCos_Q4 = (TH1D*) fCosmic -> Get("q4all");
175  fCosmic -> Close();
176  delete fCosmic;
177 
178  // Add the cosmic background to the fake data:
179  Spectrum specCos_Q1(hCos_Q1, kAna2017POT, kAna2017Livetime);
180  fakeData_Q1 += specCos_Q1;
181  fakeData_maxMix_Q1 += specCos_Q1;
182  Spectrum specCos_Q2(hCos_Q2, kAna2017POT, kAna2017Livetime);
183  fakeData_Q2 += specCos_Q2;
184  fakeData_maxMix_Q2 += specCos_Q2;
185  Spectrum specCos_Q3(hCos_Q3, kAna2017POT, kAna2017Livetime);
186  fakeData_Q3 += specCos_Q3;
187  fakeData_maxMix_Q3 += specCos_Q3;
188  Spectrum specCos_Q4(hCos_Q4, kAna2017POT, kAna2017Livetime);
189  fakeData_Q4 += specCos_Q4;
190  fakeData_maxMix_Q4 += specCos_Q4;
191 
192  // Experiments with constraints
193  // Temporary rough estimate for uncertainty on cosmic events (this will not have a visible effect on the fit)
194 
195  // SA BF
196  SingleSampleExperiment fakeExpt_Q1(prediction_Q1, fakeData_Q1, hCos_Q1, 0.1);
197  SingleSampleExperiment fakeExpt_Q2(prediction_Q2, fakeData_Q2, hCos_Q2, 0.1);
198  SingleSampleExperiment fakeExpt_Q3(prediction_Q3, fakeData_Q3, hCos_Q3, 0.1);
199  SingleSampleExperiment fakeExpt_Q4(prediction_Q4, fakeData_Q4, hCos_Q4, 0.1);
200 
201  std::vector<const IExperiment*> vConstraints;
202  vConstraints.push_back(&fakeExpt_Q1);
203  vConstraints.push_back(&fakeExpt_Q2);
204  vConstraints.push_back(&fakeExpt_Q3);
205  vConstraints.push_back(&fakeExpt_Q4);
206  vConstraints.push_back(&kSolarConstraintsPDG2017);
207  vConstraints.push_back(WorldReactorConstraint2017());
208  MultiExperiment ExptWconstr(vConstraints);
209 
210  // 2D Fits
211  FrequentistSurface* surface = new FrequentistSurface( &ExptWconstr, &calc_nonMax,
212  &kFitSinSqTheta23, 40, 0.3, 0.7,
213  &kFitDmSq32Scaled, 50, 2, 3.0, oscpar, systs);
214  FrequentistSurface* surface_stats = new FrequentistSurface( &ExptWconstr, &calc_nonMax,
215  &kFitSinSqTheta23, 40, 0.3, 0.7,
216  &kFitDmSq32Scaled, 50, 2, 3.0, oscpar, {});
217 
218 
219  // max mix
220  SingleSampleExperiment fakeExpt_maxMix_Q1(prediction_maxMix_Q1, fakeData_maxMix_Q1, hCos_Q1, 0.1);
221  SingleSampleExperiment fakeExpt_maxMix_Q2(prediction_maxMix_Q2, fakeData_maxMix_Q2, hCos_Q2, 0.1);
222  SingleSampleExperiment fakeExpt_maxMix_Q3(prediction_maxMix_Q3, fakeData_maxMix_Q3, hCos_Q3, 0.1);
223  SingleSampleExperiment fakeExpt_maxMix_Q4(prediction_maxMix_Q4, fakeData_maxMix_Q4, hCos_Q4, 0.1);
224 
225  std::vector<const IExperiment*> vConstraints_maxMix;
226  vConstraints_maxMix.push_back(&fakeExpt_maxMix_Q1);
227  vConstraints_maxMix.push_back(&fakeExpt_maxMix_Q2);
228  vConstraints_maxMix.push_back(&fakeExpt_maxMix_Q3);
229  vConstraints_maxMix.push_back(&fakeExpt_maxMix_Q4);
230  vConstraints_maxMix.push_back(&kSolarConstraintsPDG2017);
231  vConstraints_maxMix.push_back(WorldReactorConstraint2017());
232  MultiExperiment ExptWconstr_maxMix(vConstraints_maxMix);
233 
234  // 2D Fits
235  FrequentistSurface* surface_maxMix = new FrequentistSurface( &ExptWconstr_maxMix, &calc_maxMix,
236  &kFitSinSqTheta23, 40, 0.3, 0.7,
237  &kFitDmSq32Scaled, 50, 2, 3.0,
238  oscpar, systs);
239 
240  FrequentistSurface* surface_maxMix_stats = new FrequentistSurface( &ExptWconstr_maxMix, &calc_maxMix,
241  &kFitSinSqTheta23, 40, 0.3, 0.7,
242  &kFitDmSq32Scaled, 50, 2, 3.0,
243  oscpar, {});
244 
245 
246 
247  // Plot
248  TCanvas* c = new TCanvas();
249  c->SetLeftMargin(0.12);
250  c->SetBottomMargin(0.15);
251 
252  surface->DrawBestFit(kBlue);
253  surface->DrawContour( Gaussian90Percent2D(*surface), kSolid, kBlue);
254 
255  c->SaveAs("Sensitivity.root");
256  c->SaveAs("Sensitivity.pdf");
257 
258  c->Clear();
259 
260  // Get 1D marginalised plots
261  TH1* profileSinsq = Profile(&ExptWconstr, &calc_nonMax,
262  &kFitSinSqTheta23, 40, 0.3, 0.7, -1,
263  margPars_ssth23, systs );
264  TH1* profileSinsq_stat = Profile(&ExptWconstr, &calc_nonMax,
265  &kFitSinSqTheta23, 40, 0.3, 0.7, -1,
266  margPars_ssth23, {} );
267 
268  profileSinsq->SetLineColor(kBlue);
269  profileSinsq->Draw("hist");
270 
271  std::cout << "Max-mix rejection: " << profileSinsq->Interpolate(0.5) << std::endl;
272 
273  c->SaveAs("Ssth23_1D.root");
274  c->SaveAs("Ssth23_1D.pdf");
275 
276  c->Clear();
277 
278  TH1* profileDmsq = Profile(&ExptWconstr, &calc_nonMax,
279  &kFitDmSq32Scaled, 35, 2.3, 3.0, -1,
280  margPars_dmsq32, systs );
281 
282  TH1* profileDmsq_stat = Profile(&ExptWconstr, &calc_nonMax,
283  &kFitDmSq32Scaled, 35, 2.3, 3.0, -1,
284  margPars_dmsq32, {} );
285 
286  profileDmsq->SetLineColor(kBlue);
287  profileDmsq->Draw("hist");
288 
289  c->SaveAs("Dmsq32_1D.root");
290  c->SaveAs("Dmsq32_1D.pdf");
291  c->Close();
292 
293  // Save surface and 1D marg. hists
294  TFile* fSurf = new TFile("SensitivityAna2017.root","RECREATE");
295  surface -> SaveTo(fSurf, "surface");
296  surface_stats -> SaveTo(fSurf, "surface_stats_cosmics");
297  surface_maxMix -> SaveTo(fSurf, "surface_maxMix");
298  surface_maxMix_stats -> SaveTo(fSurf, "surface_maxMix_stats");
299  profileDmsq -> Write("profileDmsq");
300  profileSinsq -> Write("profileSinsq");
301  profileDmsq_stat -> Write("profileDmsq_stat");
302  profileSinsq_stat -> Write("profileSinsq_stat");
303  fSurf -> Close();
304  delete fSurf;
305 
306 
307 } // End of function
BeamSyst * GetPPFXPrincipals(int PCIdx)
Definition: BeamSysts.cxx:63
Implements systematic errors by interpolation between shifted templates.
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
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
void getSensitivity()
T sqrt(T number)
Definition: d0nt_math.hpp:156
const NOvARwgtSyst kMAQEGenieReducedSyst2017(genie::rew::GSyst::AsString(genie::rew::EGSyst(rwgt::fReweightMaCCQE))++"_reduced", genie::rew::GSyst::AsString(genie::rew::EGSyst(rwgt::fReweightMaCCQE))++"_reduced", novarwgt::kMAQEGenieReducedSyst2017)
2017 &#39;reduced&#39; M_A^QE shift. See documentation in NOvARwgt.
Definition: XSecSysts.h:47
const double kAna2017POT
Definition: Exposures.h:174
TGraph * Profile(const IExperiment *expt, osc::IOscCalcAdjustable *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
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
const CosmicBkgScaleSyst kCosmicBkgScaleSyst
const NOvARwgtSyst kMECEnuShapeSyst2017("MECEnuShape","MEC E_{#nu} shape", novarwgt::kMECEnuShapeSyst2017)
Definition: MECSysts.h:43
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const MuEScaleSyst2017 kMuEScaleSyst2017(0.0074, 0.0012)
Log-likelihood scan across two parameters.
const NOvARwgtSyst kRPACCQEEnhSyst2017("RPAShapeenh2017","RPA shape: higher-Q^{2} enhancement (2017)", novarwgt::kRPACCQEEnhSyst2017)
Definition: RPASysts.h:11
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
Definition: ISurface.cxx:32
const double kAna2017Livetime
Definition: Exposures.h:200
const RelMuEScaleSyst2017 kRelMuEScaleSyst2017(0.0045, 10.5)
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
Definition: Spectrum.cxx:349
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:44
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
const NOvARwgtSyst kRPACCQESuppSyst2017("RPAShapesupp2017","RPA shape: low-Q^{2} suppression (2017)", novarwgt::kRPACCQESuppSyst2017)
Definition: RPASysts.h:12
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
const DirectHadEScaleSyst2017 kDirectHadEScaleSyst2017(0.05)
void SetTh23(const T &th23) override
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
c1 Close()
OStream cout
Definition: OStream.cxx:6
void SetDmsq32(const T &dmsq32) override
Spectrum Predict(osc::IOscCalc *calc) const override
Combine multiple component experiments.
const DirectRelHadEScaleSyst2017 kDirectRelHadEScaleSyst2017(0.05)
const SolarConstraints kSolarConstraintsPDG2017(7.53e-5, 0.18e-5, 0.851, 0.020)
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
const NumuSummedSmallGENIESyst kNumuSummedSmallGENIESyst(FindCAFAnaDir()+"/data/numu/summed_genie_systs_nd.root", FindCAFAnaDir()+"/data/numu/summed_genie_systs_fd.root")
Definition: NumuSysts.h:62
const NOvARwgtSyst kRPARESSyst2017("RPAShapeRES2017","RPA shape: resonance events", novarwgt::kRPARESSyst2017)
Definition: RPASysts.h:20
std::vector< const ISyst * > getAllDirectNumuSysts2017()
Get a vector of all the numu group systs.
Definition: NumuSysts.cxx:168
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
const NOvARwgtSyst * GetGenieKnobSyst(rwgt::ReweightLabel_t knobIdx, std::string altName, std::string altLabel)
Convenience function to get a GENIE knob syst. (Allows using the GENIE knob name & description as the...
Definition: XSecSysts.cxx:119
enum BeamMode kBlue
const NOvARwgtSyst kMECInitStateNPFracSyst2017("MECInitStateNPFrac","MEC initial state np fraction", novarwgt::kMECInitStateNPFracSyst2017)
Definition: MECSysts.h:42
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
gm Write()
std::vector< const IFitVar * > getNumuMarginalisedOscParam()
Definition: FitVarsNumu.cxx:9