caf_numu_sensitivity.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 void caf_numu_sensitivity(std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kInputTauFileName, std::string kOutputFileName){
3  std::cout << "Sorry, you must run in compiled mode" << std::endl;
4 }
5 #else
6 
8 #include "CAFAna/Analysis/Fit.h"
10 #include "CAFAna/Analysis/Surface.h"
11 #include "CAFAna/Cuts/Cuts.h"
12 #include "CAFAna/Cuts/NumuCuts.h"
13 #include "CAFAna/Cuts/SpillCuts.h"
16 #include "CAFAna/Core/Spectrum.h"
19 #include "CAFAna/Vars/FitVars.h"
20 #include "CAFAna/Vars/FitVarsNumu.h"
21 #include "CAFAna/Vars/Vars.h"
22 #include "OscLib/func/OscCalculatorPMNSOpt.h"
23 
24 using namespace ana;
25 
26 #include "TCanvas.h"
27 #include "TFile.h"
28 #include "TH1.h"
29 #include "TH2.h"
30 #include "TGraph.h"
31 #include "TList.h"
32 #include "TObject.h"
33 #include "TStyle.h"
34 
35 const double pot = 6.e20;
36 
37 TH1D* GetXMarginalisation(TH2* h2);
38 TH1D* GetYMarginalisation(TH2* h2);
39 
40 const Cut kIsDytmanMEC({"mc.nnu", "mc.nu.mode"},
41  [](const caf::StandardRecord* sr)
42  {
43  if(sr->mc.nnu == 0) return false;
44  assert(sr->mc.nnu == 1);
45  return (sr->mc.nu[0].mode == caf::kMEC);
46  });
47 
48 const Var kSlcMaxX({"slc.boxmax.fX"},
49  [](const caf::StandardRecord *sr)
50  {
51  return sr->slc.boxmax.X()/100;
52  });
53 const Var kSlcMaxY({"slc.boxmax.fY"},
54  [](const caf::StandardRecord *sr)
55  {
56  return sr->slc.boxmax.Y()/100;
57  });
58 const Var kSlcMaxZ({"slc.boxmax.fZ"},
59  [](const caf::StandardRecord *sr)
60  {
61  return sr->slc.boxmax.Z()/100;
62  });
63 
64 void caf_numu_sensitivity(std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kInputTauFileName, std::string kOutputFileName)
65 {
66  std::cout << "\nrun : ---- Running CAF NuMu sensitivity.\n";
67  std::cout << "run : input non-swap file name: " << kInputNonSwapFileName << std::endl;
68  std::cout << "run : input swap file name: " << kInputSwapFileName << std::endl;
69  std::cout << "run : input tau file name: " << kInputTauFileName << std::endl;
70  std::cout << "run : output file name: " << kOutputFileName << std::endl;
71 
72 
73  osc::OscCalculatorPMNSOpt calc;
74  ana::ResetOscCalcToDefault(&calc); // A calculator set to default values
75 
76  SpectrumLoader loaderNonswap(kInputNonSwapFileName);
77  SpectrumLoader loaderSwap(kInputSwapFileName);
78  SpectrumLoader loaderTau(kInputTauFileName);
79 
80  loaderNonswap.SetSpillCut(kStandardSpillCuts);
81  loaderSwap.SetSpillCut(kStandardSpillCuts);
83 
84  // Taus are not required, you can use kNullLoader for fast tests or if they're not available
85  PredictionNoExtrap* prediction = new PredictionNoExtrap(loaderNonswap, loaderSwap, loaderTau,
86  "Reconstructed neutrino Energy [GeV]", Binning::Simple(20,0,5), kCCE, kNumuFD && (!kIsDytmanMEC));
87 
88  std::cout << "\nrun : --- run loaders.\n";
89  loaderNonswap.Go();
90  loaderSwap.Go();
91  loaderTau.Go();
92  std::cout << "\nrun : --- done.\n";
93 
94  std::cout << "\nrun : --- Make surfaces.\n";
95  // Create fake prediction and data, normalised to the chosen POT
96  Spectrum fakePred = prediction->Predict(&calc);
97  Spectrum fakeData = fakePred.FakeData(pot);
98 
99  // Experiment, i.e., create set of prediction vs data
100  SingleSampleExperiment fakeExpt(prediction, fakeData);
101 
102  // Fill the chi-2 (LL) 2-D histogram
103  Surface* sensitivity = new Surface( &fakeExpt, &calc,
104  &kFitSinSqTheta23, 40, 0.3, 0.7,
105  &kFitDmSq32Scaled, 50, 2.0, 3.0 );
106 
107  std::cout << "\nrun : --- save output.\n";
108  TFile* f = new TFile(kOutputFileName.c_str(),"RECREATE");
109 
110  // Save the surface
111  TH2* surf = sensitivity->ToTH2();
112  surf->SetName("reco-sensitivity_surface-numu_selection");
113  surf->Write();
114 
115  // Save the best fit
116  TGraph* bestFit = new TGraph;
117  bestFit->SetPoint(0, sensitivity->fMinX, sensitivity->fMinY);
118  bestFit->SetName("reco-sensitivity_bestfit-numu_selection");
119  bestFit->GetXaxis()->SetLimits(0.3, 0.7);
120  bestFit->GetYaxis()->SetLimits(2.0, 3.0);
121  bestFit->Write();
122 
123  // Save 1D marginalisations
124  TH1D* margX = GetXMarginalisation(surf);
125  TH1D* margY = GetYMarginalisation(surf);
126  margX->SetName("reco-sensitivity_1DmarginX-numu_selection");
127  margY->SetName("reco-sensitivity_1DmarginY-numu_selection");
128  margX->Write();
129  margY->Write();
130 
131  // Now different contours
132  TCanvas* c = new TCanvas("c","c");
133  c->SetLeftMargin(0.12);
134  c->SetBottomMargin(0.15);
135 
136  TList *ObjList = c->GetListOfPrimitives();
137  TObject* obj;
138  TIter next(ObjList);
139 
140  TGraph* g68 = new TGraph;
141  TGraph* g90 = new TGraph;
142 
143  // Draw the 68% contour
144  sensitivity->DrawContour( Gaussian68Percent2D(*sensitivity), kDashed, kRed);
145 
146  while(obj=(TObject *)next())
147  {
148  if(!obj->InheritsFrom(TGraph::Class())) continue;
149  g68 = (TGraph*) obj->Clone("reco-sensitivity_68pc_contour-numu_selection");
150  }
151 
152  g68->Write();
153 
154  // Draw the 90% contour
155  c->Clear();
156  sensitivity->DrawContour( Gaussian90Percent2D(*sensitivity), kSolid, kRed);
157 
158  next.Begin();
159 
160  while(obj=(TObject *)next())
161  {
162  if(!obj->InheritsFrom(TGraph::Class())) continue;
163  g90 = (TGraph*) obj->Clone("reco-sensitivity_90pc_contour-numu_selection");
164  }
165 
166  g90->Write();
167 
168  TH1D* TotalPOT = new TH1D("meta-TotalPOT","TotalPOT", 1, 0, 1);
169  TotalPOT->SetBinContent(1, pot);
170  TotalPOT->Write();
171 
172  f->Write();
173  f->Close();
174 
175  std::cout << "All plots finished!" << std::endl;
176 
177 } // End of function
178 
180 {
181  int nY = h2->GetYaxis()->GetNbins();
182  int nX = h2->GetXaxis()->GetNbins();
183 
184  double Xmin = h2->GetXaxis()->GetXmin();
185  double Xmax = h2->GetXaxis()->GetXmax();
186  double Ymin = h2->GetYaxis()->GetXmin();
187  double Ymax = h2->GetYaxis()->GetXmax();
188 
189  TH1D *h1Dim = new TH1D("","",nX,Xmin,Xmax);
190 
191  for(int binx = 1; binx <= nX; binx++) {
192  double minchi = 999999.;
193  for(int biny = 1; biny <= nY; biny++) {
194  double val = h2->GetBinContent(binx,biny);
195  if(val < minchi) minchi = val;
196  }
197  h1Dim->SetBinContent(binx,minchi);
198  }
199 
200  return h1Dim;
201 }
202 
204 {
205  int nY = h2->GetYaxis()->GetNbins();
206  int nX = h2->GetXaxis()->GetNbins();
207 
208  double Xmin = h2->GetXaxis()->GetXmin();
209  double Xmax = h2->GetXaxis()->GetXmax();
210  double Ymin = h2->GetYaxis()->GetXmin();
211  double Ymax = h2->GetYaxis()->GetXmax();
212 
213  TH1D *h1Dim = new TH1D("","",nY,Ymin,Ymax);
214 
215  for(int biny = 1; biny <= nY; biny++) {
216  double minchi = 999999.;
217  for(int binx = 1; binx <= nX; binx++) {
218  double val = h2->GetBinContent(binx,biny);
219  if(val < minchi) minchi = val;
220  }
221  h1Dim->SetBinContent(biny,minchi);
222  }
223 
224  return h1Dim;
225 }
226 
227 #endif
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
const Var kSlcMaxZ([](const caf::SRProxy *sr){return sr->slc.boxmax.Z()/100.;})
Definition: NumuVars.h:93
const Var kSlcMaxY([](const caf::SRProxy *sr){return sr->slc.boxmax.Y()/100.;})
Definition: NumuVars.h:92
SRVector3D boxmax
Maximum coordinates box containing all the hits [cm].
Definition: SRSlice.h:47
float Y() const
Definition: SRVector3D.h:33
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
const Cut kNumuFD
Definition: NumuCuts.h:53
void SetSpillCut(const SpillCut &cut)
osc::OscCalcDumb calc
const Var kSlcMaxX([](const caf::SRProxy *sr){return sr->slc.boxmax.X()/100.;})
Definition: NumuVars.h:91
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
int Ymin
int Ymax
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
int Xmax
float X() const
Definition: SRVector3D.h:32
TH1D * GetXMarginalisation(TH2 *h2)
Spectrum Predict(osc::IOscCalc *calc) const override
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
Definition: Spectrum.cxx:349
TH1D * GetYMarginalisation(TH2 *h2)
const Var kCCE
Definition: NumuVars.h:21
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
const double pot
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
float Z() const
Definition: SRVector3D.h:34
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
The StandardRecord is the primary top-level object in the Common Analysis File trees.
void caf_numu_sensitivity(std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kInputTauFileName, std::string kOutputFileName)
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
SRSlice slc
Slice branch: nhit, extents, time, etc.
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Prediction that just uses FD MC, with no extrapolation.
int Xmin
loaderSwap
Definition: demo4.py:9
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
void next()
Definition: show_event.C:84
std::vector< SRNeutrino > nu
implemented as a vector to maintain mc.nu structure, i.e. not a pointer, but with 0 or 1 entries...
Definition: SRTruthBranch.h:25
Compare a single data spectrum to the MC + cosmics expectation.
surf
Definition: demo5.py:35
enum BeamMode string