Functions | Variables
caf_numu_sensitivity_no_tau.C File Reference
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Fit.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Surface.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Core/ReweightableSpectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Vars/FitVars.h"
#include "CAFAna/Vars/FitVarsNumu.h"
#include "CAFAna/Vars/Vars.h"
#include "OscLib/func/OscCalculatorPMNSOpt.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TGraph.h"
#include "TList.h"
#include "TObject.h"
#include "TStyle.h"

Go to the source code of this file.

Functions

TH1D * GetXMarginalisation (TH2 *h2)
 
TH1D * GetYMarginalisation (TH2 *h2)
 
void caf_numu_sensitivity_no_tau (std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kOutputFileName)
 

Variables

const double pot = 6.e20
 
const Cut kIsDytmanMEC ({"mc.nnu","mc.nu.mode"}, [](const caf::StandardRecord *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return(sr->mc.nu[0].mode==caf::kMEC);})
 
const Var kSlcMaxX ({"slc.boxmax.fX"}, [](const caf::StandardRecord *sr){return sr->slc.boxmax.X()/100;})
 
const Var kSlcMaxY ({"slc.boxmax.fY"}, [](const caf::StandardRecord *sr){return sr->slc.boxmax.Y()/100;})
 
const Var kSlcMaxZ ({"slc.boxmax.fZ"}, [](const caf::StandardRecord *sr){return sr->slc.boxmax.Z()/100;})
 

Function Documentation

void caf_numu_sensitivity_no_tau ( std::string  kInputNonSwapFileName,
std::string  kInputSwapFileName,
std::string  kOutputFileName 
)

Definition at line 65 of file caf_numu_sensitivity_no_tau.C.

References plot_validation_datamc::c, calc, om::cout, allTimeWatchdog::endl, MakeMiniprodValidationCuts::f, ana::Spectrum::FakeData(), ana::Gaussian68Percent2D(), ana::Gaussian90Percent2D(), GetXMarginalisation(), GetYMarginalisation(), ana::SpectrumLoader::Go(), ana::kCCE, ana::kFitDmSq32Scaled, ana::kFitSinSqTheta23, ana::kIsDytmanMEC, ana::kNumuFD, ana::kStandardSpillCuts, demo4::loaderSwap, next(), pot, ana::PredictionExtrap::Predict(), ana::ResetOscCalcToDefault(), ana::SpectrumLoaderBase::SetSpillCut(), ana::Binning::Simple(), and demo5::surf.

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
const double pot
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
const Cut kIsDytmanMEC({"mc.nnu","mc.nu.mode"}, [](const caf::StandardRecord *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return(sr->mc.nu[0].mode==caf::kMEC);})
const Cut kNumuFD
Definition: NumuCuts.h:53
osc::OscCalcDumb calc
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
Spectrum Predict(osc::IOscCalc *calc) const override
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:402
const Var kCCE
Definition: NumuVars.h:21
TH1D * GetXMarginalisation(TH2 *h2)
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
OStream cout
Definition: OStream.cxx:6
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
TH1D * GetYMarginalisation(TH2 *h2)
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:46
Prediction that just uses FD MC, with no extrapolation.
loaderSwap
Definition: demo4.py:9
void next()
Definition: show_event.C:84
Compare a single data spectrum to the MC + cosmics expectation.
surf
Definition: demo5.py:35
TH1D* GetXMarginalisation ( TH2 *  h2)

Definition at line 177 of file caf_numu_sensitivity_no_tau.C.

References febshutoff_auto::val, Xmax, Xmin, Ymax, and Ymin.

Referenced by caf_numu_sensitivity_no_tau().

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 }
int Ymin
int Ymax
int Xmax
TH1F * h2
Definition: plot.C:45
int Xmin
TH1D* GetYMarginalisation ( TH2 *  h2)

Definition at line 201 of file caf_numu_sensitivity_no_tau.C.

References febshutoff_auto::val, Xmax, Xmin, Ymax, and Ymin.

Referenced by caf_numu_sensitivity_no_tau().

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 }
int Ymin
int Ymax
int Xmax
TH1F * h2
Definition: plot.C:45
int Xmin

Variable Documentation

const Cut kIsDytmanMEC({"mc.nnu","mc.nu.mode"},[](const caf::StandardRecord *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return(sr->mc.nu[0].mode==caf::kMEC);})
const Var kSlcMaxX({"slc.boxmax.fX"},[](const caf::StandardRecord *sr){return sr->slc.boxmax.X()/100;})
const Var kSlcMaxY({"slc.boxmax.fY"},[](const caf::StandardRecord *sr){return sr->slc.boxmax.Y()/100;})
const Var kSlcMaxZ({"slc.boxmax.fZ"},[](const caf::StandardRecord *sr){return sr->slc.boxmax.Z()/100;})
const double pot = 6.e20

Definition at line 36 of file caf_numu_sensitivity_no_tau.C.

Referenced by caf_numu_sensitivity_no_tau().