numu_validation_numuvars.C
Go to the documentation of this file.
1 // Uses truth information. Only run on MC
2 
3 /*
4 */
5 
6 #include "CAFAna/Core/Binning.h"
7 #include "CAFAna/Core/Spectrum.h"
10 
11 #include "CAFAna/Cuts/Cuts.h"
13 #include "CAFAna/Cuts/SpillCuts.h"
14 #include "CAFAna/Cuts/TruthCuts.h"
15 
16 #include "CAFAna/Vars/Vars.h"
20 
21 #include "CAFAna/Analysis/Plots.h"
22 #include "CAFAna/Analysis/Style.h"
23 
24 // #include "NumuEnergy/NumuEAlg.h"
25 
26 #include "TCanvas.h"
27 #include "TH1.h"
28 #include "TH2.h"
29 #include "TLine.h"
30 #include "TFile.h"
31 #include "TStyle.h"
32 #include "TLegend.h"
33 #include "TLine.h"
34 #include "TLatex.h"
35 #include "THStack.h"
36 
37 using namespace ana;
38 
39 // Put a "NOvA Preliminary" tag in the corner
40 void Preliminary(){
41  TLatex* prelim = new TLatex(.9, .95, "NO#nuA Preliminary");
42  prelim->SetTextColor(kBlue);
43  prelim->SetNDC();
44  prelim->SetTextSize(2/30.);
45  prelim->SetTextAlign(32);
46  prelim->Draw();
47 }
48 
49 template<class T> class Tangible{
50 
51 public:
52  Tangible(const T& obj, const std::string& shortName,
53  const std::string& blurb ):
54  fObj(obj),
55  fShortName(shortName),
56  fBlurb(blurb)
57  {};
58 
59  T fObj;
62 
63 };
64 
65 
68 
69 class UsefulHist{
70 
71 public:
72 
73  UsefulHist(Selection selA, TangibleAxis tanAxis, SpectrumLoader& loader, const Cut& bkg=kNoCut, const SystShifts &shift=kNoShift, const Var& wei = kTuftsWeightCC):
74  fSelA(selA),
75  fAxis(tanAxis),
76  fHist(loader, fAxis.fObj, fSelA.fObj, shift, wei),
77  fName(tanAxis.fShortName + "-" + selA.fShortName ){};
78 
80  TangibleAxis fAxis;
81  Spectrum fHist;
82  std::string fName;
83 
84 };
85 
86 
87 void numu_validation_numuvars(std::string kInputFileName, std::string kOutputFileName, bool isFD){
88 
89  std::cout << "\nrun : ---- Running CAF NuMU ND Validation.\n";
90  std::cout << "run : input file name: " << kInputFileName << std::endl;
91  std::cout << "run : output file name: " << kOutputFileName << std::endl;
92  std::cout << "Are we using FD binning/cuts? " << isFD << std::endl;
93 
94  // make loader with user input dataset
95  SpectrumLoader loader(kInputFileName);
96 
97  // make spill histogram tokeep record of spills. Going to use to count total number of events
98  TH1D* spillHist = new TH1D("reco-spills", ";Detector;Spills", 3, 0, 3);
99  SpillVar spillRun([](const caf::SRSpill* spill) {return spill->det;});
100 
101  // set spill cuts
103  loader.AddSpillHistogram(spillHist, spillRun, kStandardSpillCuts);
104 
105 
106  //////////////////// Cuts ////////////////////////
107  std::vector<Selection> selections_base;
108 
109  const Cut kNumuFD_cvn = kNumuQuality && kNumuContainFD && kCVNmSel && kNumuCosmicRej;
110  const Cut kNumuND_cvn = kNumuQuality && kNumuContainND && kCVNmSel;
111 
112  if(isFD){
113  selections_base.emplace_back(kNumuFD, "numuFD_remid",
114  "Selected events pass the (remid) numu FD selection");
115  selections_base.emplace_back(kNumuFD_cvn, "numuFD_cvn",
116  "Selected events pass the (cvn) numu FD selection");
117  }
118  else if(!isFD){
119  selections_base.emplace_back(kNumuND, "numuND_remid",
120  "Selected events pass the (remid) numu ND selection");
121  selections_base.emplace_back(kNumuND_cvn, "numuND_cvn",
122  "Selected events pass the (cvn) numu ND selection");
123  }
124 
125  ////////////////////////////////////////////////////////////////////////////////
126 
127  // variables
128  const Var kRecoMuE([](const caf::SRProxy* sr)
129  {
130  return (sr -> energy.numu.recomuonE);
131  });
132 
133  // Binning
134  const Binning kXYBins = Binning::Simple(55,-2.20,2.20); //ND
135  const Binning kZBins = Binning::Simple(32, 0,16.00); //ND
136  const Binning kFDXYBins = Binning::Simple(55,-8.0,8.0);
137  const Binning kFDZBins = Binning::Simple(55,0,60.00);
138  const Binning kEnergyBinning = Binning::Simple(50,0,5);
140 
141  std::vector<TangibleAxis> variables;
142  variables.emplace_back(kNumuCCAxis, "numuE",
143  "CC energy estimator (GeV) ");
144  variables.emplace_back( HistAxis("Muon energy (GeV)", kEnergyBinning, kRecoMuE),
145  "muE", "Muon energy (GeV)");
146  variables.emplace_back( HistAxis("Slice N_{Hit} (GeV)", Binning::Simple(50, 0, 500), kNHit),
147  "slcNHit", "Number of hits in slice. " );
148  variables.emplace_back( HistAxis("Slice Calorimetric Energy (GeV)", kEnergyBinning, kCaloE),
149  "calE", "Calorimetric energy of slice. " );
150  variables.emplace_back( HistAxis("Hadronic Energy (GeV)", kHadronicEnergyBinning, kHadE),
151  "hadE",
152  "Hadronic enery, i.e. numu energy estimate minus muon track energy. " );
153  variables.emplace_back( HistAxis("Average Hadronic Energy Per Hit (GeV)",
154  Binning::Simple(40, 0, 0.04), kHadEPerNHit),
155  "hadEPerNHit",
156  "Average energy per hit in hadronic cluster, i.e. had_E/had_n_hit. " );
157  variables.emplace_back( HistAxis("Track Energy Per Hit (GeV)", Binning::Simple(40, 0, 0.04),
158  kTrkEPerNHit),
159  "trkEPerNHit",
160  "Average energy per hit on primary track, i.e. trk_E/trk_n_hit. " );
161  variables.emplace_back( HistAxis("Hadronic N_{Hit} (GeV)", Binning::Simple(50, 0, 100), kHadNHit),
162  "hadNHit", "Number of hits in hadronic cluster. ");
163  variables.emplace_back( HistAxis("Calorimetric Hadronic Energy (GeV)",
164  kHadronicEnergyBinning, kNumuHadCalE),
165  "hadcalE", "Sum of calibrated energy deposits in hadronic cluster. " );
166  // variables.emplace_back( HistAxis("Slice Maximum Y (m)", kXYBins, kMaxY), //not in decafs
167  // "maxy", "Maximum Y position of slice. " );
168  variables.emplace_back( HistAxis("Number of Tracks in Slice", Binning::Simple(14, 1, 15), kNKalman),
169  "nkal", "Number of tracks in slice. " );
170  variables.emplace_back( HistAxis("Number of Hits in Slice", Binning::Simple(50, 0, 500), kNHit),
171  "nhit", "Number of hits in slice. " );
172  variables.emplace_back( HistAxis("Muon Track cos(#theta_{Z})", Binning::Simple(50, 0,1), kDirZ),
173  "dirZ", "Z-direction of muon track. " );
174  variables.emplace_back( HistAxis("Slice Duration [ns]", Binning::Simple(50,0,3000), kSliceDuration),
175  "sliceDuration", "Slice duration. " );
176  variables.emplace_back( HistAxis("Number of Hits in Primary Track", Binning::Simple(50,0,500), kTrkNhits),
177  "trkNhits", "Number of hits on primary track. ");
178  variables.emplace_back( HistAxis("Length of Primary Track (m)", Binning::Simple(50,0,16), kTrkLength),
179  "trkLength", "Primary track length. " );
180  variables.emplace_back( HistAxis("Scattering Log-likelihood", Binning::Simple(50,-0.5,0.5), kReMIdScatLLH),
181  "scatLL", "ReMId scattering log log-likelihood for primary track. " );
182  variables.emplace_back( HistAxis("dE/dx Log-likelihood", Binning::Simple(50,-3,1), kReMIdDEDxLLH),
183  "dedxLL", "ReMId dE/dx log log-likelihood for primary track. " );
184  variables.emplace_back( HistAxis("Non-hadronic Plane Fraction",
186  "nonHadPlaneFrac", "ReMId Non-hadronic plane fraction. " );
187  variables.emplace_back( HistAxis("Muon ID", kRemidBinning, kRemID),
188  "remid", "ReMId kNN score. " );
189 
190  // detector specific variables:
191  if(!isFD){
192  variables.emplace_back( HistAxis("Track Start X Position (m)", kXYBins, kTrkStartX),
193  "trkStartX", "Track start x position. " );
194  variables.emplace_back( HistAxis("Track Start Y Position (m)", kXYBins, kTrkStartY),
195  "trkStartY", "Track start y position. " );
196  variables.emplace_back( HistAxis("Track Start Z Position (m)", kZBins, kTrkStartZ),
197  "trkStartZ", "Track start z position. " );
198  variables.emplace_back( HistAxis("Track End X Position (m)", kXYBins, kTrkEndX),
199  "trkEndX", "Track stop x position. " );
200  variables.emplace_back( HistAxis("Track End Y Position (m)", kXYBins, kTrkEndY),
201  "trkEndY", "Track stop y position. " );
202  variables.emplace_back( HistAxis("Track End Z Position (m)", kZBins, kTrkEndZ),
203  "trkEndZ", "Track stop z position. " );
204  }
205  else if(isFD){
206  variables.emplace_back( HistAxis("Track Start X Position (m)", kFDXYBins, kTrkStartX),
207  "trkStartX", "Track start x position. " );
208  variables.emplace_back( HistAxis("Track Start Y Position (m)", kFDXYBins, kTrkStartY),
209  "trkStartY", "Track start y position. " );
210  variables.emplace_back( HistAxis("Track Start Z Position (m)", kFDZBins, kTrkStartZ),
211  "trkStartZ", "Track start z position. " );
212  variables.emplace_back( HistAxis("Track End X Position (m)", kFDXYBins, kTrkEndX),
213  "trkEndX", "Track stop x position. " );
214  variables.emplace_back( HistAxis("Track End Y Position (m)", kFDXYBins, kTrkEndY),
215  "trkEndY", "Track stop y position. " );
216  variables.emplace_back( HistAxis("Track End Z Position (m)", kFDZBins, kTrkEndZ),
217  "trkEndZ", "Track stop z position. " );
218  }
219 
220  ////////////////////////////////////////////////////////////////////////////////
221 
222  std::vector<UsefulHist> hists;
223  hists.reserve(selections_base.size() * variables.size());
224 
225  for(const auto& sel_base:selections_base){
226  for(const auto& variable:variables){
227  hists.emplace_back(sel_base, variable, loader);
228  }
229  }
230 
231 
232  std::cout << "\nrun : --- run loader.\n";
233  loader.Go();
234  std::cout << "\nrun : --- done.\n";
235  double pot = hists[0].fHist.POT();
236 
237  std::cout << "\nrun : --- save output.\n";
238  TFile inputfile(kInputFileName.c_str(), "READ");
239  TFile outputfile(kOutputFileName.c_str(), "RECREATE");
240 
241  for(const auto& hist:hists){
242  TH1 *temp = hist.fHist.ToTH1(pot);
243  std::string tempName = hist.fName;
244  temp->Write(tempName.c_str());
245  }
246 
247  std::cout << "\nrun : --- save PoT.\n";
248 
249  TH1F *pothist = new TH1F("meta-TotalPOT", "TotalPOT", 1, 0, 1);
250  TH1F *evthist = new TH1F("meta-TotalEvents", "TotalEvents", 1, 0, 1);
251  pothist->SetBinContent(1, pot);
252  evthist->SetBinContent(1, spillHist->GetEntries());
253  pothist->Write("meta-TotalPOT");
254  evthist->Write("meta-TotalEvents");
255 
256  std::cout << "\nrun : --- close output files.\n";
257  outputfile.Close();
258  inputfile.Close(); // will produce error if running over SAM dataset
259  std::cout << "\nrun : --- done.\n";
260 
261  // dictionary to translate histogram names
262  /*
263 
264  */
265 
266  return;
267 }
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
const Binning kRemidBinning
Binning for plotting remid attractively.
Definition: Binning.cxx:80
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const Cut kCVNmSel([](const caf::SRProxy *sr){return(sr->sel.cvn.numuid >0.5);})
Definition: NumuCuts.h:26
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 Cut kNumuCosmicRej([](const caf::SRProxy *sr){return(sr->sel.cosrej.anglekal > 0.5 && sr->sel.cosrej.numucontpid2019 > 0.535 && sr->slc.nhit< 400);})
Definition: NumuCuts.h:32
const Var kSliceDuration([](const caf::SRProxy *sr){return(sr->slc.endtime-sr->slc.starttime);})
Definition: NumuVars.h:35
const Cut kNumuContainFD([](const caf::SRProxy *sr){ std::pair< int, int > planes=calcFirstLastLivePlane(sr->slc.firstplane, std::bitset< 14 >(sr->hdr.dibmask));int planestofront=sr->slc.firstplane-planes.first;int planestoback=planes.second-sr->slc.lastplane;return( sr->slc.ncellsfromedge > 1 &&planestofront > 1 &&planestoback > 1 &&sr->sel.contain.kalfwdcell > 10 &&sr->sel.contain.kalbakcell > 10 &&sr->sel.contain.cosfwdcell > 0 &&sr->sel.contain.cosbakcell > 0);})
Definition: NumuCuts.h:20
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
Det_t det
Detector, ND = 1, FD = 2, NDOS = 3.
Definition: SRSpill.h:29
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
const Cut kNumuContainND([](const caf::SRProxy *sr){return( sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.ncellsfromedge > 1 &&sr->slc.firstplane > 1 &&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1150 &&( sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&( sr->energy.numu.ndhadcalcatE +sr->energy.numu.ndhadcaltranE)< 0.03 &&sr->sel.contain.kalfwdcellnd > 4 &&sr->sel.contain.kalbakcellnd > 8);})
Definition: NumuCuts.h:22
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Var kNumuHadCalE
Definition: NumuVars.cxx:538
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)
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
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
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
void Preliminary()
const SystShifts kNoShift
Definition: SystShifts.cxx:21
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
UsefulHist(Selection selA, TangibleAxis tanAxis, SpectrumLoader &loader, const Cut &bkg=kNoCut, const SystShifts &shift=kNoShift, const Var &wei=kTuftsWeightCC)
std::string fBlurb
Definition: DataMCPair.h:46
void numu_validation_numuvars(std::string kInputFileName, std::string kOutputFileName, bool isFD)
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 kTuftsWeightCC
Definition: XsecTunes.h:31
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.
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
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
const Cut kNumuQuality
Definition: NumuCuts.h:18
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
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