caf_numu_fd_validation_data.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 void caf_numu_fd_validation_data(std::string kInputFileName, std::string kOutputFileName, bool oldCAFs = false){
3  std::cout << "Sorry, you must run in compiled mode" << std::endl;
4 }
5 #else
6 
7 #include "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/Cuts/NumuCuts.h"
12 #include "CAFAna/Core/Spectrum.h"
14 #include "CAFAna/Vars/Vars.h"
15 #include "CAFAna/Vars/NumuVars.h"
16 
17 #include "TH1.h"
18 #include "TFile.h"
19 
20 using namespace ana;
21 
22 // And some binning definitions
23 const Binning kXYBins = Binning::Simple(12,-9,9);
24 const Binning kZBins = Binning::Simple(12,0,60);
26 
27 // Cut on energy
28 const Cut kCCEshiftBelow5({"energy.numusimp.trkccE","sel.remid.bestidx","energy.mutrkE.E","hdr.ismc"},
29  [](const caf::StandardRecord* sr)
30  {
31  float muE = sr->energy.mutrkE[sr->sel.remid.bestidx].E;
32  float hadE = sr->energy.numusimp.trkccE - muE;
33  return ( (sr->hdr.ismc && sr->energy.numusimp.trkccE < 5.0)
34  || (!sr->hdr.ismc && muE+1.14f*hadE < 5.0) );
35  });
36 
37 // Cut for blinding
38 const Cut kBlindedRegion({"slc.calE", "slc.meantime", "hdr.run"},
39  [](const caf::StandardRecord* sr)
40  {
41  return sr->slc.calE > 0.8 && sr->slc.calE < 3.0
42  && ( (sr->slc.meantime > 207000 && sr->slc.meantime < 237000)
43  || (sr->slc.meantime > 271000 && sr->slc.meantime < 301000 && sr->hdr.run <= 17500) );
44  });
45 
46 const Cut kIsDytmanMEC({"mc.nnu", "mc.nu.mode"},
47  [](const caf::StandardRecord* sr)
48  {
49  if(sr->mc.nnu == 0) return false;
50  assert(sr->mc.nnu == 1);
51  return (sr->mc.nu[0].mode == caf::kMEC);
52  });
53 const Var kSlcMaxX({"slc.boxmax.fX"},
54  [](const caf::StandardRecord *sr)
55  {
56  return sr->slc.boxmax.X()/100;
57  });
58 const Var kSlcMaxY({"slc.boxmax.fY"},
59  [](const caf::StandardRecord *sr)
60  {
61  return sr->slc.boxmax.Y()/100;
62  });
63 const Var kSlcMaxZ({"slc.boxmax.fZ"},
64  [](const caf::StandardRecord *sr)
65  {
66  return sr->slc.boxmax.Z()/100;
67  });
68 
69 void caf_numu_fd_validation_data(std::string kInputFileName, std::string kOutputFileName, bool oldCAFs = false)
70 {
71 
72  std::cout << "\nrun : ---- Running CAF NuMU FD Validation - data.\n";
73  std::cout << "run : input file name: " << kInputFileName << std::endl;
74  std::cout << "run : output file name: " << kOutputFileName << std::endl;
75  std::cout << "run : old caf mode: " << oldCAFs << std::endl;
76 
77  SpectrumLoader loader(kInputFileName);
79 
80  struct Selection
81  {
83  Cut cut;
84  };
85 
86  const int kNumCuts = 10;
87 
88  Selection vCut[kNumCuts] =
89  {
90  {"nocut", kNoCut},
91  {"nocut_blinded", !kBlindedRegion},
92  {"DQ", kNumuQuality},
93  {"DQ_blinded", kNumuQuality && !kBlindedRegion},
94  {"DQ_containment",kNumuQuality && kNumuContainFD},
95  {"DQ_containment_blinded",kNumuQuality && kNumuContainFD && !kBlindedRegion},
96  {"DQ_containment_NCRej",kNumuQuality && kNumuContainFD && kNumuNCRej},
97  {"DQ_containment_NCRej_blinded",kNumuQuality && kNumuContainFD && kNumuNCRej && !kBlindedRegion},
98  {"all_numu_cuts", kNumuFD && kCCEshiftBelow5},
99  {"all_numu_cuts_blinded", kNumuFD && kCCEshiftBelow5 && !kBlindedRegion}
100  };
101 
102  struct Plot
103  {
106  Binning bins;
107  Var var;
108  };
109 
110  const int kNumPlots = 34;
111 
112  Plot plots[kNumPlots] = {
113  {";Reconstructed Neutrino Energy (GeV);Events / (0.25 GeV)", "ccEshift", kEnergyBinning, kCCEshift},
114  {";#Deltat from t_{0} (#mus);Events / (12 #mus)", "tus", Binning::Simple(35,26.125,446.125), kSliceTime},
115  {";Reconstructed Neutrino Energy (GeV);Events / (0.25 GeV)", "ccE", kEnergyBinning, kCCE},
116  {";Calorimetric Energy (GeV);Events / (0.25 GeV)", "calE", kEnergyBinning, SIMPLEVAR(slc.calE)},
117  {";Hadronic Energy (GeV);Events / (0.25 GeV)", "hadEshift", kEnergyBinning, kHadEshift},
118  {";Calorimetric Hadronic Energy (GeV);Events / (0.25 GeV)", "hadcalE", kEnergyBinning, kNumuHadCalE},
119  {";Slice Maximum x (m);Events / (1.5 m)", "maxX", kXYBins, kSlcMaxX},
120  {";Slice Maximum y (m);Events / (1.5 m)", "maxY", kXYBins, kSlcMaxY},
121  {";Slice Maximum z (m);Events / (5 m)", "maxZ", kZBins, kSlcMaxZ},
122  {";Number of Tracks in Slice;Events / bin", "nkal", Binning::Simple(12, 1, 13), kNKalman},
123  {";Total Number of Hits;Events / bin", "nhit", Binning::Simple(8, 0, 400), kNHit},
124  {";Cosmic Rejection BDT Output;Events / bin", "bdt", Binning::Simple(10, 0.53, 0.73), kNumuContPID},
125  {";Track Start X Position (m);Events / (1.5 m)", "trkStartX", kXYBins, kTrkStartX},
126  {";Track Start Y Position (m);Events / (1.5 m)", "trkStartY", kXYBins, kTrkStartY},
127  {";Track Start Z Position (m);Events / (5 m)", "trkStartZ", kZBins, kTrkStartZ},
128  {";Track End X Position (m);Events / (1.5 m)", "trkEndX", kXYBins, kTrkEndX},
129  {";Track End Y Position (m);Events / (1.5 m)", "trkEndY", kXYBins, kTrkEndY},
130  {";Track End Z Position (m);Events / (5 m)", "trkEndZ", kZBins, kTrkEndZ},
131  {";Number of Planes to Front;Events / bin", "nPlanesFront", Binning::Simple(15,0,900), kNplanesToFront},
132  {";Number of Planes to Back;Events / bin", "nPlanesBack", Binning::Simple(15,0,900), kNplanesToBack},
133  {";Number of Cells to Edge;Events / bin", "nCellsFromEdge", Binning::Simple(18,0,180), kNcellsFromEdge},
134  {";Basic Track Forward Distance to Edge (m);Events / (5 m)", "cosFwdDist", kZBins, kCosFwdDist},
135  {";Basic Track Backward Distance to Edge (m);Events / (5 m)", "cosBackDist", kZBins, kCosBackDist},
136  {";Kalman Track Forward Distance to Edge (m);Events / (5 m)", "kalFwdDist", kZBins, kKalmanFwdDist},
137  {";Kalman Track Backward Distance to Edge (m);Events / (5 m)", "kalBackDist", kZBins, kKalmanBackDist},
138  {";Slice Duration (ns);Events / (200 ns)", "sliceDuration", Binning::Simple(15,0,3000), kSliceDuration},
139  {";Number of Hits in Main Track;Events / bin", "trkNhits", Binning::Simple(20,0,500), kTrkNhits},
140  {";Number of Hits not in Main Track;Events / bin", "hadNhit", Binning::Simple(8, 0, 160), kHadNHit},
141  {";Length of Main Track (m);Events / (2 m)", "trkLength", Binning::Simple(15,0,30), kTrkLength},
142  {";Scattering Log-Likelihood;Events / bin", "scatLL", Binning::Simple(11,-0.2,0.3), kReMIdScatLLH},
143  {";dE/dx Log-Likelihood;Events / bin", "dedxLL", Binning::Simple(12,-0.1,0.5), kReMIdDEDxLLH},
144  {";cos#theta_{beam};Events / bin", "angkal", Binning::Simple(12, 0.3, 1), kNumuCosRejAngleKal},
145  {";cos#theta_{Z};Events / bin", "dirZ", Binning::Simple(12, 0.3, 1), kDirZ},
146  {";Muon ID;Events / bin", "remid", Binning::Simple(14,0.75-1./160,1.+1./160), kRemID}
147  };
148 
149  Spectrum* hData[kNumPlots][kNumCuts];
150 
151  for(int j = 0; j < kNumCuts; ++j){
152 
153  Selection s = vCut[j];
154 
155  for(int i = 0; i < kNumPlots; ++i){
156 
157  Plot p = plots[i];
158 
159  hData[i][j] = new Spectrum(p.label, p.bins, loader, p.var, oldCAFs? s.cut && kInBeamSpillShifted : s.cut && kInBeamSpill);
160 
161  }
162  }
163 
164  // GO!
165  std::cout << "\nrun : --- run loader.\n";
166  loader.Go();
167  std::cout << "\nrun : --- done.\n";
168 
169  double pot = hData[0][0]->POT();
170 
171  std::cout << "\nrun : --- save output.\n";
172  TFile* f = new TFile(kOutputFileName.c_str(),"RECREATE");
173 
174  for(int j = 0; j < kNumCuts; ++j){
175  for(int i = 0; i < kNumPlots; ++i){
176 
177  std::string fullName = "reco-"+plots[i].fname+"-"+vCut[j].fname;
178  std::cout << "run : save: " << fullName << std::endl;
179 
180  TH1* Data = hData[i][j]->ToTH1(pot);
181 
182  Data->SetName(fullName.c_str());
183 
184  Data->Write();
185  } // End loop over i
186  } // End loop over j
187 
188  TH1D* TotalPOT = new TH1D("meta-TotalPOT",";;Total POT", 1, 0, 1);
189  TotalPOT->SetBinContent(1, pot);
190  TotalPOT->Write();
191 
192  f->Write();
193  f->Close();
194 
195  std::cout << "\nrun : --- done.\n";
196 
197 } // End of function
198 
199 #endif
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
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
const Var kNcellsFromEdge
Definition: NumuVars.cxx:560
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kReMIdScatLLH
Definition: NumuVars.cxx:555
SRHeader hdr
Header branch: run, subrun, etc.
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
std::string label
Definition: CutFlow_Data.C:29
const Var kSlcMaxZ([](const caf::SRProxy *sr){return sr->slc.boxmax.Z()/100.;})
Definition: NumuVars.h:93
std::string fname
Definition: getTimePeak.C:84
const Var kSlcMaxY([](const caf::SRProxy *sr){return sr->slc.boxmax.Y()/100.;})
Definition: NumuVars.h:92
const Var kNumuContPID
Definition: NumuVars.cxx:553
SRVector3D boxmax
Maximum coordinates box containing all the hits [cm].
Definition: SRSlice.h:47
const Var kNumuCosRejAngleKal
Definition: NumuVars.cxx:547
const char * p
Definition: xmltok.h:285
float Y() const
Definition: SRVector3D.h:33
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
const Var kSliceTime([](const caf::SRProxy *sr){return sr->slc.meantime/1000;})
Definition: NumuVars.h:34
const Var kCosBackDist([](const caf::SRProxy *sr){return sr->sel.contain.cosbakdist/100.;})
Definition: NumuVars.h:71
unsigned int run
run number
Definition: SRHeader.h:21
const Var kSliceDuration([](const caf::SRProxy *sr){return(sr->slc.endtime-sr->slc.starttime);})
Definition: NumuVars.h:35
const Var kCosFwdDist([](const caf::SRProxy *sr){return sr->sel.contain.cosfwddist/100.;})
Definition: NumuVars.h:72
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
void SetSpillCut(const SpillCut &cut)
bool ismc
data or MC? True if MC
Definition: SRHeader.h:27
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 kSlcMaxX([](const caf::SRProxy *sr){return sr->slc.boxmax.X()/100.;})
Definition: NumuVars.h:91
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Binning kZBins
const Var kNumuHadCalE
Definition: NumuVars.cxx:538
const char * label
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
float X() const
Definition: SRVector3D.h:32
void caf_numu_fd_validation_data(std::string kInputFileName, std::string kOutputFileName, bool oldCAFs=false)
const XML_Char * s
Definition: expat.h:262
const Cut kCCEshiftBelow5({"energy.numusimp.trkccE","sel.remid.bestidx","energy.mutrkE.E","hdr.ismc"}, [](const caf::StandardRecord *sr){float muE=sr->energy.mutrkE[sr->sel.remid.bestidx].E;float hadE=sr->energy.numusimp.trkccE-muE;return((sr->hdr.ismc &&sr->energy.numusimp.trkccE< 5.0) ||(!sr->hdr.ismc &&muE+1.14f *hadE< 5.0));})
====================================================================== ///
Definition: CutFlow_Data.C:28
const Var kKalmanFwdDist([](const caf::SRProxy *sr){return sr->sel.contain.kalfwddist/100.;})
Definition: NumuVars.h:74
const Var kRemID
PID
Definition: Vars.cxx:81
const int kNumPlots
Definition: GetSpectra.h:1
const Var kNHit
Definition: Vars.cxx:71
const std::vector< Plot > plots
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
float calE
Calorimetric energy of the cluster [GeV].
Definition: SRSlice.h:38
const Var kCCE
Definition: NumuVars.h:21
const Binning kEnergyBinning
#define pot
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
const double j
Definition: BetheBloch.cxx:29
const Cut kNumuNCRej([](const caf::SRProxy *sr){return(sr->sel.remid.pid >0.75);})
Definition: NumuCuts.h:24
loader
Definition: demo0.py:10
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
const Var kKalmanBackDist([](const caf::SRProxy *sr){return sr->sel.contain.kalbakdist/100.;})
Definition: NumuVars.h:75
float Z() const
Definition: SRVector3D.h:34
std::vector< float > Spectrum
Definition: Constants.h:610
double POT() const
Definition: Spectrum.h:227
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
const Binning bins
Definition: NumuCC_CPiBin.h:8
std::string fname
Definition: CutFlow_Data.C:30
The StandardRecord is the primary top-level object in the Common Analysis File trees.
const Cut cut
Definition: exporter_fd.C:30
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.
const Var kNplanesToBack
Definition: NumuVars.cxx:559
const Cut kBlindedRegion({"slc.calE","slc.meantime","hdr.run"}, [](const caf::StandardRecord *sr){return sr->slc.calE > 0.8 &&sr->slc.calE< 3.0 &&((sr->slc.meantime > 207000 &&sr->slc.meantime< 237000) ||(sr->slc.meantime > 271000 &&sr->slc.meantime< 301000 &&sr->hdr.run<=17500));})
const Binning kXYBins
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
SRIDBranch sel
Selector (PID) branch.
const Cut kInBeamSpill([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return kInBeamSpill_main(sr);else return kInBeamSpill_main(sr)||kInBeamSpill_shifted(sr);}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return kInBeamSpill_main.Livetime(spill);else return kInBeamSpill_main.Livetime(spill)+kInBeamSpill_shifted.Livetime(spill);}, [](const caf::SRSpillProxy *spill) -> double{return spill->spillpot;})
Does the event fall inside the window we call the beam spill?
Definition: TimingCuts.h:8
assert(nhit_max >=nhit_nbins)
SRSlice slc
Slice branch: nhit, extents, time, etc.
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
const Var kReMIdDEDxLLH
Definition: NumuVars.cxx:556
const Var kNplanesToFront
Definition: NumuVars.cxx:558
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Binning bins
Definition: CutFlow_Data.C:31
const Cut kNumuQuality
Definition: NumuCuts.h:18
Var var
Definition: CutFlow_Data.C:32
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
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
float meantime
mean time, weighted by charge [ns]
Definition: SRSlice.h:41
enum BeamMode string