caf_numu_vars.C
Go to the documentation of this file.
1 // Uses truth information. Only run on MC
2 
3 /*
4 */
5 #ifdef __CINT__
6 void caf_numu_vars(std::string kInputFileName, std::string kOutputFileName, bool isFD)
7 {
8  std::cout << "Sorry, you must run in compiled mode" << std::endl;
9 }
10 #else
11 
12 #include "CAFAna/Core/Binning.h"
13 #include "CAFAna/Core/Spectrum.h"
16 
17 #include "CAFAna/Cuts/Cuts.h"
18 #include "CAFAna/Cuts/NumuCuts.h"
19 #include "CAFAna/Cuts/SpillCuts.h"
20 #include "CAFAna/Cuts/TruthCuts.h"
21 
22 #include "CAFAna/Vars/Vars.h"
23 #include "CAFAna/Vars/NumuVars.h"
24 #include "CAFAna/Vars/HistAxes.h"
25 
26 #include "CAFAna/Analysis/Plots.h"
27 #include "CAFAna/Analysis/Style.h"
28 
29 // #include "NumuEnergy/NumuEAlg.h"
30 
31 #include "TCanvas.h"
32 #include "TH1.h"
33 #include "TH2.h"
34 #include "TLine.h"
35 #include "TFile.h"
36 #include "TStyle.h"
37 #include "TLegend.h"
38 #include "TLine.h"
39 #include "TLatex.h"
40 #include "THStack.h"
41 
42 using namespace ana;
43 
44 // Put a "NOvA Preliminary" tag in the corner
45 void Preliminary(){
46  TLatex* prelim = new TLatex(.9, .95, "NO#nuA Preliminary");
47  prelim->SetTextColor(kBlue);
48  prelim->SetNDC();
49  prelim->SetTextSize(2/30.);
50  prelim->SetTextAlign(32);
51  prelim->Draw();
52 }
53 
54 template<class T> class Tangible{
55 
56 public:
57  Tangible(const T& obj, const std::string& shortName,
58  const std::string& blurb ):
59  fObj(obj),
60  fShortName(shortName),
61  fBlurb(blurb)
62  {};
63 
64  T fObj;
67 
68 };
69 
70 
73 
74 const Cut kIsDytmanMEC({"mc.nnu", "mc.nu.mode"},
75  [](const caf::StandardRecord* sr)
76  {
77  if(sr->mc.nnu == 0) return false;
78  assert(sr->mc.nnu == 1);
79  return (sr->mc.nu[0].mode == caf::kMEC);
80  });
81 
82 class UsefulHist{
83 
84 public:
85 
87  fSelA(selA),
88  fAxis(tanAxis),
89  fHist(loader, fAxis.fObj, fSelA.fObj && !kIsDytmanMEC),
90  fName(tanAxis.fShortName + "-" + selA.fShortName ){};
91 
92  Selection fSelA;
93  TangibleAxis fAxis;
94  Spectrum fHist;
95  std::string fName;
96 
97 };
98 
99 
100 void caf_numu_vars(std::string kInputFileName, std::string kOutputFileName, bool isFD){
101 
102  std::cout << "\nrun : ---- Running CAF NuMU ND Validation.\n";
103  std::cout << "run : input file name: " << kInputFileName << std::endl;
104  std::cout << "run : output file name: " << kOutputFileName << std::endl;
105  std::cout << "Are we using FD binning/cuts? " << isFD << std::endl;
106 
107  // make loader with user input dataset
108  SpectrumLoader loader(kInputFileName);
109 
110  // make spill histogram tokeep record of spills. Going to use to count total number of events
111  TH1D* spillHist = new TH1D("reco-spills", ";Detector;Spills", 3, 0, 3);
112  SpillVar spillRun({"det"}, [](const caf::SRSpill* spill) {return spill->det;});
113 
114  // set spill cuts
116  loader.AddSpillHistogram(spillHist, spillRun, kStandardSpillCuts);
117 
118 
119  //////////////////// Cuts ////////////////////////
120  std::vector<Selection> selections_base;
121 
122  if(isFD) selections_base.emplace_back(kNumuFD, "numuFD",
123  "Selected events pass the numu FD selection");
124  else selections_base.emplace_back(kNumuND, "numuND",
125  "Selected events pass the numu ND selection");
126  ////////////////////////////////////////////////////////////////////////////////
127 
128  // variables
129  const Var kRecoMuE({"energy.mutrkE.E","sel.remid.bestidx","energy.mutrkE"},
130  [](const caf::StandardRecord* sr)
131  {
132  if(sr->sel.remid.bestidx < 0.f) return 0.f;
133  if(sr->sel.remid.bestidx == 999) return 0.f;
134  if(sr->energy.mutrkE.size() == 0) return 0.f;
135  return (sr -> energy.mutrkE[sr->sel.remid.bestidx].E);
136  });
137 
138 
139  // Binning
140  const Binning kXYBins = Binning::Simple(55,-2.20,2.20); //ND
141  const Binning kZBins = Binning::Simple(32, 0,16.00); //ND
142  const Binning kFDXYBins = Binning::Simple(55,-8.0,8.0);
143  const Binning kFDZBins = Binning::Simple(55,0,60.00);
144  const Binning kEnergyBinning = Binning::Simple(50,0,5);
146 
147  // Axes
148  const HistAxis axRecoMuE("Muon E. (GeV)", kEnergyBinning, kRecoMuE);
149  /*
150  Variables not included yet:
151  kSliceTime,
152  kCosBackDist, kCosFwdDist, kKalmanFwdDist, kKalmanBackDist,
153  kNKalman, kNumuCosRejAngleKal, kNumuContPID, kNplanesToFront, kNplanesToBack, kNcellsFromEdge
154  */
155 
156  std::vector<TangibleAxis> variables;
157  variables.emplace_back(kNumuCCShiftAxis, "reco-numuEShift",
158  "Shifted CC energy estimator (GeV) ");
159  variables.emplace_back(kNumuCCAxis, "reco-numuE",
160  "CC energy estimator (GeV) ");
161  variables.emplace_back(axRecoMuE, "reco-muE",
162  "Muon energy (GeV)");
163  variables.emplace_back( HistAxis("Slice N_{Hit} (GeV)", Binning::Simple(50, 0, 500), kNHit),
164  "reco-slcNHit", "Number of hits in slice. " );
165  variables.emplace_back( HistAxis("Slice Calorimetric Energy (GeV)", kEnergyBinning, kCaloE),
166  "reco-calE", "Calorimetric energy of slice. " );
167  variables.emplace_back( HistAxis("Hadronic Energy (GeV)", kHadronicEnergyBinning, kHadE),
168  "reco-hadE",
169  "Hadronic enery, i.e. numu energy estimate minus muon track energy. " );
170  variables.emplace_back( HistAxis("Hadronic Energy (GeV)", kHadronicEnergyBinning, kHadEshift),
171  "reco-hadEShift",
172  "Hadronic enery, i.e. numu energy estimate minus muon track energy. " );
173  variables.emplace_back( HistAxis("Average Hadronic Energy Per Hit (GeV)",
174  Binning::Simple(40, 0, 0.04), kHadEPerNHit),
175  "reco-hadEPerNHit",
176  "Average energy per hit in hadronic cluster, i.e. had_E/had_n_hit. " );
177  variables.emplace_back( HistAxis("Track Energy Per Hit (GeV)", Binning::Simple(40, 0, 0.04),
178  kTrkEPerNHit),
179  "reco-trkEPerNHit",
180  "Average energy per hit on primary track, i.e. trk_E/trk_n_hit. " );
181  variables.emplace_back( HistAxis("Hadronic N_{Hit} (GeV)", Binning::Simple(50, 0, 100), kHadNHit),
182  "reco-hadNHit", "Number of hits in hadronic cluster. ");
183  variables.emplace_back( HistAxis("Calorimetric Hadronic Energy (GeV)",
184  kHadronicEnergyBinning, kNumuHadCalE),
185  "reco-hadcalE", "Sum of calibrated energy deposits in hadronic cluster. " );
186  variables.emplace_back( HistAxis("Slice Maximum Y (m)", kXYBins, kMaxY),
187  "reco-maxy", "Maximum Y position of slice. " );
188  variables.emplace_back( HistAxis("Number of Tracks in Slice", Binning::Simple(14, 1, 15), kNKalman),
189  "reco-nkal", "Number of tracks in slice. " );
190  variables.emplace_back( HistAxis("Number of Hits in Slice", Binning::Simple(50, 0, 500), kNHit),
191  "reco-nhit", "Number of hits in slice. " );
192  variables.emplace_back( HistAxis("Muon Track cos(#theta_{Z})", Binning::Simple(50, 0,1), kDirZ),
193  "reco-dirZ", "Z-direction of muon track. " );
194  variables.emplace_back( HistAxis("Track Start X Position (m)", kXYBins, kTrkStartX),
195  "reco-trkStartX", "Track start x position. " );
196  variables.emplace_back( HistAxis("Track Start Y Position (m)", kXYBins, kTrkStartY),
197  "reco-trkStartY", "Track start y position. " );
198  variables.emplace_back( HistAxis("Track Start Z Position (m)", kZBins, kTrkStartZ),
199  "reco-trkStartZ", "Track start z position. " );
200  variables.emplace_back( HistAxis("Track End X Position (m)", kXYBins, kTrkEndX),
201  "reco-trkEndX", "Track stop x position. " );
202  variables.emplace_back( HistAxis("Track End Y Position (m)", kXYBins, kTrkEndY),
203  "reco-trkEndY", "Track stop y position. " );
204  variables.emplace_back( HistAxis("Track End Z Position (m)", kZBins, kTrkEndZ),
205  "reco-trkEndZ", "Track stop z position. " );
206  variables.emplace_back( HistAxis("Slice Duration [ns]", Binning::Simple(50,0,3000), kSliceDuration),
207  "reco-sliceDuration", "Slice duration. " );
208  variables.emplace_back( HistAxis("Number of Hits in Primary Track", Binning::Simple(50,0,500), kTrkNhits),
209  "reco-trkNhits", "Number of hits on primary track. ");
210  variables.emplace_back( HistAxis("Length of Primary Track (m)", Binning::Simple(50,0,16), kTrkLength),
211  "reco-trkLength", "Primary track length. " );
212  variables.emplace_back( HistAxis("Scattering Log-likelihood", Binning::Simple(50,-0.5,0.5), kReMIdScatLLH),
213  "reco-scatLL", "ReMId scattering log log-likelihood for primary track. " );
214  variables.emplace_back( HistAxis("dE/dx Log-likelihood", Binning::Simple(50,-3,1), kReMIdDEDxLLH),
215  "reco-dedxLL", "ReMId dE/dx log log-likelihood for primary track. " );
216  variables.emplace_back( HistAxis("Non-hadronic Plane Fraction",
218  "reco-nonHadPlaneFrac", "ReMId Non-hadronic plane fraction. " );
219  variables.emplace_back( HistAxis("Muon ID", kRemidBinning, kRemID),
220  "reco-remid", "ReMId kNN score. " );
221  ////////////////////////////////////////////////////////////////////////////////
222 
223  std::vector<UsefulHist> hists;
224  hists.reserve(selections_base.size() * variables.size());
225 
226  for(const auto& sel_base:selections_base){
227  for(const auto& variable:variables){
228  hists.emplace_back(sel_base, variable, loader);
229  }
230  }
231 
232 
233  std::cout << "\nrun : --- run loader.\n";
234  loader.Go();
235  std::cout << "\nrun : --- done.\n";
236  double pot = hists[0].fHist.POT();
237 
238  std::cout << "\nrun : --- save output.\n";
239  // TFile inputfile(kInputFileName.c_str(), "READ");
240  TFile outputfile(kOutputFileName.c_str(), "RECREATE");
241 
242  for(const auto& hist:hists){
243  TH1 *temp = hist.fHist.ToTH1(pot);
244  std::string tempName = hist.fName;
245  temp->Write(tempName.c_str());
246  }
247 
248  std::cout << "\nrun : --- save PoT.\n";
249 
250  TH1F *pothist = new TH1F("meta-TotalPOT", "TotalPOT", 1, 0, 1);
251  TH1F *evthist = new TH1F("meta-TotalEvents", "TotalEvents", 1, 0, 1);
252  pothist->SetBinContent(1, pot);
253  evthist->SetBinContent(1, spillHist->GetEntries());
254  pothist->Write("meta-TotalPOT");
255  evthist->Write("meta-TotalEvents");
256 
257  std::cout << "\nrun : --- close output files.\n";
258  outputfile.Close();
259  // inputfile.Close(); // will produce error if running over SAM dataset
260  std::cout << "\nrun : --- done.\n";
261 
262  // dictionary to translate histogram names
263  /*
264 
265  */
266 
267  return;
268 }
269 
270 #endif
const Var kHadE
Definition: NumuVars.h:23
const Var kHadNHit([](const caf::SRProxy *sr){unsigned int nought=0;if(sr->trk.kalman.ntracks< 1) return nought;return sr->slc.nhit-sr->trk.kalman.tracks[0].nhit;})
Definition: NumuVars.h:61
const Var kNKalman
Definition: NumuVars.cxx:540
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
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
UsefulHist(Selection selA, TangibleAxis tanAxis, SpectrumLoader &loader, const Cut &bkg=kNoCut)
Definition: caf_numu_vars.C:86
const Binning kRemidBinning
Binning for plotting remid attractively.
Definition: Binning.cxx:80
const Var kTrkStartY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Y()/100;})
Definition: NumuVars.h:52
std::string fShortName
Definition: DataMCPair.h:47
Tangible< Cut > Selection
Definition: DataMCPair.h:52
const Var kSliceDuration([](const caf::SRProxy *sr){return(sr->slc.endtime-sr->slc.starttime);})
Definition: NumuVars.h:35
const Cut kNumuFD
Definition: NumuCuts.h:53
const Var kRecoMuE
void SetSpillCut(const SpillCut &cut)
const Var kTrkEPerNHit([](const caf::SRProxy *sr){return(sr->trk.kalman.tracks[0].calE/sr->trk.kalman.tracks[0].nhit);})
TString hists[nhists]
Definition: bdt_com.C:3
const Var kReMIdMeasFrac
Definition: NumuVars.cxx:557
const Var kDirZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].dir.Z();})
Definition: NumuVars.h:39
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
const Var kHadEPerNHit([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 0.0f;int nHit=sr->slc.nhit-sr->trk.kalman.tracks[0].nhit;if(nHit<=0) return 0.0f;float hadE=sr->energy.numu.hadcalE;return hadE/nHit;})
Definition: NumuVars.h:63
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Var kNumuHadCalE
Definition: NumuVars.cxx:538
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
void Preliminary()
Put NOvA Preliminary on plots.
Definition: caf_numu_vars.C:45
const Var kMaxY([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1000.0f;if(sr->vtx.elastic.fuzzyk.nshwlid==0) return-1000.0f;float maxyall=-999.0;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;i++){maxyall=std::max(std::max(sr->vtx.elastic.fuzzyk.png[i].shwlid.start.Y(), sr->vtx.elastic.fuzzyk.png[i].shwlid.stop.Y()), maxyall);}return maxyall;})
Definition: NueVars.h:60
const Var kRemID
PID
Definition: Vars.cxx:81
const HistAxis kNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNumuEnergyBinning, kCCE)
Definition: HistAxes.h:8
virtual void AddSpillHistogram(TH1 *h, const SpillVar &var, const SpillCut &cut, const SpillVar &wei=kSpillUnweighted)
Uses include counting the total POT or spills in a run.
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
Tangible(const T &obj, const std::string &shortName, const std::string &blurb)
Definition: caf_numu_vars.C:57
const Binning kHadronicEnergyBinning
const Var kNHit
Definition: Vars.cxx:71
const Var kTrkStartZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Z()/100;})
Definition: NumuVars.h:53
SRRemid remid
Output from RecoMuonID (ReMId) package.
Definition: SRIDBranch.h:39
Tangible< HistAxis > TangibleAxis
Definition: DataMCPair.h:53
const Cut kNumuND
Definition: NumuCuts.h:55
const Binning kEnergyBinning
double energy
Definition: plottest35.C:25
#define pot
const Binning kXYBins
Definition: VarsAndCuts.h:95
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
TLatex * prelim
Definition: Xsec_final.C:133
const Var kTrkStartX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.X()/100;})
Definition: NumuVars.h:51
const Var kTrkEndZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.Z()/100;})
Definition: NumuVars.h:57
OStream cout
Definition: OStream.cxx:6
std::string fBlurb
Definition: DataMCPair.h:46
The StandardRecord is the primary top-level object in the Common Analysis File trees.
const Var kTrkEndY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.Y()/100;})
Definition: NumuVars.h:56
const Var kTrkEndX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.X()/100;})
Definition: NumuVars.h:55
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
SRIDBranch sel
Selector (PID) branch.
assert(nhit_max >=nhit_nbins)
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
const Var kReMIdDEDxLLH
Definition: NumuVars.cxx:556
string tempName
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
double T
Definition: Xdiff_gwt.C:5
void caf_numu_vars(std::string kInputFileName, std::string kOutputFileName, bool isFD)
Template for Var and SpillVar.
SREnergyBranch energy
Energy estimator branch.
const Binning kZBins
Definition: VarsAndCuts.h:96
enum BeamMode kBlue
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
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
const Var kTrkNhits([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 65535;return int(sr->trk.kalman.tracks[0].nhit);})
Definition: NumuVars.h:59
enum BeamMode string