caf_numu_sensitivity_no_tau.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 void caf_numu_sensitivity_no_tau(std::string kInputNonSwapFileName, std::string kInputSwapFileName, 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"
14 #include "CAFAna/Cuts/TruthCuts.h"
17 #include "CAFAna/Core/Spectrum.h"
20 #include "CAFAna/Vars/FitVars.h"
21 #include "CAFAna/Vars/FitVarsNumu.h"
22 #include "CAFAna/Vars/Vars.h"
23 #include "OscLib/func/OscCalculatorPMNSOpt.h"
24 
25 using namespace ana;
26 
27 #include "TCanvas.h"
28 #include "TFile.h"
29 #include "TH1.h"
30 #include "TH2.h"
31 #include "TGraph.h"
32 #include "TList.h"
33 #include "TObject.h"
34 #include "TStyle.h"
35 
36 const double pot = 6.e20;
37 
38 TH1D* GetXMarginalisation(TH2* h2);
39 TH1D* GetYMarginalisation(TH2* h2);
40 
41 const Cut kIsDytmanMEC({"mc.nnu", "mc.nu.mode"},
42  [](const caf::StandardRecord* sr)
43  {
44  if(sr->mc.nnu == 0) return false;
45  assert(sr->mc.nnu == 1);
46  return (sr->mc.nu[0].mode == caf::kMEC);
47  });
48 
49 const Var kSlcMaxX({"slc.boxmax.fX"},
50  [](const caf::StandardRecord *sr)
51  {
52  return sr->slc.boxmax.X()/100;
53  });
54 const Var kSlcMaxY({"slc.boxmax.fY"},
55  [](const caf::StandardRecord *sr)
56  {
57  return sr->slc.boxmax.Y()/100;
58  });
59 const Var kSlcMaxZ({"slc.boxmax.fZ"},
60  [](const caf::StandardRecord *sr)
61  {
62  return sr->slc.boxmax.Z()/100;
63  });
64 
65 void caf_numu_sensitivity_no_tau(std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kOutputFileName)
66 {
67  std::cout << "\nrun : ---- Running CAF NuMu sensitivity (no tau).\n";
68  std::cout << "run : input non-swap file name: " << kInputNonSwapFileName << std::endl;
69  std::cout << "run : input swap file name: " << kInputSwapFileName << 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 
79  loaderNonswap.SetSpillCut(kStandardSpillCuts);
80  loaderSwap.SetSpillCut(kStandardSpillCuts);
81 
82  // Taus are not required, you can use kNullLoader for fast tests or if they're not available
83  PredictionNoExtrap* prediction = new PredictionNoExtrap(loaderNonswap, loaderSwap,
84  // "Reconstructed neutrino Energy [GeV]", Binning::Simple(20,0,5), kCCE, kNumuFD);
85  "Reconstructed neutrino Energy [GeV]", Binning::Simple(20,0,5), kCCE, kNumuFD && (!kIsDytmanMEC));
86 
87  std::cout << "\nrun : --- run loaders.\n";
88  loaderNonswap.Go();
89  loaderSwap.Go();
90  std::cout << "\nrun : --- done.\n";
91 
92  std::cout << "\nrun : --- Make surfaces.\n";
93  // Create fake prediction and data, normalised to the chosen POT
94  Spectrum fakePred = prediction->Predict(&calc);
95  Spectrum fakeData = fakePred.FakeData(pot);
96 
97  // Experiment, i.e., create set of prediction vs data
98  SingleSampleExperiment fakeExpt(prediction, fakeData);
99 
100  // Fill the chi-2 (LL) 2-D histogram
101  Surface* sensitivity = new Surface( &fakeExpt, &calc,
102  &kFitSinSqTheta23, 40, 0.3, 0.7,
103  &kFitDmSq32Scaled, 50, 2.0, 3.0 );
104 
105  std::cout << "\nrun : --- save output.\n";
106  TFile* f = new TFile(kOutputFileName.c_str(),"RECREATE");
107 
108  // Save the surface
109  TH2* surf = sensitivity->ToTH2();
110  surf->SetName("reco-sensitivity_surface-numu_selection");
111  surf->Write();
112 
113  // Save the best fit
114  TGraph* bestFit = new TGraph;
115  bestFit->SetPoint(0, sensitivity->fMinX, sensitivity->fMinY);
116  bestFit->SetName("reco-sensitivity_bestfit-numu_selection");
117  bestFit->GetXaxis()->SetLimits(0.3, 0.7);
118  bestFit->GetYaxis()->SetLimits(2.0, 3.0);
119  bestFit->Write();
120 
121  // Save 1D marginalisations
122  TH1D* margX = GetXMarginalisation(surf);
123  TH1D* margY = GetYMarginalisation(surf);
124  margX->SetName("reco-sensitivity_1DmarginX-numu_selection");
125  margY->SetName("reco-sensitivity_1DmarginY-numu_selection");
126  margX->Write();
127  margY->Write();
128 
129  // Now different contours
130  TCanvas* c = new TCanvas("c","c");
131  c->SetLeftMargin(0.12);
132  c->SetBottomMargin(0.15);
133 
134  TList *ObjList = c->GetListOfPrimitives();
135  TObject* obj;
136  TIter next(ObjList);
137 
138  TGraph* g68 = new TGraph;
139  TGraph* g90 = new TGraph;
140 
141  // Draw the 68% contour
142  sensitivity->DrawContour( Gaussian68Percent2D(*sensitivity), kDashed, kRed);
143 
144  while(obj=(TObject *)next())
145  {
146  if(!obj->InheritsFrom(TGraph::Class())) continue;
147  g68 = (TGraph*) obj->Clone("reco-sensitivity_68pc_contour-numu_selection");
148  }
149 
150  g68->Write();
151 
152  // Draw the 90% contour
153  c->Clear();
154  sensitivity->DrawContour( Gaussian90Percent2D(*sensitivity), kSolid, kRed);
155 
156  next.Begin();
157 
158  while(obj=(TObject *)next())
159  {
160  if(!obj->InheritsFrom(TGraph::Class())) continue;
161  g90 = (TGraph*) obj->Clone("reco-sensitivity_90pc_contour-numu_selection");
162  }
163 
164  g90->Write();
165 
166  TH1D* TotalPOT = new TH1D("meta-TotalPOT","TotalPOT", 1, 0, 1);
167  TotalPOT->SetBinContent(1, pot);
168  TotalPOT->Write();
169 
170  f->Write();
171  f->Close();
172 
173  std::cout << "All plots finished!" << std::endl;
174 
175 } // End of function
176 
178 {
179  int nY = h2->GetYaxis()->GetNbins();
180  int nX = h2->GetXaxis()->GetNbins();
181 
182  double Xmin = h2->GetXaxis()->GetXmin();
183  double Xmax = h2->GetXaxis()->GetXmax();
184  double Ymin = h2->GetYaxis()->GetXmin();
185  double Ymax = h2->GetYaxis()->GetXmax();
186 
187  TH1D *h1Dim = new TH1D("","",nX,Xmin,Xmax);
188 
189  for(int binx = 1; binx <= nX; binx++) {
190  double minchi = 999999.;
191  for(int biny = 1; biny <= nY; biny++) {
192  double val = h2->GetBinContent(binx,biny);
193  if(val < minchi) minchi = val;
194  }
195  h1Dim->SetBinContent(binx,minchi);
196  }
197 
198  return h1Dim;
199 }
200 
202 {
203  int nY = h2->GetYaxis()->GetNbins();
204  int nX = h2->GetXaxis()->GetNbins();
205 
206  double Xmin = h2->GetXaxis()->GetXmin();
207  double Xmax = h2->GetXaxis()->GetXmax();
208  double Ymin = h2->GetYaxis()->GetXmin();
209  double Ymax = h2->GetYaxis()->GetXmax();
210 
211  TH1D *h1Dim = new TH1D("","",nY,Ymin,Ymax);
212 
213  for(int biny = 1; biny <= nY; biny++) {
214  double minchi = 999999.;
215  for(int binx = 1; binx <= nX; binx++) {
216  double val = h2->GetBinContent(binx,biny);
217  if(val < minchi) minchi = val;
218  }
219  h1Dim->SetBinContent(biny,minchi);
220  }
221 
222  return h1Dim;
223 }
224 
225 #endif
const double pot
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
Spectrum Predict(osc::IOscCalc *calc) const override
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
Definition: Spectrum.cxx:349
const Var kCCE
Definition: NumuVars.h:21
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
TH1D * GetXMarginalisation(TH2 *h2)
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
void caf_numu_sensitivity_no_tau(std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kOutputFileName)
The StandardRecord is the primary top-level object in the Common Analysis File trees.
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)
TH1D * GetYMarginalisation(TH2 *h2)
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