hyperon_nom_macro.C
Go to the documentation of this file.
1 #pragma once
2 
3 #include "CAFAna/Core/Var.h"
4 #include "CAFAna/Core/MultiVar.h"
5 #include "CAFAna/Core/Binning.h"
6 #include "CAFAna/Core/Spectrum.h"
9 
10 #include "CAFAna/Cuts/Cuts.h"
16 #include "CAFAna/Cuts/SpillCuts.h"
17 
20 #include "CAFAna/Vars/Vars.h"
23 
27 
29 
30 #include "Utilities/rootlogon.C"
31 
32 #include "TCanvas.h"
33 #include "TH2.h"
34 #include "TProfile.h"
35 #include "TSystem.h"
36 
37 using namespace ana;
38 
39 void hyperon_nom_macro(bool IsFarDet = true)
40 {
41 
42  std::string sFD_ND = ( IsFarDet == true ? "FD" : "ND" );
43  std::string sOutFile = "nom_"+sFD_ND+".root";
44 
45  // --- Loaders ---
46  std::string nom = "";
47  if (IsFarDet) {
48  nom = "prod_caf_R17-11-14-prod4reco.neutron-respin.c_fd_genie_nonswap_rhc_nova_v08_full_v1"; //
49  }
50  else {
51  nom = "prod_caf_R17-11-14-prod4reco.neutron-respin.c_nd_genie_nonswap_rhc_nova_v08_full_v1";
52  }
53 
54  SpectrumLoader ldr(nom);
56 
57  // --- Weights ---
59 
60  // --- Interaction types ---
61  const std::vector<Cut> GENIECuts = { kNoCut, kIsDytmanMEC, kIsDIS, kIsRes, kIsCoh, kIsQE };
62  const std::vector<std::string> GENIEStr = { "All", "MEC", "DIS", "RES", "COH", "QE"};
63  const unsigned int nGENIE = GENIECuts.size();
64 
65  // --- Selection Cuts ---
66  std::vector<std::string> cutNames; std::vector<Cut> kDetCuts;
67  if (IsFarDet) {
68  cutNames.emplace_back("kNoCut"); kDetCuts.emplace_back(kNoCut);
69  cutNames.emplace_back("qual"); kDetCuts.emplace_back(kNumuQuality);
70  cutNames.emplace_back("qual+cont"); kDetCuts.emplace_back(kNumuQuality && kNumuContainFD2017);
71  cutNames.emplace_back("qual+cont+PID"); kDetCuts.emplace_back(kNumuQuality && kNumuContainFD2017 && kNumuPID2018);
72  cutNames.emplace_back("full"); kDetCuts.emplace_back(kNumuCutFD2018);
73  }
74  else {
75  cutNames.emplace_back("kNoCut"); kDetCuts.emplace_back(kNoCut);
76  cutNames.emplace_back("qual"); kDetCuts.emplace_back(kNumuQuality);
77  cutNames.emplace_back("qual+cont"); kDetCuts.emplace_back(kNumuQuality && kNumuContainND2017);
78  cutNames.emplace_back("full"); kDetCuts.emplace_back(kNumuQuality && kNumuContainND2017 && kNumuPID2018);
79  }
80  unsigned int nDet = kDetCuts.size();
81 
82  // --- Quantiles ---
83  auto cvmfs_dir = std::getenv("NUMUDATA_DIR");
84  if (!cvmfs_dir) { std::cerr << "Couldn't find UPS dir for numu prediction!" << std::endl; return; }
85  std::string fdspecfile = std::string(cvmfs_dir) + "/ana2018/Quantiles/quantiles__rhc_full__numu2018.root";
86 
87  TFile* inFile = TFile::Open(fdspecfile.c_str());
88  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny( "FDSpec2D" );
89  const int NHadEFracQuantiles = 4;
90  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
91  HadEFracQuantCuts.push_back( kNoCut );
92  std::vector<std::string> quantNames = { "Quant1", "Quant2", "Quant3", "Quant4", "AllQuant" };
93  const unsigned int nQ = quantNames.size();
94 
95  // --- Additional Vars ---
96  // reco - true
97  const Var kRmT([](const caf::SRProxy* sr){
98  if (sr->mc.nnu == 0) return -5.;
99  return kCCE(sr) - kTrueE(sr);
100  });
101 
102  // (reco - true) / true
103  const Var kRmToT([](const caf::SRProxy* sr){
104  if (sr->mc.nnu == 0) return -5.;
105  if(kTrueE(sr) == 0) return -5.;
106  return (kCCE(sr) - kTrueE(sr))/kTrueE(sr);
107  });
108 
109  // number of prongs
110  const Var kNprongs([](const caf::SRProxy* sr){
111  double npng = -5;
112  if( !sr->vtx.elastic.IsValid ) return npng;
113  npng = sr->vtx.elastic.fuzzyk.npng;
114  return npng;
115  });
116 
117  // W
118  const Var kW([](const caf::SRProxy* sr){
119  double W = -5;
120  if(sr->mc.nnu == 0) return W;
121  if(sr->mc.nu[0].W2<0) return W;
122  W = sqrt(sr->mc.nu[0].W2);
123  return W;
124  });
125 
126  // vertex X
127  const Var kVx([](const caf::SRProxy* sr){
128  double pos = -5;
129  if(sr->mc.nnu == 0) return pos;
130  pos = sr->mc.nu[0].vtx.x;
131  return pos;
132  });
133 
134  // vertex Y
135  const Var kVy([](const caf::SRProxy* sr){
136  double pos = -5;
137  if(sr->mc.nnu == 0) return pos;
138  pos = sr->mc.nu[0].vtx.y;
139  return pos;
140  });
141 
142  // vertex Z
143  const Var kVz([](const caf::SRProxy* sr){
144  double pos = -5;
145  if(sr->mc.nnu == 0) return pos;
146  pos = sr->mc.nu[0].vtx.z;
147  return pos;
148  });
149 
150  // --- Binning ---
151  const Binning RmT_bins = Binning::Simple(80, -2, 2);
152  const Binning RmToT_bins = Binning::Simple(200, -1, 1);
153  const Binning mu_bins = Binning::Simple(100, 0, 5);
154  const Binning had_bins = Binning::Simple(40, 0, 2);
155  const Binning png_bins = Binning::Simple(10, 0, 9);
156  const Binning W_bins = Binning::Simple(60, 0, 3);
157  const Binning id_bins = Binning::Simple(120, -0.1, 1.1);
158  const Binning xybins = Binning::Simple(160, -800, 800);
159  const Binning zbins = Binning::Simple(610, 0, 6100);
160  const Binning TLbins = Binning::Simple(600,0.0,60.0);
161  const Binning Scatbins = Binning::Simple(200,-0.5,0.5);
162  const Binning dEdxbins = Binning::Simple(200,-3.0,2.0);
163  const Binning MFbins = Binning::Simple(120,-0.1,1.1);
164 
165  // --- Spectra ---
166  Spectrum *RmT_nom [nDet][nQ][nGENIE];
167  Spectrum *RmToT_nom [nDet][nQ][nGENIE];
168  Spectrum *trueHadE_nom [nDet][nQ][nGENIE];
169  Spectrum *recoHadE_nom [nDet][nQ][nGENIE];
170  Spectrum *png_nom [nDet][nQ][nGENIE];
171  Spectrum *trueNuE_nom [nDet][nQ][nGENIE];
172  Spectrum *recoNuE_nom [nDet][nQ][nGENIE];
173  Spectrum *trueMuonE_nom [nDet][nQ][nGENIE];
174  Spectrum *recoMuonE_nom [nDet][nQ][nGENIE];
175  Spectrum *remid_nom [nDet][nQ][nGENIE];
176  Spectrum *cvnprod3_nom [nDet][nQ][nGENIE];
177  Spectrum *cvn2017_nom [nDet][nQ][nGENIE];
178  Spectrum *Vx_nom [nDet][nQ][nGENIE];
179  Spectrum *Vy_nom [nDet][nQ][nGENIE];
180  Spectrum *Vz_nom [nDet][nQ][nGENIE];
181  Spectrum *trklen_nom [nDet][nQ][nGENIE];
182  Spectrum *scatLL_nom [nDet][nQ][nGENIE];
183  Spectrum *dedxLL_nom [nDet][nQ][nGENIE];
184  Spectrum *MF_nom [nDet][nQ][nGENIE];
185 
186  for (unsigned int det=0; det<nDet; ++det) {
187  for (unsigned int quant=0; quant<nQ; ++quant) {
188  for (unsigned int gen=0; gen<nGENIE; ++gen){
189  Cut nom_cut = GENIECuts[gen] && kDetCuts[det] && HadEFracQuantCuts[quant];
190  std::string nomStr = " for Cut "+cutNames[det]+", "+quantNames[quant] +", "+GENIEStr[gen];
191 
192  RmT_nom [det][quant][gen] = new Spectrum("Reco - True (GeV)"+nomStr, RmT_bins, ldr, kRmT, nom_cut, kNoShift, weight);
193  RmToT_nom [det][quant][gen] = new Spectrum("(Reco - True)/True"+nomStr, RmToT_bins, ldr, kRmToT, nom_cut, kNoShift, weight);
194  trueHadE_nom [det][quant][gen] = new Spectrum("true nu - reco muon"+nomStr, had_bins, ldr, kTrueE-kMuE, nom_cut);
195  recoHadE_nom [det][quant][gen] = new Spectrum("Hadronic Energy"+nomStr, had_bins, ldr, kHadE, nom_cut);
196  png_nom [det][quant][gen] = new Spectrum("Number of Prongs"+nomStr, png_bins, ldr, kNprongs, nom_cut);
197  trueNuE_nom [det][quant][gen] = new Spectrum("kTrueE"+nomStr, mu_bins, ldr, kTrueE, nom_cut);
198  recoNuE_nom [det][quant][gen] = new Spectrum("kCCE"+nomStr, mu_bins, ldr, kCCE, nom_cut);
199  trueMuonE_nom [det][quant][gen] = new Spectrum("kTrueMuonE"+nomStr, mu_bins, ldr, kTrueMuonE, nom_cut);
200  recoMuonE_nom [det][quant][gen] = new Spectrum("kMuE"+nomStr, mu_bins, ldr, kMuE, nom_cut);
201  remid_nom [det][quant][gen] = new Spectrum("remid"+nomStr, id_bins, ldr, SIMPLEVAR(sel.remid.pid), nom_cut);
202  cvnprod3_nom [det][quant][gen] = new Spectrum("cvnprod3"+nomStr, id_bins, ldr, SIMPLEVAR(sel.cvnProd3Train.numuid), nom_cut);
203  cvn2017_nom [det][quant][gen] = new Spectrum("cvn2017"+nomStr, id_bins, ldr, SIMPLEVAR(sel.cvn2017.numuid), nom_cut);
204  Vx_nom [det][quant][gen] = new Spectrum("true vertex x"+nomStr, xybins, ldr, kVx, nom_cut);
205  Vy_nom [det][quant][gen] = new Spectrum("true vertex y"+nomStr, xybins, ldr, kVy, nom_cut);
206  Vz_nom [det][quant][gen] = new Spectrum("true vertex z"+nomStr, zbins, ldr, kVz, nom_cut);
207  trklen_nom [det][quant][gen] = new Spectrum("trklen"+nomStr, TLbins, ldr, kTrkLength, nom_cut);
208  scatLL_nom [det][quant][gen] = new Spectrum("scatLL"+nomStr, Scatbins, ldr, kReMIdScatLLH, nom_cut);
209  dedxLL_nom [det][quant][gen] = new Spectrum("dEdxLL"+nomStr, dEdxbins, ldr, kReMIdDEDxLLH, nom_cut);
210  MF_nom [det][quant][gen] = new Spectrum("meas. frac."+nomStr, MFbins, ldr, kReMIdMeasFrac, nom_cut);
211  }
212  }
213  }
214 
215  // --- Do all the things!
216  ldr.Go();
217 
218  // --- Save all the things!
219  TFile *OutFile = new TFile(sOutFile.c_str(),"RECREATE");
220  for (unsigned int det=0; det<nDet; ++det) {
221  for (unsigned int quant=0; quant<nQ; ++quant) {
222  for (unsigned int gen=0; gen<nGENIE; ++gen) {
223  std::string nomStr = cutNames[det]+"_"+quantNames[quant]+"_"+GENIEStr[gen];
224  RmT_nom [det][quant][gen] -> SaveTo(OutFile, TString("RmT_nom_")+TString(nomStr));
225  RmToT_nom [det][quant][gen] -> SaveTo(OutFile, TString("RmToT_nom_")+TString(nomStr));
226  trueHadE_nom [det][quant][gen] -> SaveTo(OutFile, TString("trueHadE_nom_")+TString(nomStr));
227  recoHadE_nom [det][quant][gen] -> SaveTo(OutFile, TString("recoHadE_nom_")+TString(nomStr));
228  png_nom [det][quant][gen] -> SaveTo(OutFile, TString("png_nom_")+TString(nomStr));
229  trueNuE_nom [det][quant][gen] -> SaveTo(OutFile, TString("trueNuE_nom_")+TString(nomStr));
230  recoNuE_nom [det][quant][gen] -> SaveTo(OutFile, TString("recoNuE_nom_")+TString(nomStr));
231  trueMuonE_nom [det][quant][gen] -> SaveTo(OutFile, TString("trueMuonE_nom_")+TString(nomStr));
232  recoMuonE_nom [det][quant][gen] -> SaveTo(OutFile, TString("recoMuonE_nom_")+TString(nomStr));
233  remid_nom [det][quant][gen] -> SaveTo(OutFile, TString("remid_nom_")+TString(nomStr));
234  cvnprod3_nom [det][quant][gen] -> SaveTo(OutFile, TString("cvnprod3_nom_")+TString(nomStr));
235  cvn2017_nom [det][quant][gen] -> SaveTo(OutFile, TString("cvn2017_nom_")+TString(nomStr));
236  Vx_nom [det][quant][gen] -> SaveTo(OutFile, TString("Vx_nom_")+TString(nomStr));
237  Vy_nom [det][quant][gen] -> SaveTo(OutFile, TString("Vy_nom_")+TString(nomStr));
238  Vz_nom [det][quant][gen] -> SaveTo(OutFile, TString("Vz_nom_")+TString(nomStr));
239  trklen_nom [det][quant][gen] -> SaveTo(OutFile, TString("trklen_nom_")+TString(nomStr));
240  scatLL_nom [det][quant][gen] -> SaveTo(OutFile, TString("scatLL_nom_")+TString(nomStr));
241  dedxLL_nom [det][quant][gen] -> SaveTo(OutFile, TString("dedxLL_nom_")+TString(nomStr));
242  MF_nom [det][quant][gen] -> SaveTo(OutFile, TString("MF_nom_")+TString(nomStr));
243  }
244  }
245  }
246  OutFile->Close();
247 }
const Var kHadE
Definition: NumuVars.h:23
const Cut kIsQE
Definition: TruthCuts.cxx:104
caf::Proxy< size_t > npng
Definition: SRProxy.h:2038
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
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
const Var kReMIdScatLLH
Definition: NumuVars.cxx:555
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Binning zbins
Definition: NumuCCIncBins.h:21
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
void hyperon_nom_macro(bool IsFarDet=true)
const Var weight
const Cut kNumuContainND2017([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;++i){TVector3 start=sr->vtx.elastic.fuzzyk.png[i].shwlid.start;TVector3 stop=sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -180.0) return false;if(std::max(start.X(), stop.X()) > 180.0) return false;if(std::min(start.Y(), stop.Y())< -180.0) return false;if(std::max(start.Y(), stop.Y()) > 180.0) return false;if(std::min(start.Z(), stop.Z())< 20.0) return false;if(std::max(start.Z(), stop.Z()) > 1525.0) return false;}if(sr->trk.kalman.ntracks< 1) return false;for(unsigned int i=0;i< sr->trk.kalman.ntracks;++i){if(i==sr->trk.kalman.idxremid) continue;else if(sr->trk.kalman.tracks[i].start.Z() > 1275||sr->trk.kalman.tracks[i].stop.Z() > 1275) return false;}return(sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.firstplane > 1 &&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1100 &&(sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&sr->sel.contain.kalfwdcellnd > 5 &&sr->sel.contain.kalbakcellnd > 10);})
Definition: NumuCuts2017.h:11
const Cut kIsRes
Definition: TruthCuts.cxx:111
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const Binning xybins
Definition: NumuCCIncBins.h:18
T sqrt(T number)
Definition: d0nt_math.hpp:156
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
OStream cerr
Definition: OStream.cxx:7
const unsigned int nQ
Definition: hyperon_plot.C:61
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
const Cut kIsDIS
Definition: TruthCuts.cxx:118
void SetSpillCut(const SpillCut &cut)
const Var kReMIdMeasFrac
Definition: NumuVars.cxx:557
const unsigned int nDet
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
Definition: NumuVars.h:65
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
ifstream inFile
Definition: AnaPlotMaker.h:34
const Var kTrueE([](const caf::SRProxy *sr){assert(sr->mc.nnu==1);return sr->mc.nu[0].E;})
Definition: Vars.cxx:85
const Cut kNumuContainFD2017
Definition: NumuCuts2017.h:21
const unsigned int nGENIE
std::string getenv(std::string const &name)
const Var kCCE
Definition: NumuVars.h:21
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
std::vector< float > Spectrum
Definition: Constants.h:610
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
const Var kW
const SystShifts kNoShift
Definition: SystShifts.cxx:21
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kTrueMuonE([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return 0.f;if(sr->mc.nu[0].prim.empty()) return 0.f;if(std::abs(sr->mc.nu[0].prim[0].pdg)!=13) return 0.f;return float(sr->mc.nu[0].prim[0].p.E);})
Definition: NumuVars.h:107
const Cut kNumuPID2018([](const caf::SRProxy *sr){std::cout<< "ERROR::kNumuPID2018, cutting on both cvnProd3Train and cvn2017."<< " Neither branch exists anymore. Returning False."<< std::endl;abort();return false;})
Definition: NumuCuts2018.h:22
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
const std::vector< std::string > GENIEStr
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
const Cut kNumuCutFD2018
Definition: NumuCuts2018.h:39
const Var kReMIdDEDxLLH
Definition: NumuVars.cxx:556
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
TFile * OutFile
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
const Cut kNumuQuality
Definition: NumuCuts.h:18
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
#define W(x)
const Cut kIsCoh
Definition: TruthCuts.cxx:133
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
enum BeamMode string