test_xsecanalysis.C
Go to the documentation of this file.
4 
7 //#include "CAFAna/Core/LoadFromFile.h"
8 
11 
12 #include "CAFAna/Vars/Vars.h"
13 
14 #include "TVector3.h"
15 #include "TCanvas.h"
16 #include "TFile.h"
17 #include "TH1.h"
18 #include "TH2.h"
20 
21 #include <cmath>
22 #include <iostream>
23 
24 //#include "Utilities/func/MathUtil.h"
25 
26 using namespace ana;
27 
28 // Demonstrative In-n-Out test macro that puts together the MC prediction
29 // for the flux-averaged muon energy distribution for all numu CC events
30 
31 // Entirely based on properly defining the inputs and setting up a
32 // CrossSectionAnalysis. That class does all the algebra that you need, and
33 // once you have one filled and loaded, it's easy to make interesting plots
34 
35 // Just a pedagogical shell, a lot of work will go into thinking about the
36 // specifics of each analysis, and the final CrossSectionAnalysis you put
37 // together will be more involved.
38 
39 
40 
41 // Need to consider all simulated events, not just those sliced + reco'd.
42 // The true muon energy variable has to come from the neutrino
44  if (sr->pdg != 14) return float(-1); // Require event be numu
45  if (!sr->iscc) return float(-1); // Require event be CC
46  if (sr->prim.size() == 0) return float(-1);
47  return float(sr->prim[0].p.T());
48  });
50 
51 // Signal definition - all numu CC
52 const NuTruthCut kIsNuMuCC_NT([](const caf::SRNeutrinoProxy* sr){
53  if (sr->pdg != 14) return false;
54  return bool(sr->iscc);
55  });
57 
58 // Analysis will have fiducial cut - define one here
59 const Cut kIsFiducial([](const caf::SRProxy* sr){
60  if (!sr->vtx.elastic.IsValid) return false;
61  TVector3 vtx(sr->vtx.elastic.vtx);
62  if (vtx[0] > 100) return false;
63  if (vtx[1] > 100) return false;
64  if (vtx[2] < 100 || vtx[2] > 1000) return false;
65  return true;
66  });
67 
68 
69 void test_xsecanalysis(bool makePlots=false)
70 {
71  if (!makePlots){
72  // Set up spectrum loaders - for an ND analysis, you need to have ND
73  // data and MC
74 
75  // Faster to copy a file to working directory for testing
76  // cp /pnfs/nova/production/caf/R17-11-14-prod4reco.d/genie/nd/000105/10565/neardet_genie_nonswap_genierw_fhc_v08_2157_r00010565_s07_c003_R17-11-14-prod4reco.d_v1_20170322_204739_sim.caf.root .
77  std::string test_file = "/pnfs/nova/production/caf/R17-11-14-prod4reco.d/genie/nd/000105/10565/neardet_genie_nonswap_genierw_fhc_v08_2157_r00010565_s07_c003_R17-11-14-prod4reco.d_v1_20170322_204739_sim.caf.root";
78  SpectrumLoader lMC(test_file);
79  SpectrumLoader lDa(test_file); // In-n-Out test, use MC fake-data until you're ready for data
82 
83  // For estimating data - you need an axis in your reco'd variable
84  HistAxis axis("Muon Energy [GeV]",Binning::Simple(20,0,5),kMuE);
85  // For evaluating the GENIE xsec, you need an axis in the true variable
86  NuTruthHistAxis axisNuTruth("Muon Energy [GeV]",Binning::Simple(20,0,5),
87  kTrueMuE_NT);
88  Binning EBins = Binning::Simple(40,0,10); // For estimating flux
89  Cut sel = kNumuCutND2017 && kIsFiducial; // Analysis cuts
90 
91  // You need some way to estimate your background. For a real analysis,
92  // the specifics will be different, so put thought into the procedure.
93  // For an in-n-out test, this class will return the MC bkgd estimate.
94  TrivialBkgdEstimator *bkgdest =
95  new TrivialBkgdEstimator(lMC, axis, sel, {!kIsNuMuCC});
96 
97  // Analysis will have a fiducial volume, need to define for flux estimate
98  TVector3 fidMin(-100,-100,100);
99  TVector3 fidMax(100,100,1000);
100 
101  // Unfolding won't matter for an in-n-out test, but in a real analysis,
102  // studying this will be important
103  double unfoldReg = 2;
104  UnfoldMethod_t unfoldType = kIterative;
105 
106  // And we're ready to make the xsec class
107  // Pretty simple, no SystShifts / MC weighting. Real analysis will likely
108  // include these
109  CrossSectionAnalysis xsec(lMC, lDa, axis, axisNuTruth, EBins,
110  sel, kIsNuMuCC_NT, bkgdest, fidMin, fidMax,
111  unfoldReg, unfoldType);
112 
113  // Save to a file to analyize
114  TFile *out = new TFile("out_xsectest.root","recreate");
115  xsec.SaveTo(out, "xsec");
116  out->Close();
117 
118  }
119  else {
120  TFile *f = TFile::Open("out_xsectest.root");
121  f->cd("xsec");
122 
123  // Load from your file
125  CrossSectionAnalysis::LoadFrom(f, "xsec").release();
126 
127  // Now we can start plotting
128  new TCanvas();
129  xsec->Result()->Draw("hist");
130 
131  new TCanvas();
132  xsec->PlotEfficiency()->Draw("hist");
133 
134  new TCanvas();
135  xsec->PlotFluxEstimate()->Draw("hist");
136 
137  new TCanvas();
138  xsec->PlotUnfoldedSignal()->Draw();
139 
140  new TCanvas();
141  xsec->PlotRecoToTrueMatrix()->Draw("colz");
142 
143  }
144 
145 }
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
void makePlots()
Definition: makePlots.C:17
void SaveTo(TDirectory *dir, const std::string &name) const
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
void SetSpillCut(const SpillCut &cut)
const NuTruthCut kIsNuMuCC_NT([](const caf::SRNeutrinoProxy *sr){if(sr->pdg!=14) return false;return bool(sr->iscc);})
const Var kTrueMuE
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
const NuTruthVar kTrueMuE_NT([](const caf::SRNeutrinoProxy *sr){if(sr->pdg!=14) return float(-1);if(!sr->iscc) return float(-1);if(sr->prim.size()==0) return float(-1);return float(sr->prim[0].p.T());})
if(dump)
static std::unique_ptr< CrossSectionAnalysis > LoadFrom(TDirectory *dir, const std::string &name, UnfoldMethod_t unfoldmethod=kIterative)
caf::Proxy< bool > iscc
Definition: SRProxy.h:538
caf::StandardRecord * sr
const Cut kNumuCutND2017
Definition: NumuCuts2017.h:41
void test_xsecanalysis(bool makePlots=false)
const Cut kIsFiducial([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return false;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return false;bool isfid=(sr->trk.kalman.tracks[ibesttrk].start.X()< vtxmax.X()&&sr->trk.kalman.tracks[ibesttrk].start.X() > vtxmin.X()&&sr->trk.kalman.tracks[ibesttrk].start.Y() > vtxmin.Y()&&sr->trk.kalman.tracks[ibesttrk].start.Y()< vtxmax.Y()&&sr->trk.kalman.tracks[ibesttrk].start.Z() > vtxmin.Z()&&sr->trk.kalman.tracks[ibesttrk].start.Z()< vtxmax.Z());return isfid;})
Definition: NumuCCIncCuts.h:27
Double_t xsec[nknots]
Definition: testXsec.C:47
UnfoldMethod_t
Enumerator for unfolding methods. TODO: Allow swapping to RooUnfold, TUnfold, etc?
const Cut kIsNuMuCC
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:7
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2073
Generic organizational class for a cross section analysis.
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
caf::Proxy< short int > pdg
Definition: SRProxy.h:552
Template for Cut and SpillCut.
Definition: Cut.h:15
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
Template for Var and SpillVar.
const Var kMuE
Definition: NumuVars.h:22
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
Just return the MC prediction for the background.
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:7
caf::Proxy< std::vector< caf::SRTrueParticle > > prim
Definition: SRProxy.h:555
enum BeamMode string