mre_example_test.C
Go to the documentation of this file.
1 //Short Example Script Utilizing the MRECorrection Class
2 
4 #include "CAFAna/Core/Spectrum.h"
5 #include "CAFAna/Core/ISyst.h"
10 #include "3FlavorAna/Cuts/NueCutsSecondAna.h"
17 #include "CAFAna/Vars/XsecTunes.h"
22 
23 using namespace ana;
24 using namespace nueccinc;
25 
26 //Define MRE Specific Vars & Cuts
29  if(parent.eff< 0.9||parent.pur< 0.9) return false;
30  if(parent.remid< 0.75) return false;
31  if(parent.numuE< 0|| parent.nhit< 20|| parent.contplanes< 5) return false;
32  if(parent.firstplane< 2 || parent.lastplane > 211||
33  parent.muonstop.Z() > 1275|| parent.muonfwdcell< 5 ||
34  parent.muonbkcell< 9 || parent.hadEinmucat > 0.03) return false;
35  return true;
36  });
37 
38 const ana::Var kTrueNu([](const caf::SRProxy *sr){
39  if(sr->mc.nu.size() == 0) return -5.0f;
40  if(!sr->mc.nu[0].iscc) return -5.0f;
41 
42  else return (float)sr->mc.nu[0].E;
43  });
44 
45 const ana::Var kTrueParentMuonE([](const caf::SRProxy* sr){
46  if(sr->mc.nu.size() == 0) return -5.0f;
47  if(!(sr->mc.nu[0].iscc)) return -5.0f;
48  if(!(sr->mc.nu[0].pdg == 14)) return -5.0f;
49 
50  float ee = -5;
51 
52  int nprims = sr->mc.nu[0].prim.size();
53  for(int iprim = 0; iprim < nprims; iprim++){
54  if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
55  float E= (float)sr->mc.nu[0].prim[iprim].p.E;
56  ee = E;
57  break;
58  }
59  }
60  return (float)ee;
61  });
62 
63 const ana::Var kTrueParentMuonCos([](const caf::SRProxy* sr){
64  if(sr->mc.nu.size() == 0) return -5.0f;
65  if(!sr->mc.nu[0].iscc)return -5.0f;
66  if(!(sr->mc.nu[0].pdg == 14)) return -5.0f;
67  int nprims = sr->mc.nu[0].prim.size();
68  for(int iprim = 0; iprim < nprims; iprim++){
69  if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
70  TVector3 edir =
71  sr->mc.nu[0].prim[iprim].p.Vect();
72  TVector3 beamdir =
74  return
75  (float)edir.Unit().Dot(beamdir.Unit());
76  }
77  }
78  return -5.0f;
79  });
80 
81 
82 
87 
88 const int kvMREAnalysisVars = 2;
89 
91  {
92  {"true_e", {"True Neutrino Energy (GeV)",enubins, kTrueNu}},
93  {"true_muon_double", {"; True E_{#mu} (GeV); cos #{theta}",
94  double_diff_bins, kTrueMuonEVsCos}},
95 
96  };
97 
98 void mre_example_test(std::string outFileName = "fMRECorrection_test.root")
99 {
100  std::vector<ana::SpectrumLoader*> lStandard;
101  std::vector<std::string> dataset_names;
102  std::vector<std::string> dirNames;
103 
105 
106  dataset_names.push_back(ana::nueccinc::MRE_NominalMC);
107  dirNames.push_back("Nominal");
108 
109  //Set Up Loader(s)
110  for(auto name : dataset_names){
111  lStandard.push_back(new ana::SpectrumLoader(name));
112  lStandard.back()->SetSpillCut(ana::kStandardSpillCuts);
113  }
114 
115  //Holder for MRE Correction Class
116  std::map<std::string,NueCCIncMRECorrection*> sMREAnalysis;
117 
118  //Analysis Weights
120 
121 
122  //Initialize MRE Correction Spectra
123  for(uint iLoad = 0; iLoad < lStandard.size(); iLoad++){
124 
125  std::string outName = dirNames[iLoad];
126 
127  sMREAnalysis.insert(std::make_pair
128  (
129  outName,
131  (
132  *lStandard[iLoad], //Loader
133  lNDMC, //Loader we want to compare too
134  analysis_vars[0].axis, //1D hist axis
135  kRecoElecEVsCosStandardAxis, //2D hist axis
136  mre_analysis_vars[0].axis, //True muon 1D axis
137  mre_analysis_vars[1].axis, //True muon 2D axis
138  kcDQ &&
140  kMREParentSliceCut, //"Fiducially" selected
142  kRecoElectronPhaseSpaceCut, //"Fully selected"
143  kNoShift, //loader shifts
144  kNoShift, //lNDMC shifts
145  kCVWeight, //loader weights
146  kCVWeight) //lNDMC weights
147  ));
148 
149 
150  }//loop over loaders
151 
152  //Class Creates a matrix which relates the Reconstructed Electron
153  //Kinematics to the True Muon Kinematics in MC-only.
154  //This is used to "unfold" the reconstructed space of lStandard[iLoad]
155  //to the true lepton kinematic space.
156 
157 
158  //Let's Go
159  for(uint iLoad = 0; iLoad < lStandard.size(); iLoad++){
160  lStandard[iLoad]->Go();
161  }
162  lNDMC.Go();
163 
164  //Save Output
165  TFile* fout =
166  new TFile(outFileName.c_str(),"recreate");
167 
168  for( std::pair<std::string, NueCCIncMRECorrection*> ele : sMREAnalysis)
169  {
170  ele.second->SaveTo(fout->mkdir(ele.first.c_str()));
171  }
172 
173  fout->Close();
174 }
Near Detector underground.
Definition: SREnums.h:10
const Binning enubins
Definition: NueCCIncBins.h:18
const XML_Char * name
Definition: expat.h:151
const ana::nueccinc::mHistAxisDef mre_analysis_vars[kvMREAnalysisVars]
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const mHistAxisDef analysis_vars[kvAnalysisVars]
Definition: NueCCIncVars.h:478
caf::Proxy< float > hadEinmucat
Definition: SRProxy.h:739
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Cut kcPresel
Definition: NueCCIncCuts.h:139
caf::Proxy< caf::SRParentBranch > parent
Definition: SRProxy.h:2140
caf::Proxy< int > muonfwdcell
Definition: SRProxy.h:743
Proxy for caf::SRMRCCParent.
Definition: SRProxy.h:725
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
_Var< T > Var2D(const _Var< T > &a, const Binning &binsa, const _Var< T > &b, const Binning &binsb)
Variable formed from two input variables.
Definition: Var.cxx:247
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
void abs(TH1 *hist)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
const ana::Var kTrueNu([](const caf::SRProxy *sr){if(sr->mc.nu.size()==0) return-5.0f;if(!sr->mc.nu[0].iscc) return-5.0f;else return(float) sr->mc.nu[0].E;})
const Binning eelecbins
Definition: NueCCIncBins.h:23
const std::string MRE_NominalMC
Definition: NueCCIncExtra.h:87
const int kvMREAnalysisVars
caf::Proxy< int > lastplane
Definition: SRProxy.h:740
const Cut kcDQ
Definition: NueCCIncCuts.h:114
Float_t E
Definition: plot.C:20
const ana::Cut kMREParentSliceCut([](const caf::SRProxy *sr){const caf::SRMRCCParentProxy &parent=sr->parent.mrccpar;if(parent.eff< 0.9||parent.pur< 0.9) return false;if(parent.remid< 0.75) return false;if(parent.numuE< 0||parent.nhit< 20||parent.contplanes< 5) return false;if(parent.firstplane< 2||parent.lastplane > 211|| parent.muonstop.Z() > 1275||parent.muonfwdcell< 5|| parent.muonbkcell< 9||parent.hadEinmucat > 0.03) return false;return true;})
caf::Proxy< int > contplanes
Definition: SRProxy.h:735
const ana::Var kTrueParentMuonE([](const caf::SRProxy *sr){if(sr->mc.nu.size()==0) return-5.0f;if(!(sr->mc.nu[0].iscc)) return-5.0f;if(!(sr->mc.nu[0].pdg==14)) return-5.0f;float ee=-5;int nprims=sr->mc.nu[0].prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(sr->mc.nu[0].prim[iprim].pdg)==13){float E=(float) sr->mc.nu[0].prim[iprim].p.E;ee=E;break;}}return(float) ee;})
const Binning double_diff_bins
Definition: NueCCIncBins.h:63
void mre_example_test(std::string outFileName="fMRECorrection_test.root")
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
const Binning costhetabins
Definition: NueCCIncBins.h:27
caf::Proxy< int > nhit
Definition: SRProxy.h:748
const SystShifts kNoShift
Definition: SystShifts.cxx:22
caf::Proxy< float > eff
Definition: SRProxy.h:737
const HistAxis kRecoElecEVsCosStandardAxis("Reconstructed E_{e} vs cos #theta; Elec_{e} (GeV); cos #{theta}", double_diff_bins, kRecoElecEVsCos)
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
caf::Proxy< caf::SRMRCCParent > mrccpar
Definition: SRProxy.h:766
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
caf::Proxy< float > pur
Definition: SRProxy.h:750
caf::Proxy< caf::SRVector3D > muonstop
Definition: SRProxy.h:745
const ana::Var kTrueMuonEVsCos
TVector3 NuMIBeamDirection(caf::Det_t det)
Average direction of NuMI neutrinos in a given detector This function is a copy of geo::GeometryBase:...
Definition: Utilities.cxx:329
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
caf::Proxy< float > numuE
Definition: SRProxy.h:749
caf::Proxy< int > firstplane
Definition: SRProxy.h:738
const ana::Var kTrueParentMuonCos([](const caf::SRProxy *sr){if(sr->mc.nu.size()==0) return-5.0f;if(!sr->mc.nu[0].iscc) return-5.0f;if(!(sr->mc.nu[0].pdg==14)) return-5.0f;int nprims=sr->mc.nu[0].prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(sr->mc.nu[0].prim[iprim].pdg)==13){TVector3 edir=sr->mc.nu[0].prim[iprim].p.Vect();TVector3 beamdir=NuMIBeamDirection(caf::kNEARDET);return(float) edir.Unit().Dot(beamdir.Unit());}}return-5.0f;})
caf::Proxy< int > muonbkcell
Definition: SRProxy.h:742
const Cut kRecoElectronPhaseSpaceCut([](const caf::SRProxy *sr){float costheta=kvRecoCos(sr);float elece=kvElecE(sr);if((costheta >=0.97)&&(elece >=1.4)&&(elece< 6.0)) return true;if((costheta >=0.94)&&(costheta< 0.97)&& (elece >=1.4)&&(elece< 4.1)) return true;if((costheta >=0.90)&&(costheta< 0.94)&& (elece >=1.0)&&(elece< 2.5)) return true;if((costheta >=0.85)&&(costheta< 0.90)&& (elece >=1.0)&&(elece< 2.0)) return true;return false;})
caf::Proxy< float > remid
Definition: SRProxy.h:751
enum BeamMode string
unsigned int uint
const Cut kcNueCCIncFiducial([](const caf::SRProxy *sr){assert(sr->vtx.nelastic > 0 &&"Must apply DQ cuts");if(sr->vtx.elastic[0].vtx.X()< -130.0) return false;if(sr->vtx.elastic[0].vtx.X() > 150.0) return false;if(sr->vtx.elastic[0].vtx.Y()< -140.0) return false;if(sr->vtx.elastic[0].vtx.Y() > 140.0) return false;if(sr->vtx.elastic[0].vtx.Z()< 150.0) return false;if(sr->vtx.elastic[0].vtx.Z() > 800.0) return false;return true;})