caf_numu_fd_validation_MC_no_tau.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 void caf_numu_fd_validation_MC_no_tau(std::string kInputNonSwapFileName, std::string kInputSwapFileName, 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/TruthCuts.h"
11 #include "CAFAna/Cuts/NumuCuts.h"
13 #include "CAFAna/Core/Spectrum.h"
15 #include "CAFAna/Vars/Vars.h"
16 #include "CAFAna/Vars/NumuVars.h"
17 
18 #include "TH1.h"
19 #include "TFile.h"
20 
21 using namespace ana;
22 
23 const double pot = 36e20;
24 
25 // And some binning definitions
26 const Binning kXYBins = Binning::Simple(12,-9,9);
27 const Binning kZBins = Binning::Simple(12,0,60);
29 
30 // Cut on energy
31 const Cut kCCEshiftBelow5({"energy.numusimp.trkccE","sel.remid.bestidx","energy.mutrkE.E","hdr.ismc"},
32  [](const caf::StandardRecord* sr)
33  {
34  float muE = sr->energy.mutrkE[sr->sel.remid.bestidx].E;
35  float hadE = sr->energy.numusimp.trkccE - muE;
36  return ( (sr->hdr.ismc && sr->energy.numusimp.trkccE < 5.0)
37  || (!sr->hdr.ismc && muE+1.14f*hadE < 5.0) );
38  });
39 
40 // Cut for blinding
41 const Cut kBlindedRegion({"slc.calE", "slc.meantime", "hdr.run"},
42  [](const caf::StandardRecord* sr)
43  {
44  return sr->slc.calE > 0.8 && sr->slc.calE < 3.0
45  && ( (sr->slc.meantime > 207000 && sr->slc.meantime < 237000)
46  || (sr->slc.meantime > 271000 && sr->slc.meantime < 301000 && sr->hdr.run <= 17500) );
47  });
48 
49 // Dumb oscillations
50 const Var kWOscDumb({"mc.nnu", "mc.nu.woscdumb"},
51  [](const caf::StandardRecord* sr)
52  {
53  if(sr->mc.nnu == 0) return 0.f;
54  return sr->mc.nu[0].woscdumb;
55  });
56 
57 const Cut kIsDytmanMEC({"mc.nnu", "mc.nu.mode"},
58  [](const caf::StandardRecord* sr)
59  {
60  if(sr->mc.nnu == 0) return false;
61  assert(sr->mc.nnu == 1);
62  return (sr->mc.nu[0].mode == caf::kMEC);
63  });
64 
65 const Var kSlcMaxX({"slc.boxmax.fX"},
66  [](const caf::StandardRecord *sr)
67  {
68  return sr->slc.boxmax.X()/100;
69  });
70 const Var kSlcMaxY({"slc.boxmax.fY"},
71  [](const caf::StandardRecord *sr)
72  {
73  return sr->slc.boxmax.Y()/100;
74  });
75 const Var kSlcMaxZ({"slc.boxmax.fZ"},
76  [](const caf::StandardRecord *sr)
77  {
78  return sr->slc.boxmax.Z()/100;
79  });
80 
81 void caf_numu_fd_validation_MC_no_tau(std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kOutputFileName)
82 {
83  std::cout << "\nrun : ---- Running CAF NuMU FD Validation - MC (no Tau).\n";
84  std::cout << "run : input non-swap file name: " << kInputNonSwapFileName << std::endl;
85  std::cout << "run : input swap file name: " << kInputSwapFileName << std::endl;
86  std::cout << "run : output file name: " << kOutputFileName << std::endl;
87 
88  SpectrumLoader loaderNonswap(kInputNonSwapFileName);
89  SpectrumLoader loaderSwap(kInputSwapFileName);
90 
91  loaderNonswap.SetSpillCut(kStandardSpillCuts);
92  loaderSwap.SetSpillCut(kStandardSpillCuts);
93 
94  struct Selection
95  {
97  Cut cut;
98  };
99 
100  const int kNumCuts = 10;
101 
102  Selection vCut[kNumCuts] =
103  {
104  {"nocut", (!kIsDytmanMEC) && kNoCut},
105  // {"nocut", 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 
167  for(int j = 0; j < kNumCuts; ++j){
168 
169  Selection s = vCut[j];
170 
171  for(int i = 0; i < kNumPlots; ++i){
172 
173  Plot p = plots[i];
174 
175  hNonswap[i][j] = new Spectrum(p.label, p.bins, loaderNonswap, p.var, s.cut, kNoShift, kWOscDumb);
176  hSwap[i][j] = new Spectrum(p.label, p.bins, loaderSwap, p.var, s.cut, kNoShift, kWOscDumb);
177 
178  }
179  }
180 
181  // GO!
182  std::cout << "\nrun : --- run loaders.\n";
183  loaderNonswap.Go();
184  loaderSwap.Go();
185  std::cout << "\nrun : --- done.\n";
186 
187  std::cout << "\nrun : --- save output.\n";
188  TFile* f = new TFile(kOutputFileName.c_str(),"RECREATE");
189 
190  for(int j = 0; j < kNumCuts; ++j){
191  for(int i = 0; i < kNumPlots; ++i){
192 
193  std::string fullName = "reco-"+plots[i].fname+"-"+vCut[j].fname;
194 
195  TH1* NonSwap = hNonswap[i][j]->ToTH1(pot);
196  TH1* FluxSwap = hSwap[i][j]->ToTH1(pot);
197 
198  NonSwap->SetName(fullName.c_str());
199  NonSwap->Add(FluxSwap);
200 
201  NonSwap->Write();
202 
203  } // End loop over i
204  } // End loop over j
205 
206  TH1D* TotalPOT = new TH1D("meta-TotalPOT",";;Total POT", 1, 0, 1);
207  TotalPOT->SetBinContent(1, pot);
208  TotalPOT->Write();
209 
210  f->Write();
211  f->Close();
212 
213  std::cout << "\nrun : --- done.\n";
214 
215 } // End of function
216 
217 #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 Var kNumuHadCalE
Definition: NumuVars.cxx:538
const char * label
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
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 Var kCCE
Definition: NumuVars.h:21
const Binning kEnergyBinning
const double pot
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
const double j
Definition: BetheBloch.cxx:29
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 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 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
const Binning kZBins
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.
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));})
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
void caf_numu_fd_validation_MC_no_tau(std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kOutputFileName)
const Cut kNumuQuality
Definition: NumuCuts.h:18
Var var
Definition: CutFlow_Data.C:32
loaderSwap
Definition: demo4.py:9
SREnergyBranch energy
Energy estimator branch.
const Binning kXYBins
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