DataMCNDLoad_nus17.C
Go to the documentation of this file.
1 // A macro for ND data/MC comparisons
2 // Gavin S. Davies
3 
4 #include "CAFAna/Core/HistAxis.h"
5 #include "CAFAna/Core/Loaders.h"
6 #include "CAFAna/Core/Spectrum.h"
12 #include "CAFAna/Vars/HistAxes.h"
14 #include "NuXAna/Vars/NusVars.h"
15 #include "CAFAna/Vars/Vars.h"
16 
17 
18 using namespace ana;
19 
20 std::pair<Spectrum*, CheatDecomp*> make_pair(
21  SpectrumLoaderBase& loader_data,
22  SpectrumLoaderBase& loader_mc,
23  HistAxis* axis,
24  Cut* cut,
25  const SystShifts& shift,
26  const Var& wei
27 );
28 
29 std::map<std::string, std::pair<Spectrum*, CheatDecomp*> > MakeMap(
30  SpectrumLoaderBase& loader_data,
31  SpectrumLoaderBase& loader_mc,
32  std::map<std::string, HistAxis* > axes,
33  std::map<std::string, Cut* > cuts,
34  const SystShifts& shift,
35  const Var& wei
36 );
37 
38 std::map<std::string, std::pair<Spectrum*, CheatDecomp*> > MakeMap(
39  Loaders& loader,
40  std::map<std::string, HistAxis*> axes,
41  std::map<std::string, Cut*> cuts,
42  const SystShifts& shift,
43  const Var& wei
44 );
45 
47 {
48  TH1::AddDirectory(0);
49 
50  // Strings to file paths
51  const std::string fnearcaffull = "prod_decaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
52 
53  const std::string fneardatafull = "prod_decaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_nue_or_numu_or_nus_contain_v1_goodruns";
54 
55  // Set up a Loaders object
56  // For full SA data and MC filesets, use 'fneardatafull' and 'fnamenearcaf' respectively
60 
61  const Var kFluxMECWgt = kPPFXFluxCVWgt*kXSecCVWgt2017;
62 
63  const Var kHpP = SIMPLEVAR(sel.nuecosrej.hitsperplane);
64 
65  const Var kRecoVtxX(
66  [](const caf::SRProxy* sr)
67  {
68  if( !sr->vtx.elastic.IsValid ) { return 9000.f; }
69  return sr->vtx.elastic.vtx.X();
70  });
71  const Var kRecoVtxY(
72  [](const caf::SRProxy* sr)
73  {
74  if( !sr->vtx.elastic.IsValid) { return 9000.f; }
75  return sr->vtx.elastic.vtx.Y();
76  });
77  const Var kRecoVtxZ(
78  [](const caf::SRProxy* sr)
79  {
80  if( !sr->vtx.elastic.IsValid) { return 9000.f; }
81  return sr->vtx.elastic.vtx.Z();
82  });
83 
84  const Cut kNus17NDPtpCut(
85  [](const caf::SRProxy* sr)
86  {
87  double partptp = sr->sel.nuecosrej.partptp;
88  if(partptp > 0.8) return false;
89  return true;
90  });
91 
92  // Axis definitions
93  std::string labelHpP = "Hits per Plane";
94  std::string labelVtxX = "X Vertex Position (cm)";
95  std::string labelVtxY = "Y Vertex Position (cm)";
96  std::string labelVtxZ = "Z Vertex Position (cm)";
97  std::string labelTop = "Minimum Distance to Top Detector Face (cm)";
98  std::string labelBottom = "Minimum Distance to Bottom Detector Face (cm)";
99  std::string labelEast = "Minimum Distance to East Detector Face (cm)";
100  std::string labelWest = "Minimum Distance to West Detector Face (cm)";
101  std::string labelFront = "Minimum Distance to Front Detector Face (cm)";
102  std::string labelBack = "Minimum Distance to Back Detector Face (cm)";
103  std::string labelHits = "Number of Hits";
104  std::string labelCVN = "CVN NC Classifier";
105  std::string labelPTP = "p_{T}/p";
106  std::string labelRecoE = "Visible Energy (GeV)";
107 
108  Binning kNusEBins = Binning::Simple(100, 0, 10);
109  Binning kNCAllEBins = Binning::Simple(200, 0, 20);
110  Binning kXYBinsND = Binning::Simple(100, -250., 250.);
111  Binning kZBinsND = Binning::Simple(410, -50., 2000.);
112  Binning kNearEdgeBins = Binning::Simple(160, 0., 800.);
113  Binning kFinePIDBins = Binning::Simple(100, 0., 1.);
114  Binning kHitBins = Binning::Simple(800, 0., 800.);
115  Binning kHpPBins = Binning::Simple(10, 0., 10.);
116 
117  const HistAxis kNus17EAxis (labelRecoE, kNusEBins, kNus17Energy);
118  const HistAxis kNCAllEAxis (labelRecoE, kNCAllEBins, kNus17Energy);
119  const HistAxis kAxisHpP (labelHpP, kHpPBins, kHpP);
120  const HistAxis kAxisVtxX (labelVtxX, kXYBinsND, kRecoVtxX);
121  const HistAxis kAxisVtxY (labelVtxY, kXYBinsND, kRecoVtxY);
122  const HistAxis kAxisVtxZ (labelVtxZ, kZBinsND, kRecoVtxZ);
123  const HistAxis kAxisTop (labelTop, kNearEdgeBins, kDistAllTop);
124  const HistAxis kAxisBottom (labelBottom, kNearEdgeBins, kDistAllBottom);
125  const HistAxis kAxisEast (labelEast, kNearEdgeBins, kDistAllEast);
126  const HistAxis kAxisWest (labelWest, kNearEdgeBins, kDistAllWest);
127  const HistAxis kAxisFront (labelFront, kNearEdgeBins, kDistAllFront);
128  const HistAxis kAxisBack (labelBack, kNearEdgeBins, kDistAllBack);
129  const HistAxis kAxisCVN (labelCVN, kFinePIDBins, kCVNnc);
130  const HistAxis kAxisNHit (labelHits, kHitBins, kNHit);
131  const HistAxis kAxisPTP (labelPTP, kFinePIDBins, kPartPtp);
132 
133  // Two intermediate cuts, created from cuts in NusCuts.h
134  const Cut kNusEQFid = kNus17EventQuality && kNus17NDFiducial;
135 
136  const Cut kNM1Fid =
138  kNus17CVNnc &&
139  kNus17EnergyCut &&
140  kNus17NDPtpCut;
141 
142  const Cut kNM1Cont =
143  kNusEQFid &&
144  kNus17CVNnc &&
145  kNus17EnergyCut &&
146  kNus17NDPtpCut;
147 
148 
149  const Cut kNM1CVNnc =
150  kNus17NDPresel &&
151  kNus17EnergyCut &&
152  kNus17NDPtpCut;
153 
154  const Cut kNM1Ptp =
155  kNus17NDPresel &&
156  kNus17CVNnc &&
158 
159  const Cut kNM1Energy =
160  kNus17NDPresel &&
161  kNus17CVNnc &&
162  kNus17NDPtpCut;
163 
164 
165  const Cut kFidMinusX(
166  [](const caf::SRProxy* sr)
167  {
168  if( !sr->vtx.elastic.IsValid ) return false;
169 
170  const TVector3 vtx = sr->vtx.elastic.vtx;
171  if(vtx.Y() < -100.0) return false;
172  if(vtx.Y() > 100.0) return false;
173  if(vtx.Z() < 150.0) return false;
174  if(vtx.Z() > 1000.0) return false;
175  return true;
176  });
177  const Cut kFidMinusY(
178  [](const caf::SRProxy* sr)
179  {
180  if(!sr->vtx.elastic.IsValid ) return false;
181 
182  const TVector3 vtx = sr->vtx.elastic.vtx;
183  if(vtx.X() < -100.0) return false;
184  if(vtx.X() > 100.0) return false;
185  if(vtx.Z() < 150.0) return false;
186  if(vtx.Z() > 1000.0) return false;
187  return true;
188  });
189  const Cut kFidMinusZ(
190  [](const caf::SRProxy* sr)
191  {
192  if( !sr->vtx.elastic.IsValid ) return false;
193 
194  const TVector3 vtx = sr->vtx.elastic.vtx;
195  if(vtx.X() < -100.0) return false;
196  if(vtx.X() > 100.0) return false;
197  if(vtx.Y() < -100.0) return false;
198  if(vtx.Y() > 100.0) return false;
199  return true;
200  });
201 
202  const Cut kMinusTop =
203  kDistAllBottom > 25 &&
204  kDistAllEast > 25 &&
205  kDistAllWest > 25 &&
206  kDistAllFront > 25 &&
207  kDistAllBack > 25;
208 
209  const Cut kMinusBottom =
210  kDistAllTop > 25 &&
211  kDistAllEast > 25 &&
212  kDistAllWest > 25 &&
213  kDistAllFront > 25 &&
214  kDistAllBack > 25;
215 
216  const Cut kMinusEast =
217  kDistAllTop > 25 &&
218  kDistAllBottom > 25 &&
219  kDistAllWest > 25 &&
220  kDistAllFront > 25 &&
221  kDistAllBack > 25;
222 
223  const Cut kMinusWest =
224  kDistAllTop > 25 &&
225  kDistAllBottom > 25 &&
226  kDistAllEast > 25 &&
227  kDistAllFront > 25 &&
228  kDistAllBack > 25;
229 
230  const Cut kMinusFront =
231  kDistAllTop > 25 &&
232  kDistAllBottom > 25 &&
233  kDistAllEast > 25 &&
234  kDistAllWest > 25 &&
235  kDistAllBack > 25;
236 
237  const Cut kMinusBack =
238  kDistAllTop > 25 &&
239  kDistAllBottom > 25 &&
240  kDistAllEast > 25 &&
241  kDistAllWest > 25 &&
242  kDistAllFront > 25;
243 
244  // set up the axes map
245  std::map<std::string, HistAxis*> axes;
246  axes["RecoE"] = new HistAxis(kNus17EAxis);
247  axes["RecoELong"] = new HistAxis(kNCAllEAxis);
248  axes["HitpPlane"] = new HistAxis(kAxisHpP);
249  axes["VtxX"] = new HistAxis(kAxisVtxX);
250  axes["VtxY"] = new HistAxis(kAxisVtxY);
251  axes["VtxZ"] = new HistAxis(kAxisVtxZ);
252  axes["Top"] = new HistAxis(kAxisTop);
253  axes["Bottom"] = new HistAxis(kAxisBottom);
254  axes["East"] = new HistAxis(kAxisEast);
255  axes["West"] = new HistAxis(kAxisWest);
256  axes["Front"] = new HistAxis(kAxisFront);
257  axes["Back"] = new HistAxis(kAxisBack);
258  axes["CVN"] = new HistAxis(kAxisCVN);
259  axes["NHit"] = new HistAxis(kAxisNHit);
260  axes["PTP"] = new HistAxis(kAxisPTP);
261 
262  // set up the cuts map
263  std::map<std::string, Cut*> cuts;
264  cuts["DQ"] = new Cut(kNoCut);
265  cuts["EQ"] = new Cut(kNus17EventQuality);
266  cuts["EQFid"] = new Cut(kNusEQFid);
267  cuts["Presel"] = new Cut(kNus17NDPresel);
268  cuts["PreselCVN"] = new Cut(kNus17NDPresel && kNus17CVNnc);
269  cuts["PreselEnergy"] = new Cut(kNus17NDPresel && kNus17CVNnc && kNus17EnergyCut);
270  cuts["PreselPtp"] = new Cut(kNus17NDPresel && kNus17CVNnc && kNus17EnergyCut && kNus17NDPtpCut);
271 
272 
273  // make a map containing data spectra and MC cheat decomps for each axis/cut pair
274  std::map<std::string, std::pair<ana::Spectrum*, ana::CheatDecomp*> > loaded_specs =
275  MakeMap(loaders, axes, cuts, kNoShift, kFluxMECWgt);
276 
277  std::map<std::string, Cut*> nm1s;
278  nm1s["VtxX"] = new Cut(kNM1Fid && kFidMinusX);
279  nm1s["VtxY"] = new Cut(kNM1Fid && kFidMinusY);
280  nm1s["VtxZ"] = new Cut(kNM1Fid && kFidMinusZ);
281  nm1s["Top"] = new Cut(kNM1Cont && kMinusTop);
282  nm1s["Bottom"] = new Cut(kNM1Cont && kMinusBottom);
283  nm1s["East"] = new Cut(kNM1Cont && kMinusEast);
284  nm1s["West"] = new Cut(kNM1Cont && kMinusWest);
285  nm1s["Front"] = new Cut(kNM1Cont && kMinusFront);
286  nm1s["Back"] = new Cut(kNM1Cont && kMinusBack);
287  nm1s["CVN"] = new Cut(kNM1CVNnc);
288  nm1s["PTP"] = new Cut(kNM1Ptp);
289  nm1s["RecoE"] = new Cut(kNM1Energy);
290  nm1s["RecoELong"] = new Cut(kNM1Energy);
291 
292  for(const auto& nm1 : nm1s) {
293  std::string strcut = nm1.first;
294  std::string strNM1 = "NM1";
295 
296  loaded_specs[strcut + strNM1] =
297  make_pair(
300  axes[strcut], nm1.second, kNoShift, kFluxMECWgt
301  );
302  }
303 
304  // Go loaders
306  loaders.Go();
307 
308  // Open output file
309  TFile* rootOut = new TFile(outfile.c_str(), "RECREATE");
310 
311  for(const auto& loaded_spec : loaded_specs) {
312  loaded_spec.second.first->SaveTo(rootOut, (loaded_spec.first + "Data").c_str());
313  loaded_spec.second.second->SaveTo(rootOut, (loaded_spec.first + "MC").c_str());
314  }
315 
316  rootOut->Close();
317 }
318 
319 //..............................................................................
320 std::pair<Spectrum*, CheatDecomp*> make_pair(SpectrumLoaderBase& loader_data,
321  SpectrumLoaderBase& loader_mc,
322  HistAxis* axis,
323  Cut* cut,
324  const SystShifts& shift,
325  const Var& wei)
326 {
327  std::pair<Spectrum*, CheatDecomp*> output;
328  output.first = new Spectrum(
329  loader_data,
330  *axis, *cut,
331  shift, wei
332  );
333  output.second = new CheatDecomp(
334  loader_mc,
335  *axis, *cut,
336  shift, wei
337  );
338 
339  return output;
340 }
341 
342 //..............................................................................
343 std::map<std::string, std::pair<Spectrum*, CheatDecomp*> > MakeMap(
344  SpectrumLoaderBase& loader_data,
345  SpectrumLoaderBase& loader_mc,
346  std::map<std::string, HistAxis* > axes,
347  std::map<std::string, Cut* > cuts,
348  const SystShifts& shift,
349  const Var& wei
350 )
351 {
352  std::map<std::string, std::pair<Spectrum*, CheatDecomp*> > output;
353  for(const auto& cut : cuts) {
354  for(const auto& axis : axes) {
355  output[axis.first + cut.first] = make_pair(
356  loader_data, loader_mc,
357  axis.second, cut.second, shift, wei
358  );
359  }
360  }
361 
362  return output;
363 };
364 
365 //..............................................................................
366 std::map<std::string, std::pair<Spectrum*, CheatDecomp*> > MakeMap(
367  Loaders& loader,
368  std::map<std::string, HistAxis*> axes,
369  std::map<std::string, Cut*> cuts,
370  const SystShifts& shift,
371  const Var& wei
372 )
373 {
374  return MakeMap(
377  axes, cuts, shift, wei
378  );
379 };
Near Detector underground.
Definition: SREnums.h:10
ofstream output
const Cut kNus17EnergyCut([](const caf::SRProxy *sr){double energy=0.0;energy=kNus17Energy(sr);if(energy< 0.25) return false;if(energy > 10.0) return false;return true;})
Definition: NusCuts17.h:51
const Var kPartPtp([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.partptp)|| std::isinf(1.*sr->sel.nuecosrej.partptp)) return-5.f;return float(sr->sel.nuecosrej.partptp);})
Definition: NusVars.h:74
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
float partptp
Definition: NusVarsTemp.cxx:50
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kDistAllBottom([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngbottom)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngbottom);})
Distance of all showers in slice from the bottom edge of detector.
Definition: NueVars.h:33
const Var kCVNnc
PID
Definition: Vars.cxx:44
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var kDistAllWest([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngwest)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngwest);})
Distance of all showers in slice from the west edge of detector.
Definition: NueVars.h:36
const Var kDistAllTop([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngtop)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngtop);})
Distance of all showers in slice from the top edge of detector.
Definition: NueVars.h:30
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Var kDistAllBack([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngback)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngback);})
Distance of all showers in slice from the back edge of detector.
Definition: NueVars.h:42
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
const Cut kNus17NDFiducial([](const caf::SRProxy *sr){assert(sr->vtx.elastic.IsValid &&"Must apply EQ cuts");if(sr->vtx.elastic.vtx.X()< -100.0) return false;if(sr->vtx.elastic.vtx.X() > 100.0) return false;if(sr->vtx.elastic.vtx.Y()< -100.0) return false;if(sr->vtx.elastic.vtx.Y() > 100.0) return false;if(sr->vtx.elastic.vtx.Z()< 150.0) return false;if(sr->vtx.elastic.vtx.Z() > 1000.0) return false;return true;})
Definition: NusCuts17.h:103
const Cut kNus17CVNnc([](const caf::SRProxy *sr){if(sr->slc.nhit < 25) return false;std::cout<< "ERROR::kNus17CVNnc. Looking for cvnProd3Train. Branch no longer exists."<< std::endl;abort();return true;})
Definition: NusCuts17.h:49
const Binning kZBinsND
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
const Var kNus17Energy([](const caf::SRProxy *sr){double cale=sr->slc.calE;double recoE=0.;if(sr->hdr.det==caf::kFARDET) recoE=FDscaleCalE17 *cale;if(sr->hdr.det==caf::kNEARDET) recoE=NDscaleCalE17 *cale;return recoE;})
Definition: NusVars.h:62
const Var kDistAllEast([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngeast)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngeast);})
Distance of all showers in slice from the east edge of detector.
Definition: NueVars.h:39
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
caf::Proxy< caf::SRNueCosRej > nuecosrej
Definition: SRProxy.h:1265
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
const Cut kNus17NDPresel
Definition: NusCuts17.h:114
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
Definition: Loaders.cxx:129
const Var kNHit
Definition: Vars.cxx:71
const ana::Var kRecoVtxY
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
void DataMCNDLoad_nus17(std::string outfile)
caf::StandardRecord * sr
caf::Proxy< float > partptp
Definition: SRProxy.h:1057
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:610
std::map< std::string, std::pair< Spectrum *, CheatDecomp * > > MakeMap(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, std::map< std::string, HistAxis * > axes, std::map< std::string, Cut * > cuts, const SystShifts &shift, const Var &wei)
const SystShifts kNoShift
Definition: SystShifts.cxx:21
const Cut kNus17NDContain
Definition: NusCuts17.h:105
Base class for the various types of spectrum loader.
const Cut cut
Definition: exporter_fd.C:30
const Binning kXYBinsND
const Var kDistAllFront([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngfront)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngfront);})
Distance of all showers in slice from the front edge of detector.
Definition: NueVars.h:45
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
std::vector< Loaders * > loaders
Definition: syst_header.h:386
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2073
const Binning kFinePIDBins
Definition: NusLoadProd3.h:102
const ana::Var kRecoVtxZ
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
Definition: Loaders.cxx:25
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Cut kNus17EventQuality([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->sel.nuecosrej.hitsperplane >=8) return false;if(sr->vtx.elastic.fuzzyk.npng==0) return false;if(sr->slc.ncontplanes<=2) return false;return true;})
Data Quality cuts from docdb 21113.
Definition: NusCuts17.h:21
Just return the ND truth spectra as the decomposition.
Definition: CheatDecomp.h:10
const ana::Var kRecoVtxX
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
FILE * outfile
Definition: dump_event.C:13
enum BeamMode string