caf_numu_fd_validation_MC.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 void caf_numu_fd_validation_MC(std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kInputTauFileName, std::string kOutputFileName){
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 const double pot = 36e20;
23 
24 // And some binning definitions
25 const Binning kXYBins = Binning::Simple(12,-9,9);
26 const Binning kZBins = Binning::Simple(12,0,60);
28 
29 // Cut on energy
30 const Cut kCCEshiftBelow5({"energy.numusimp.trkccE","sel.remid.bestidx","energy.mutrkE.E","hdr.ismc"},
31  [](const caf::StandardRecord* sr)
32  {
33  float muE = sr->energy.mutrkE[sr->sel.remid.bestidx].E;
34  float hadE = sr->energy.numusimp.trkccE - muE;
35  return ( (sr->hdr.ismc && sr->energy.numusimp.trkccE < 5.0)
36  || (!sr->hdr.ismc && muE+1.14f*hadE < 5.0) );
37  });
38 
39 // Cut for blinding
40 const Cut kBlindedRegion({"slc.calE", "slc.meantime", "hdr.run"},
41  [](const caf::StandardRecord* sr)
42  {
43  return sr->slc.calE > 0.8 && sr->slc.calE < 3.0
44  && ( (sr->slc.meantime > 207000 && sr->slc.meantime < 237000)
45  || (sr->slc.meantime > 271000 && sr->slc.meantime < 301000 && sr->hdr.run <= 17500) );
46  });
47 
48 // Dumb oscillations
49 const Var kWOscDumb({"mc.nnu", "mc.nu.woscdumb"},
50  [](const caf::StandardRecord* sr)
51  {
52  if(sr->mc.nnu == 0) return 0.f;
53  return sr->mc.nu[0].woscdumb;
54  });
55 
56 const Cut kIsDytmanMEC({"mc.nnu", "mc.nu.mode"},
57  [](const caf::StandardRecord* sr)
58  {
59  if(sr->mc.nnu == 0) return false;
60  assert(sr->mc.nnu == 1);
61  return (sr->mc.nu[0].mode == caf::kMEC);
62  });
63 const Var kSlcMaxX({"slc.boxmax.fX"},
64  [](const caf::StandardRecord *sr)
65  {
66  return sr->slc.boxmax.X()/100;
67  });
68 const Var kSlcMaxY({"slc.boxmax.fY"},
69  [](const caf::StandardRecord *sr)
70  {
71  return sr->slc.boxmax.Y()/100;
72  });
73 const Var kSlcMaxZ({"slc.boxmax.fZ"},
74  [](const caf::StandardRecord *sr)
75  {
76  return sr->slc.boxmax.Z()/100;
77  });
78 
79 void caf_numu_fd_validation_MC(std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kInputTauFileName, std::string kOutputFileName)
80 {
81  std::cout << "\nrun : ---- Running CAF NuMU FD Validation - MC.\n";
82  std::cout << "run : input non-swap file name: " << kInputNonSwapFileName << std::endl;
83  std::cout << "run : input swap file name: " << kInputSwapFileName << std::endl;
84  std::cout << "run : input tau file name: " << kInputTauFileName << std::endl;
85  std::cout << "run : output file name: " << kOutputFileName << std::endl;
86 
87  SpectrumLoader loaderNonswap(kInputNonSwapFileName);
88  SpectrumLoader loaderSwap(kInputSwapFileName);
89  SpectrumLoader loaderTau(kInputTauFileName);
90 
91  loaderNonswap.SetSpillCut(kStandardSpillCuts);
92  loaderSwap.SetSpillCut(kStandardSpillCuts);
94 
95  struct Selection
96  {
98  Cut cut;
99  };
100 
101  const int kNumCuts = 10;
102 
103  Selection vCut[kNumCuts] =
104  {
105  {"nocut", (!kIsDytmanMEC) && kNoCut},
106  {"nocut_blinded", (!kIsDytmanMEC) && !kBlindedRegion},
107  {"DQ", (!kIsDytmanMEC) && kNumuQuality},
108  {"DQ_blinded", (!kIsDytmanMEC) && kNumuQuality && !kBlindedRegion},
109  {"DQ_containment", (!kIsDytmanMEC) && kNumuQuality && kNumuContainFD},
110  {"DQ_containment_blinded", (!kIsDytmanMEC) && kNumuQuality && kNumuContainFD && !kBlindedRegion},
111  {"DQ_containment_NCRej", (!kIsDytmanMEC) && kNumuQuality && kNumuContainFD && kNumuNCRej},
112  {"DQ_containment_NCRej_blinded", (!kIsDytmanMEC) && kNumuQuality && kNumuContainFD && kNumuNCRej && !kBlindedRegion},
113  {"all_numu_cuts", (!kIsDytmanMEC) && kNumuFD && kCCEshiftBelow5},
114  {"all_numu_cuts_blinded", (!kIsDytmanMEC) && kNumuFD && kCCEshiftBelow5 && !kBlindedRegion}
115  };
116 
117  struct Plot
118  {
121  Binning bins;
122  Var var;
123  };
124 
125  const int kNumPlots = 34;
126 
127  Plot plots[kNumPlots] = {
128  {";Reconstructed Neutrino Energy (GeV);Events / (0.25 GeV)", "ccEshift", kEnergyBinning, kCCEshift},
129  {";#Deltat from t_{0} (#mus);Events / (12 #mus)", "tus", Binning::Simple(35,26.125,446.125), kSliceTime},
130  {";Reconstructed Neutrino Energy (GeV);Events / (0.25 GeV)", "ccE", kEnergyBinning, kCCE},
131  {";Calorimetric Energy (GeV);Events / (0.25 GeV)", "calE", kEnergyBinning, SIMPLEVAR(slc.calE)},
132  {";Hadronic Energy (GeV);Events / (0.25 GeV)", "hadEshift", kEnergyBinning, kHadEshift},
133  {";Calorimetric Hadronic Energy (GeV);Events / (0.25 GeV)", "hadcalE", kEnergyBinning, kNumuHadCalE},
134  {";Slice Maximum x (m);Events / (1.5 m)", "maxX", kXYBins, kSlcMaxX},
135  {";Slice Maximum y (m);Events / (1.5 m)", "maxY", kXYBins, kSlcMaxY},
136  {";Slice Maximum z (m);Events / (5 m)", "maxZ", kZBins, kSlcMaxZ},
137  {";Number of Tracks in Slice;Events / bin", "nkal", Binning::Simple(12, 1, 13), kNKalman},
138  {";Total Number of Hits;Events / bin", "nhit", Binning::Simple(8, 0, 400), kNHit},
139  {";Cosmic Rejection BDT Output;Events / bin", "bdt", Binning::Simple(10, 0.53, 0.73), kNumuContPID},
140  {";Track Start X Position (m);Events / (1.5 m)", "trkStartX", kXYBins, kTrkStartX},
141  {";Track Start Y Position (m);Events / (1.5 m)", "trkStartY", kXYBins, kTrkStartY},
142  {";Track Start Z Position (m);Events / (5 m)", "trkStartZ", kZBins, kTrkStartZ},
143  {";Track End X Position (m);Events / (1.5 m)", "trkEndX", kXYBins, kTrkEndX},
144  {";Track End Y Position (m);Events / (1.5 m)", "trkEndY", kXYBins, kTrkEndY},
145  {";Track End Z Position (m);Events / (5 m)", "trkEndZ", kZBins, kTrkEndZ},
146  {";Number of Planes to Front;Events / bin", "nPlanesFront", Binning::Simple(15,0,900), kNplanesToFront},
147  {";Number of Planes to Back;Events / bin", "nPlanesBack", Binning::Simple(15,0,900), kNplanesToBack},
148  {";Number of Cells to Edge;Events / bin", "nCellsFromEdge", Binning::Simple(18,0,180), kNcellsFromEdge},
149  {";Basic Track Forward Distance to Edge (m);Events / (5 m)", "cosFwdDist", kZBins, kCosFwdDist},
150  {";Basic Track Backward Distance to Edge (m);Events / (5 m)", "cosBackDist", kZBins, kCosBackDist},
151  {";Kalman Track Forward Distance to Edge (m);Events / (5 m)", "kalFwdDist", kZBins, kKalmanFwdDist},
152  {";Kalman Track Backward Distance to Edge (m);Events / (5 m)", "kalBackDist", kZBins, kKalmanBackDist},
153  {";Slice Duration (ns);Events / (200 ns)", "sliceDuration", Binning::Simple(15,0,3000), kSliceDuration},
154  {";Number of Hits in Main Track;Events / bin", "trkNhits", Binning::Simple(20,0,500), kTrkNhits},
155  {";Number of Hits not in Main Track;Events / bin", "hadNhit", Binning::Simple(8, 0, 160), kHadNHit},
156  {";Length of Main Track (m);Events / (2 m)", "trkLength", Binning::Simple(15,0,30), kTrkLength},
157  {";Scattering Log-Likelihood;Events / bin", "scatLL", Binning::Simple(11,-0.2,0.3), kReMIdScatLLH},
158  {";dE/dx Log-Likelihood;Events / bin", "dedxLL", Binning::Simple(12,-0.1,0.5), kReMIdDEDxLLH},
159  {";cos#theta_{beam};Events / bin", "angkal", Binning::Simple(12, 0.3, 1), kNumuCosRejAngleKal},
160  {";cos#theta_{Z};Events / bin", "dirZ", Binning::Simple(12, 0.3, 1), kDirZ},
161  {";Muon ID;Events / bin", "remid", Binning::Simple(14,0.75-1./160,1.+1./160), kRemID}
162  };
163 
164  Spectrum* hNonswap[kNumPlots][kNumCuts];
165  Spectrum* hSwap[kNumPlots][kNumCuts];
166  Spectrum* hTau[kNumPlots][kNumCuts];
167 
168  for(int j = 0; j < kNumCuts; ++j){
169 
170  Selection s = vCut[j];
171 
172  for(int i = 0; i < kNumPlots; ++i){
173 
174  Plot p = plots[i];
175 
176  hNonswap[i][j] = new Spectrum(p.label, p.bins, loaderNonswap, p.var, s.cut, kNoShift, kWOscDumb);
177  hSwap[i][j] = new Spectrum(p.label, p.bins, loaderSwap, p.var, s.cut, kNoShift, kWOscDumb);
178  hTau[i][j] = new Spectrum(p.label, p.bins, loaderTau, p.var, s.cut, kNoShift, kWOscDumb);
179 
180  }
181  }
182 
183  // GO!
184  std::cout << "\nrun : --- run loaders.\n";
185  loaderNonswap.Go();
186  loaderSwap.Go();
187  loaderTau.Go();
188  std::cout << "\nrun : --- done.\n";
189 
190  std::cout << "\nrun : --- save output.\n";
191  TFile* f = new TFile(kOutputFileName.c_str(),"RECREATE");
192 
193  for(int j = 0; j < kNumCuts; ++j){
194  for(int i = 0; i < kNumPlots; ++i){
195 
196  std::string fullName = "reco-"+plots[i].fname+"-"+vCut[j].fname;
197 
198  TH1* NonSwap = hNonswap[i][j]->ToTH1(pot);
199  TH1* FluxSwap = hSwap[i][j]->ToTH1(pot);
200  TH1* Tau = hTau[i][j]->ToTH1(pot);
201 
202  NonSwap->SetName(fullName.c_str());
203  NonSwap->Add(FluxSwap);
204  NonSwap->Add(Tau);
205 
206  NonSwap->Write();
207 
208  } // End loop over i
209  } // End loop over j
210 
211  TH1D* TotalPOT = new TH1D("meta-TotalPOT",";;Total POT", 1, 0, 1);
212  TotalPOT->SetBinContent(1, pot);
213  TotalPOT->Write();
214 
215  f->Write();
216  f->Close();
217 
218  std::cout << "\nrun : --- done.\n";
219 
220 } // End of function
221 
222 #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 double pot
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 Binning kZBins
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 Var kNumuHadCalE
Definition: NumuVars.cxx:538
const char * label
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
void caf_numu_fd_validation_MC(std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kInputTauFileName, std::string kOutputFileName)
float X() const
Definition: SRVector3D.h:32
const XML_Char * s
Definition: expat.h:262
====================================================================== ///
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 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));})
const Var kCCE
Definition: NumuVars.h:21
const Binning kEnergyBinning
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
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
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 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 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
const Binning bins
Definition: NumuCC_CPiBin.h:8
const Var kWOscDumb([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return 0.f;return float(sr->mc.nu[0].woscdumb);})
Definition: TruthVars.h:10
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
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
SRIDBranch sel
Selector (PID) branch.
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
loaderSwap
Definition: demo4.py:9
const Binning kXYBins
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