GetHistsND.C
Go to the documentation of this file.
1 // Fills spectra for systematic error band
2 #ifdef __CINT__
3 void GetHistsND()
4 {
5  std::cout << "Sorry, you must run in compiled mode" << std::endl;
6 }
7 #else
8 
9 #include "CAFAna/Core/Binning.h"
10 #include "CAFAna/Core/Loaders.h"
11 #include "CAFAna/Core/Spectrum.h"
13 #include "CAFAna/Core/Var.h"
14 #include "CAFAna/Cuts/Cuts.h"
15 #include "CAFAna/Cuts/NumuCuts.h"
16 #include "CAFAna/Cuts/SpillCuts.h"
18 #include "CAFAna/Vars/NumuVars.h"
20 
21 #include "TCanvas.h"
22 #include "TFile.h"
23 #include "TH1.h"
24 #include "TH2.h"
25 #include "TStyle.h"
26 
27 using namespace ana;
28 
29 void GetHistsND()
30 {
31 
32  std::string nd_nonswap = "/pnfs/nova/persistent/production/concat/R16-03-03-prod2reco.d/prod_decaf_R16-03-03-prod2reco.d_nd_genie_nonswap_nogenierw_fhc_nova_v08_full_numu_contain_v1/*.root";
33  SpectrumLoader loader(nd_nonswap);
35 
36  struct Plot2D
37  {
39  std::string labelx;
40  Binning binsx;
41  Var varx;
42  std::string labely;
43  Binning binsy;
44  Var vary;
45  Cut cut;
46  };
47 
48  const Cut kTrueEbelow5({"mc.nnu","mc.nu.E"},
49  [](const caf::StandardRecord* sr)
50  {
51  if (sr->mc.nnu < 1) return false;
52  else return (sr->mc.nu[0].E <= 5.0);
53  });
54 
55  const Cut kHasTrueMuon({"energy.numusimp.mc.truegoodmuon","energy.numusimp.mc.truemuonE"},
56  [](const caf::StandardRecord* sr)
57  {
58  if (sr->energy.numusimp.mc.truegoodmuon < 1) return false;
59  else return (sr->energy.numusimp.mc.truemuonE > 0);
60  });
61 
62  const Cut kCut = kTrueEbelow5 && kIsNumuCC && kHasTrueMuon && kNumuND;
63 
64  //// ------ Now the binnings and variables ------ /////
65 
66  const Binning kMuonEnergyBinning = Binning::Simple(150,0,5);
67  const Binning kHadEnergyBinning = Binning::Simple(150,0,5);
68  const Binning kTrackLengthBinning = Binning::Simple(150, 0, 15);
69 
70  std::vector<double> hadronBins;
71  double hadronAxis = 0.0;
72  for(int i = 0; i < 116; ++i){
73  hadronBins.push_back(hadronAxis);
74  if (hadronAxis < 1.0){hadronAxis = hadronAxis + 0.01;}
75  else if (hadronAxis < 1.5){hadronAxis = hadronAxis + 0.05;}
76  else {hadronAxis = hadronAxis + 0.1;}
77  }
78 
79  const Binning kHadEBinning = Binning::Custom(hadronBins);
80 
81  const Var kTrueHadE = kTrueE - kTrueMuonE;
82 
83  const Var kTrueCatcherE = SIMPLEVAR(energy.numusimp.mc.truemuoncatcherE);
84 
85  const Var kActive1 = SIMPLEVAR(energy.numusimp.ndhadcalactE);
86  const Var kActive2 = SIMPLEVAR(energy.numusimp.ndhadtrkactE);
87 
88  const Var kTran1 = SIMPLEVAR(energy.numusimp.ndhadcaltranE);
89  const Var kTran2 = SIMPLEVAR(energy.numusimp.ndhadtrktranE);
90 
91  const Var kCatcher1 = SIMPLEVAR(energy.numusimp.ndhadcalcatE);
92  const Var kCatcher2 = SIMPLEVAR(energy.numusimp.ndhadtrkcatE);
93 
94  const Var kHadActive = kActive1 + kActive2;
95  const Var kHadTran = kTran1 + kTran2;
96  const Var kHadCatcher = kCatcher1 + kCatcher2;
97 
98  const Var kHadAll = kHadActive + kHadTran + kHadCatcher;
99 
100  const Var kTrkCalAct = SIMPLEVAR(energy.numusimp.ndtrkcalactE);
101  const Var kTrkCalTran = SIMPLEVAR(energy.numusimp.ndtrkcaltranE);
102  const Var kTrkCalCat = SIMPLEVAR(energy.numusimp.ndtrkcalcatE);
103 
104  const Var kTrkLenAct ({"energy.numusimp.ndtrklenact"},
105  [](const caf::StandardRecord* sr)
106  {
107  return (sr->energy.numusimp.ndtrklenact / 100); // in m
108  });
109 
110  const Var kTrkLenCat ({"energy.numusimp.ndtrklencat"},
111  [](const caf::StandardRecord* sr)
112  {
113  return (sr->energy.numusimp.ndtrklencat / 100); // in m
114  });
115 
116  const Cut kAllActive = (kTrkCalAct + kTrkCalTran > 0) && (kTrkCalCat == 0);
117  const Cut kAllCatcher = (kTrkCalAct + kTrkCalTran == 0) && (kTrkCalCat > 0);
118  const Cut kActiveAndCatcher = (kTrkCalAct + kTrkCalTran > 0) && (kTrkCalCat > 0);
119 
120  const int kNumPlots2D = 4;
121 
122  Plot2D plots2D[kNumPlots2D] = {
123 
124  {"MuonE_hist_active","Reco muon track length (m)", kTrackLengthBinning, kTrkLenAct, "True muon energy (GeV)" , kMuonEnergyBinning, kTrueMuonE, kCut && kAllActive},
125  {"MuonE_hist_catcher","Reco muon track length (m)", kTrackLengthBinning, kTrkLenCat, "True muon energy (GeV)" , kMuonEnergyBinning, kTrueMuonE, kCut && kAllCatcher},
126  {"MuonE_hist_activeAndCatcher","Reco muon track length (m)", kTrackLengthBinning, kTrkLenCat, "True muon energy in catcher (GeV)" , kMuonEnergyBinning, kTrueCatcherE, kCut && kActiveAndCatcher},
127  {"HadE_hist_DIS","Visible hadronic energy (GeV)", kHadEBinning, kHadAll, "True hadronic energy (GeV)" , kHadEnergyBinning, kTrueHadE, kCut}
128  };
129 
130  Spectrum* sSpect2D[kNumPlots2D];
131 
132  for(int i = 0; i < kNumPlots2D; i++)
133  {
134  Plot2D p = plots2D[i];
135  // sSpect2D[i] = new Spectrum(p.labelx, p.labely, loader, p.binsx, p.varx, p.binsy, p.vary, p.cut, kNoShift, kTuftsWeightCC * kPPFXFluxCVWgt ); // for miniproduction
136  sSpect2D[i] = new Spectrum(p.labelx, p.labely, loader, p.binsx, p.varx, p.binsy, p.vary, p.cut, kNoShift, kTuftsWeightCC );
137  } // loop over 1D plots
138 
139  loader.Go();
140 
141  TFile f("2DPlotsForFittingND.root","RECREATE");
142 
143  for(int i = 0; i < kNumPlots2D; i++)
144  {
145  Plot2D p = plots2D[i];
146  TH2* h = sSpect2D[i]->ToTH2(sSpect2D[i]->POT());
147  h->SetName(p.name.c_str());
148  h->Write();
149  } // loop over 1D plots
150 
151 }
152 
153 #endif
const XML_Char * name
Definition: expat.h:151
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
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 char * p
Definition: xmltok.h:285
void Plot2D(TH1 *total, std::map< std::string, TH1 * > systs, std::vector< std::string > to_plot, std::string basename)
void SetSpillCut(const SpillCut &cut)
const Binning kHadEnergyBinning
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Var kTrkLenCat([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f; if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat< 0) return 0.f; if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat > 0) return float(sr->trk.kalman.tracks[ibesttrk].lenincat/100.); if(sr->trk.kalman.tracks[ibesttrk].leninact< 0 && sr->trk.kalman.tracks[ibesttrk].lenincat > 0) return float((sr->trk.kalman.tracks[ibesttrk].leninact/100.) +(sr->trk.kalman.tracks[ibesttrk].lenincat/100.));return-1000.f;})
Definition: NumuCCIncVars.h:75
const Var kTrueE([](const caf::SRProxy *sr){assert(sr->mc.nnu==1);return sr->mc.nu[0].E;})
Definition: Vars.cxx:85
const Cut kNumuND
Definition: NumuCuts.h:55
double energy
Definition: plottest35.C:25
const Var kTrkLenAct([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f; if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat< 0) return float((sr->trk.kalman.tracks[ibesttrk].leninact/100.) +(sr->trk.kalman.tracks[ibesttrk].lenincat/100.)); if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat > 0) return float(sr->trk.kalman.tracks[ibesttrk].leninact/100.); if(sr->trk.kalman.tracks[ibesttrk].leninact< 0 && sr->trk.kalman.tracks[ibesttrk].lenincat > 0) return 0.f;return-1000.f;})
Definition: NumuCCIncVars.h:73
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
loader
Definition: demo0.py:10
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
static Binning Custom(const std::vector< double > &edges)
Definition: Binning.cxx:145
std::vector< float > Spectrum
Definition: Constants.h:759
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
The StandardRecord is the primary top-level object in the Common Analysis File trees.
const Cut cut
Definition: exporter_fd.C:30
const Var kTuftsWeightCC
Definition: XsecTunes.h:30
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 Var kTrueHadE
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
const Cut kCut
Definition: VarsAndCuts.h:36
const Var kTrueCatcherE
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
const Cut kActiveAndCatcher([](const caf::SRProxy *sr){int ibesttrk=muonid_classifier::kBestMuonTrack(sr);if(sr->trk.kalman.ntracks< 1||ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return false;const caf::SRKalmanTrackProxy &bestmuon=sr->trk.kalman.tracks[ibesttrk];return(bestmuon.leninact > 0 && bestmuon.lenincat > 0);})
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
void GetHistsND()
Definition: GetHistsND.C:29
SREnergyBranch energy
Energy estimator branch.
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
std::vector< SRNeutrino > nu
implemented as a vector to maintain mc.nu structure, i.e. not a pointer, but with 0 or 1 entries...
Definition: SRTruthBranch.h:25
enum BeamMode string