sterile_demo.C
Go to the documentation of this file.
3 
4 #include "CAFAna/Core/Binning.h"
6 #include "CAFAna/Cuts/Cuts.h"
9 #include "CAFAna/Core/HistAxis.h"
10 #include "CAFAna/Core/Loaders.h"
18 #include "CAFAna/Core/Spectrum.h"
19 #include "CAFAna/Vars/Vars.h"
20 
21 #include "TCanvas.h"
22 #include "TFile.h"
23 #include "TH1.h"
24 #include "TH2.h"
26 
27 using namespace ana;
28 
30 {
31  TH1::AddDirectory(0);
32 
33  // Strings to file paths
34  const std::string fnamenear = "$NOVA_PROD/mc/FA14-11-25/genie/nd/caf/*/*/*_s*0_FA*";
35  const std::string fnamefar = "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/*/*/*_nonswap_*_r*00_s??_FA*";
36  const std::string fnameswap = "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/*/*/*_swap_*_r*00_s??_FA*";
37  const std::string fnametau = "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/*/*/*_tau_*_r*00_s??_FA*";
38  const std::string fnamedata = "$NOVA_PROD/mc/FA14-11-25/genie/nd/caf/*/*/*_s*0_FA*"; // note that this is the same as the near MC!!!
39 
40  // Set up a Loaders object
41  Loaders loaders("$NOVA_PROD/mc/FA14-10-28", Loaders::kFHC);
47 
48  // Create HistAxis objects for the extrapolation
49  const HistAxis axisNC(
50  "Calorimetric Energy (GeV)", Binning::Simple(100, 0., 10.), kCaloE);
51  const HistAxis axisNumu(
52  "Numu CC Energy Estimator (GeV)", Binning::Simple(100, 0., 10.), kCCE);
53 
54  // Create Cuts for NC Selection
55  const Cut kND([](const caf::SRProxy* sr)
56  {
57  if( !sr->vtx.elastic.IsValid) { return false; }
58  return (fabs(sr->vtx.elastic.vtx.X()) < 160. &&
59  fabs(sr->vtx.elastic.vtx.Y()) < 160. &&
60  sr->vtx.elastic.vtx.Z() > 15. && sr->vtx.elastic.vtx.Z() < 1250.);
61  });
62 
63  const Cut kFD([](const caf::SRProxy* sr)
64  {
65  if(!sr->vtx.elastic.IsValid ) { return false; }
66  return (fabs(sr->vtx.elastic.vtx.X()) < 700. &&
67  fabs(sr->vtx.elastic.vtx.Y()) < 700. &&
68  sr->vtx.elastic.vtx.Z() > 15. && sr->vtx.elastic.vtx.Z() < 5800.);
69  });
70 
71  const Cut kGoodFracMIP([](const caf::SRProxy* sr)
72  {
73  float fracmip = ((sr->slc.ncalhit > 0)
74  ? sr->slc.nmiphit/float(sr->slc.ncalhit)
75  : -1);
76  return (fracmip > 0. && fracmip < 0.5);
77  });
78 
79  const Cut kIsNotMuon([](const caf::SRProxy* sr)
80  {
81  if(sr->vtx.elastic.fuzzyk.nshwlid < 1) { return true; }
82  return sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.ismuon < 1;
83  });
84 
85  // This is the important one, ultimately
86  const Cut kNCPresel = kGoodFracMIP &&
87  kElecID < 0.5 && kLEM < 0.3 && kRemID < 0.7 &&
88  kIsNotMuon;
89 
90  // Package the cuts in a slightly nicer way
91  const Cut cutNDNumu(kNumuQEND);
92  const Cut cutND(kND && kNCPresel);
93  const Cut cutFD(kFD && kNCPresel);
94 
95  // Setup an NC and a numu decomposition
96  CheatDecomp decompNC(
98  axisNC, cutND
99  );
100  CheatDecomp decompNumu(
102  axisNumu, cutNDNumu
103  );
104 
105  // Setup the extrapolation
107  loaders, decompNC, decompNumu, axisNC, axisNumu, cutFD, cutND, cutNDNumu
108  );
109  /*ModularExtrapSterile extrap = ModularExtrapSterile::TrivialExtrapNC(
110  loaders.GetLoader(caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kFluxSwap),
111  loaders.GetLoader(caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
112  loaders.GetLoader(caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kTauSwap),
113  axisNC, cutFD
114  );*/
115 
116  // Make a sterile style prediction based on the extrapolation
117  PredictionSterile pred(&extrap);
118  // Make a standard 3 flavor prediction based on the extrapolation
119  PredictionExtrap pred_3flav(&extrap);
120  // Make a prediction based on FD MC directly from the loaders
125  axisNC.label, axisNC.bins, axisNC.var, cutFD
126  );
127 
128  // Define these spectra directly for individual component comparisons
130  axisNC, cutFD && !kIsAntiNu && kIsSig);
132  axisNC, cutFD && kIsAntiNu && kIsSig);
134  axisNC, cutFD && !kIsAntiNu && kIsNumuApp);
136  axisNC, cutFD && kIsAntiNu && kIsNumuApp);
137 
138  // Load from files
139  loaders.Go();
140 
141  // Open output file
142  TFile* rootOut = new TFile("extrapsterile_demo.root", "RECREATE");
143  // Save things in output file
144  extrap.SaveTo(rootOut, "NCextrap");
145 
146  // A 'calculator' that applies no oscillations
148 
149  // Set up oscillation calculator with default values
151 
152  // Set up oscillation calculator that uses default 3 flavor parameters
154  calcst->SetNFlavors(3);
155  calcst->SetL(810);
156  calcst->SetRho(2.84);
157  calcst->SetAngle(1, 2, asin(sqrt(.846))/2);
158  calcst->SetAngle(1, 3, asin(sqrt(.084))/2);
159  calcst->SetAngle(2, 3, M_PI/4);
160  calcst->SetDm(2, 7.53e-5); //Dm21
161  calcst->SetDm(3, 2.40e-3 + 7.53e-5); //Dm31 = Dm32 + Dm21
162  calcst->SetDelta(1, 3, 0);
163 
164  double scale = 18e20;
165 
166 /* These plots show the fraction of NC events
167  that come from neutrinos of a given starting flavor/sign */
168  new TCanvas;
169  TH1* hncnue = extrap.NCNueProportion().ToTH1();
170  hncnue->SetTitle("Proportion of NC events originating from #nu_{e}");
171  hncnue->Draw("hist");
172  new TCanvas;
173  TH1* hncantinue = extrap.NCAntiNueProportion().ToTH1();
174  hncantinue->SetTitle("Proportion of NC events originating from #bar{#nu}_{e}");
175  hncantinue->Draw("hist");
176  new TCanvas;
177  TH1* hncnumu = extrap.NCNumuProportion().ToTH1();
178  hncnumu->SetTitle("Proportion of NC events originating from #nu_{#mu}");
179  hncnumu->Draw("hist");
180  new TCanvas;
181  TH1* hncantinumu = extrap.NCAntiNumuProportion().ToTH1();
182  hncantinumu->SetTitle("Proportion of NC events originating from #bar{#nu}_{#mu}");
183  hncantinumu->Draw("hist");
184 
185 /* This plot is a check that all components add to 1 */
186 /* new TCanvas;
187  TH1* hTotal = (TH1*)hncnue->Clone();
188  hTotal->Add(hncantinue);
189  hTotal->Add(hncnumu);
190  hTotal->Add(hncantinumu);
191  hTotal->SetTitle("Sum of all Proportions");
192  hTotal->Draw("hist");
193 */
194 
195 /* These plots show each extrapolated component */
196  new TCanvas;
197  extrap.OscNCComponent().ToTH2(scale)->Draw("colz");
198  new TCanvas;
199  extrap.NueSurvComponent().ToTH2(scale)->Draw("colz");
200  new TCanvas;
201  extrap.AntiNueSurvComponent().ToTH2(scale)->Draw("colz");
202  new TCanvas;
203  extrap.NumuSurvComponent().ToTH2(scale)->Draw("colz");
204  new TCanvas;
205  extrap.AntiNumuSurvComponent().ToTH2(scale)->Draw("colz");
206  new TCanvas;
207  extrap.NueAppComponent().ToTH2(scale)->Draw("colz");
208  new TCanvas;
209  extrap.AntiNueAppComponent().ToTH2(scale)->Draw("colz");
210  new TCanvas;
211  extrap.NumuAppComponent().ToTH2(scale)->Draw("colz");
212  new TCanvas;
213  extrap.AntiNumuAppComponent().ToTH2(scale)->Draw("colz");
214 
215 /* These plots are ratios of extrapolated components to FD MC
216  The extrapolated components are collapsed to the extrapolated variable,
217  i.e., calorimetric energy
218  The FD MC are spectra that select a given component by truth
219  The same cuts are applied to both spectra, so the ratio should be 1 */
220 /* new TCanvas;
221  TH1* hNCSurv = extrap.OscNCComponent().Unoscillated().ToTH1(scale);
222  Spectrum sNCSurv = base.PredictComponent(&noosc, Flavors::kAll, Current::kNC, Sign::kBoth);
223  hNCSurv->Divide(sNCSurv.ToTH1(scale));
224  hNCSurv->SetTitle("Ratio of Extrapolation to FD MC: NC Surv");
225  hNCSurv->Draw("hist");
226  new TCanvas;
227  TH1* hNueSurv = extrap.NueSurvComponent().Unoscillated().ToTH1(scale);
228  Spectrum sNueSurv = base.PredictComponent(&noosc, Flavors::kNuEToNuE, Current::kCC, Sign::kNu);
229  hNueSurv->Divide(sNueSurv.ToTH1(scale));
230  hNueSurv->SetTitle("Ratio of Extrapolation to FD MC: #nu_{e} Surv");
231  hNueSurv->Draw("hist");
232  new TCanvas;
233  TH1* hAntiNueSurv = extrap.AntiNueSurvComponent().Unoscillated().ToTH1(scale);
234  Spectrum sAntiNueSurv = base.PredictComponent(&noosc, Flavors::kNuEToNuE, Current::kCC, Sign::kAntiNu);
235  hAntiNueSurv->Divide(sAntiNueSurv.ToTH1(scale));
236  hAntiNueSurv->SetTitle("Ratio of Extrapolation to FD MC: #bar{#nu}_{e} Surv");
237  hAntiNueSurv->Draw("hist");
238  new TCanvas;
239  TH1* hNumuSurv = extrap.NumuSurvComponent().Unoscillated().ToTH1(scale);
240  Spectrum sNumuSurv = base.PredictComponent(&noosc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kNu);
241  hNumuSurv->Divide(sNumuSurv.ToTH1(scale));
242  hNumuSurv->SetTitle("Ratio of Extrapolation to FD MC: #nu_{#mu} Surv");
243  hNumuSurv->Draw("hist");
244  new TCanvas;
245  TH1* hAntiNumuSurv = extrap.AntiNumuSurvComponent().Unoscillated().ToTH1(scale);
246  Spectrum sAntiNumuSurv = base.PredictComponent(&noosc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kAntiNu);
247  hAntiNumuSurv->Divide(sAntiNumuSurv.ToTH1(scale));
248  hAntiNumuSurv->SetTitle("Ratio of Extrapolation to FD MC: #bar{#nu}_{#mu} Surv");
249  hAntiNumuSurv->Draw("hist");
250  new TCanvas;
251  TH1* hNueApp = extrap.NueAppComponent().Unoscillated().ToTH1(scale);
252  hNueApp->Divide(sNueApp.ToTH1(scale));
253  hNueApp->SetTitle("Ratio of Extrapolation to FD MC: #nu_{e} App");
254  hNueApp->Draw("hist");
255  new TCanvas;
256  TH1* hAntiNueApp = extrap.AntiNueAppComponent().Unoscillated().ToTH1(scale);
257  hAntiNueApp->Divide(sAntiNueApp.ToTH1(scale));
258  hAntiNueApp->SetTitle("Ratio of Extrapolation to FD MC: #bar{#nu}_{e} App");
259  hAntiNueApp->Draw("hist");
260  new TCanvas;
261  TH1* hNumuApp = extrap.NumuAppComponent().Unoscillated().ToTH1(scale);
262  hNumuApp->Divide(sNumuApp.ToTH1(scale));
263  hNumuApp->SetTitle("Ratio of Extrapolation to FD MC: #nu_{#mu} App");
264  hNumuApp->Draw("hist");
265  new TCanvas;
266  TH1* hAntiNumuApp = extrap.AntiNumuAppComponent().Unoscillated().ToTH1(scale);
267  hAntiNumuApp->Divide(sAntiNumuApp.ToTH1(scale));
268  hAntiNumuApp->SetTitle("Ratio of Extrapolation to FD MC: #bar{#nu}_{#mu} App");
269  hAntiNumuApp->Draw("hist");
270 */
271 
272 /* This plot takes the extrapolation and base MC,
273  applies null oscillations, and forms a ratio.
274  If the ND Data and MC inputs to the extrapolation are identical,
275  or if the extrapolation is the TrivialExtrapNC,
276  the ratio should be 1. */
277  new TCanvas;
278  Spectrum predNoOsc = pred.PredictUnoscillated();
279  Spectrum baseNoOsc = base.PredictUnoscillated();
280  TH1* hratioNull = predNoOsc.ToTH1(scale);
281  hratioNull->Divide(baseNoOsc.ToTH1(scale));
282  hratioNull->SetTitle("Ratio of Extrapolation to FD MC (No Oscillations)");
283  hratioNull->Draw("hist");
284 
285 /* This plot takes the sterile prediction and standard 3 flavor prediction,
286  applies the same oscillations, and forms a ratio.
287  Since both predictions use the same extrapolation, the ratio should be 1. */
288  new TCanvas;
289  Spectrum predOsc = pred.Predict(calcst);
290  Spectrum pred3FlavOsc = base.Predict(calc);
291  TH1* hratioOsc = predOsc.ToTH1(scale);
292  hratioOsc->Divide(pred3FlavOsc.ToTH1(scale));
293  hratioOsc->SetTitle("Ratio of Sterile to 3 Flavor Prediction (Same Oscillations)");
294  hratioOsc->Draw("hist");
295 
296  rootOut->Close();
297 }
caf::Proxy< unsigned int > nshwlid
Definition: SRProxy.h:1999
Near Detector underground.
Definition: SREnums.h:10
Pass neutrinos through unchanged.
Definition: IOscCalc.h:40
const Var kLEM
PID
Definition: Vars.cxx:26
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2018
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetNFlavors(int nflavors)
const std::string fnamefar
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:148
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2085
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
T sqrt(T number)
Definition: d0nt_math.hpp:156
void SetDelta(int i, int j, double delta)
const Cut kIsAntiNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return sr->mc.nu[0].pdg< 0;})
Is this truly an antineutrino?
Definition: TruthCuts.h:53
const Cut kIsNumuApp(CCFlavSel(14, 12))
Select CC .
Adapt the PMNS_Sterile calculator to standard interface.
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
#define M_PI
Definition: SbMath.h:34
caf::Proxy< unsigned int > ncalhit
Definition: SRProxy.h:1268
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
const std::string fnamenear
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2077
const std::string fnameswap
const Cut kNumuQEND
Definition: NumuCuts.h:66
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2002
Double_t scale
Definition: plot.C:25
const Var kRemID
PID
Definition: Vars.cxx:81
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
Spectrum Predict(osc::IOscCalc *calc) const override
caf::Proxy< unsigned int > nmiphit
Definition: SRProxy.h:1272
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
Definition: Loaders.cxx:129
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
const Var kCCE
Definition: NumuVars.h:21
const Cut kIsSig(CCFlavSel(12, 14))
Select CC .
caf::StandardRecord * sr
osc::NoOscillations noosc
Definition: NueSystFuncs.h:44
void sterile_demo()
Definition: sterile_demo.C:29
void SetAngle(int i, int j, double th)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2017
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2101
void SetDm(int i, double dm)
virtual void SetL(double L) override
const std::string fnametau
virtual void SetRho(double rho) override
std::vector< Loaders * > loaders
Definition: syst_header.h:386
virtual Spectrum PredictUnoscillated() const override
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2032
A prediction object compatible with sterile oscillations.
Prediction that just uses FD MC, with no extrapolation.
const HistAxis axisNumu
Take the output of an extrapolation and oscillate it as required.
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2105
Float_t e
Definition: plot.C:35
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
Definition: Loaders.cxx:25
static ModularExtrapSterile NCDisappearance(Loaders &loaders, const IDecomp &NCSurvDecomp, const IDecomp &NumuOscDecomp, const HistAxis &axis, const HistAxis &axisNumuND, const Cut &fdcut, const Cut &NCNDcut, const Cut &NumuNDcut, const SystShifts &shiftMC=kNoShift, const Var &weight=kUnweighted)
Creates a NC disappearance extrapolation.
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
Just return the ND truth spectra as the decomposition.
Definition: CheatDecomp.h:10
T asin(T number)
Definition: d0nt_math.hpp:60