datamc_ND_numu_kinematics.C
Go to the documentation of this file.
4 
5 // Core
8 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Core/HistAxis.h"
11 #include "CAFAna/Core/Binning.h"
12 #include "CAFAna/Core/Cut.h"
13 
14 // Cuts
15 #include "CAFAna/Cuts/Cuts.h"
16 #include "CAFAna/Cuts/SpillCuts.h"
20 #include "CAFAna/Cuts/NumuCuts.h"
21 #include "CAFAna/Cuts/TruthCuts.h"
22 
23 // Vars
24 #include "CAFAna/Vars/Vars.h"
27 #include "CAFAna/Vars/NueVars.h"
28 #include "CAFAna/Vars/NumuVars.h"
30 
31 // Systs
32 #include "CAFAna/Systs/Systs.h"
33 #include "CAFAna/Analysis/Style.h"
34 #include "CAFAna/Decomp/IDecomp.h"
36 
38 
39 
40 #include "CAFAna/Core/SystShifts.h"
41 #include "CAFAna/Cuts/NumuCuts.h"
47 #include "TFile.h"
50 #include "CAFAna/Core/Spectrum.h"
51 #include "TH1.h"
53 #include "TCanvas.h"
55 #include "CAFAna/Analysis/Plots.h"
56 #include "CAFAna/Vars/FitVars.h"
57 #include "CAFAna/Cuts/Cuts.h"
58 #include "CAFAna/Cuts/NumuCuts.h"
59 
60 #include "TStopwatch.h"
61 #include "CAFAna/Core/Utilities.h"
62 
63 
64 
65 #include <fstream>
66 #include <iostream>
67 #include <cmath>
68 #include <string>
69 
70 using namespace ana;
71 
73 
74 // ROOT
75 #include "TCanvas.h"
76 #include "TFile.h"
77 #include "TH1.h"
78 #include "TLegend.h"
79 #include "TPad.h"
80 #include "TMath.h"
81 #include "TH2.h"
82 #include "TLine.h"
83 #include "TText.h"
84 #include "TString.h"
85 
86 #include "Utilities/func/MathUtil.h"
87 
88 
90 {
91 
92  std::string nameData = "defname: prod_decaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_nue_or_numu_or_nus_contain_v1_goodruns";
93  std::string nameMC = "defname: prod_decaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
94 
95  //choose your samweb datasets or concat files:
96  SpectrumLoader loaderData(nameData);
97  SpectrumLoader loaderMC(nameMC);
98 
99  loaderData.SetSpillCut(kStandardSpillCuts);
100  loaderMC.SetSpillCut(kStandardDQCuts);
101 
102  //which variables do you want to look at?
103  const int kNumVars = 7;
104 
105  const HistDef defs[kNumVars] = {
106  {"recoE", {"E_{reco} (GeV)", Binning::Simple(100, 0, 5), kCCE}},
107  {"trueQ2", {"true Q^{2} (GeV^{2})", Binning::Simple(500, 0, 5), kTrueQ2}},
108  {"recoQ2", {"reco Q^{2} (GeV^{2})", Binning::Simple(500, 0, 5), kRecoQ2}},
109  {"trueW2", {"true W^{2} (GeV^{2})", Binning::Simple(100, 0, 2), kTrueW}},
110  {"trueE", {"true E (GeV)", Binning::Simple(100, 0, 5), kTrueE}},
111  {"PtP", {"p_{t}/p", Binning::Simple(500,0,1), kPtP}},
112  {"CosNumi", {"CosNumi", Binning::Simple(500,0,1), kCosNumi}}
113  };
114 
115  const Cut kIsRock(
116  [](const caf::SRProxy* sr)
117  {
118  if (sr->mc.nnu == 0)
119  return false;
120  if (fabs(sr->mc.nu[0].vtx.X())>180 ||
121  fabs(sr->mc.nu[0].vtx.Y())>180 ||
122  sr->mc.nu[0].vtx.Z()<0 ||
123  sr->mc.nu[0].vtx.Z()>1250)
124  return true;
125  return false;
126  });
127 
128  const int kNumSels = 1;
129  const Cut sels[kNumSels] = {
130  kNumuNDCvn
131  };
132  const std::string selNames[kNumSels] = {
133  "numuND"
134  };
135 
137  IPrediction* preds[kNumSels][kNumVars];
138 
139  //loop over the selectors and variables to create a full selection
140  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
141  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
142  const HistAxis& axis = defs[varIdx].axis;
143  spects[selIdx][varIdx] = new Spectrum(loaderData, axis, sels[selIdx], kNoShift);
144  preds[selIdx][varIdx] = new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
145  axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0],
146  sels[selIdx], kNoShift, kXSecCVWgt2017*kPPFXFluxCVWgt);
147  }
148  }
149 
150  loaderMC.Go();
151  loaderData.Go();
152 
153  TString fname = "datamc_ND_numu_kinematics.root";
154  TFile* fout = new TFile(fname, "RECREATE");
155  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
156  TDirectory* d = fout->mkdir(selNames[selIdx].c_str());
157  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
158  const char* name = defs[varIdx].name.c_str();
159  spects[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("spect_%s", name)));
160  preds[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("pred_%s", name)));
161  }
162  }
163 }
const XML_Char * name
Definition: expat.h:151
const int kNumVars
Definition: vars.h:14
Oscillation analysis framework, runs over CAF files outside of ART.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
VectorProxy< SRNeutrinoProxy > nu
Definition: SRProxy.h:1890
virtual void SaveTo(TDirectory *dir) const
std::vector< T > GetVars() const
Definition: HistAxis.h:79
void datamc_ND_numu_kinematics()
HistAxis axis
Definition: NuePlotLists.h:13
Proxy for StandardRecord.
Definition: SRProxy.h:2237
void SetSpillCut(const SpillCut &cut)
const Var kTrueQ2
Definition: TruthVars.h:27
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
const Var kPtP
Transverse momentum fraction in slice.
Definition: NueVars.cxx:90
Proxy< short int > nnu
Definition: SRProxy.h:1899
const Var kTrueE([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].E);})
Definition: Vars.cxx:87
if(dump)
const Cut sels[kNumSels]
Definition: vars.h:44
const Cut kNumuNDCvn
Definition: NumuCuts.h:61
const SpillCut kStandardDQCuts([](const caf::SRSpillProxy *spill){if(spill->dcmedgematchfrac==0 &&spill->fracdcm3hits==0 &&spill->nmissingdcmslg==0) return bool(spill->isgoodspill); if(spill->det==caf::kNEARDET && (spill->fracdcm3hits > 0.45|| spill->nmissingdcms > 0)) return false; if(spill->eventincomplete) return false; if(spill->det==caf::kFARDET && spill->nmissingdcmslg > 0) return false; if(spill->det==caf::kFARDET && !spill->ismc && spill->dcmedgematchfrac<=0.2) return false;return true;})
Cut out events with a noisy detector or with parts missing.
Definition: SpillCuts.h:12
const Var kCCE
Definition: Vars.h:90
Float_t d
Definition: plot.C:236
virtual void Go() override
Load all the registered spectra.
const HistDef defs[kNumVars]
Definition: vars.h:15
const int kNumSels
Definition: vars.h:43
std::string name
Definition: NuePlotLists.h:12
const SystShifts kNoShift
Definition: SystShifts.h:112
const ana::Var kRecoQ2([](const caf::SRProxy *sr){const double M_mu_sqrd=util::sqr(0.1056);double E_mu=kNumuTrackE(sr);double p_mu=sqrt(util::sqr(E_mu)-M_mu_sqrd);return 2 *kCCE(sr)*(E_mu-p_mu *kCosNumi(sr))-M_mu_sqrd;})
Reconstructed four-momentum transfer invariant (Q^2)
Definition: NumuVars.h:128
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
SRTruthBranchProxy mc
Definition: SRProxy.h:2253
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:28
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void spects(osc::IOscCalculatorAdjustable *calc, IPrediction *predFHC, IPrediction *predRHC)
Definition: sensitivity.C:100
Prediction that just uses FD MC, with no extrapolation.
const std::string selNames[kNumSels]
Definition: vars.h:46
const Var kXSecCVWgt2017
Definition: XsecTunes.h:38
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
const Var kTrueW
Definition: TruthVars.h:22
std::vector< std::string > GetLabels() const
Definition: HistAxis.h:77
const Var kCosNumi([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999){if(sr->hdr.det==1){return sr->trk.kalman.tracks[0].dir.Dot(beamDirND);}if(sr->hdr.det==2){return sr->trk.kalman.tracks[0].dir.Dot(beamDirFD);}}return-5.f;})
Definition: NumuVars.h:27
static NullLoader kNullLoader
Dummy loader that doesn&#39;t load any files.
std::vector< Binning > GetBinnings() const
Definition: HistAxis.h:78
static constexpr Double_t sr
Definition: Munits.h:164
void SaveTo(TDirectory *dir) const
Definition: Spectrum.cxx:1029