MakeSelectionPlots.C
Go to the documentation of this file.
1 // script to plot the variables for the NuS20 preselection,
2 // along with the correlation between each variable and
3 // each other variable.
4 //
5 // list of preselection cuts and variables:
6 // kIsElastic --- is there an elastic vertex? 2 , 0, 2
7 // kHitsPerPlane --- number of hits per plane 100, 0, 20
8 // kNFuzzyProng --- number of fuzzyk prongs 10, 0, 10
9 // kVtxX/kVtxY/kVtxZ --- vertex positions 100, -300, 300 . 100, -300, 300 . 100, -50, 1600
10 // kNSliceHits --- number of hits in the slice 200, 0, 200
11 // kNContPlanes --- number of continuouse planes 120, 0, 120
12 // kDistAllTop 100, 0, 600
13 // kDistAllBottom 100, 0, 600
14 // kDistAllEast 100, 0, 600
15 // kDistAllWest 100, 0, 600
16 // kDistAllFront 100, 0, 1600
17 // kDistAllBack 100, 0, 1600
18 // kPartPtp 100, 0, 1
19 //
20 // author: Adam Lister
21 // date : 2020-02-19
22 // email : adam.lister@wisc.edu
23 
24 // framework includes
26 #include "CAFAna/Core/Binning.h"
27 #include "CAFAna/Core/Spectrum.h"
29 #include "CAFAna/Core/EventList.h"
30 #include "CAFAna/Core/Utilities.h"
31 #include "CAFAna/Cuts/Cuts.h"
32 #include "CAFAna/Cuts/NusCuts20.h"
33 #include "CAFAna/Cuts/NusCuts.h"
34 #include "CAFAna/Vars/Vars.h"
35 #include "CAFAna/Vars/NusVars.h"
36 #include "CAFAna/Vars/NueVarsExtra.h"
37 #include "CAFAna/Vars/NumuVarsExtra.h"
38 #include "CAFAna/Cuts/NueCuts2020.h"
39 #include "CAFAna/Cuts/NumuCuts2020.h"
40 
41 // selection includes
42 #include "CAFAna/nus/Nus20/Selection/BinningDefinitions.h"
43 
44 // root includes
45 #include "TFile.h"
46 #include "TCanvas.h"
47 #include "TH1.h"
48 #include "TH2.h"
49 
50 using namespace ana;
51 
52 const Cut kCCNC([](const caf::SRProxy* sr) {
53  if (sr->mc.nnu == 0) return false;
54  return (sr->mc.nu[0].iscc) ? false : true;
55  });
56 
58 
59 void MakeSelectionPlots(const std::string simName = "",
60  const std::string datName = ""){
61 
62  // setup output file
63  std::string fileOutName = "PreSelectionPlots.root";
64  TFile fileOut(fileOutName.c_str(), "RECREATE");
65 
66  std::vector< ana::Cut > cutLevels = {
70  kNus20NDPreSelection && kNus20NDQualityCuts,
71  kNus20NDPreSelection && kNus20NDQualityCuts && kTruthSelection,
72  kNus20NDPreSelection && kNus20NDQualityCuts && kNus20NDTrueVtxOutOfDet,
73  kNus20NDPreSelection && kNus20NDQualityCuts && kNus20NDVtxPosCut,
74  kNus20NDPreSelection && kNus20NDQualityCuts && kNus20NDVtxPosCut && kTruthSelection,
75  kNus20NDPreSelection && kNus20NDQualityCuts && kNus20NDVtxPosCut && kNus20NDTrueVtxOutOfDet,
76  kNus20NDPreSelection && kNus20NDQualityCuts && kNus20NDFiducialCuts,
77  kNus20NDPreSelection && kNus20NDQualityCuts && kNus20NDFiducialCuts && kTruthSelection,
78  kNus20NDPreSelection && kNus20NDQualityCuts && kNus20NDFiducialCuts &&kNus20NDTrueVtxOutOfDet,
80  kNus20NDPreSelection && kNus20NDCuts && kTruthSelection,
81  kNus20NDPreSelection && kNus20NDCuts && kNus20NDTrueVtxOutOfDet
82  };
83 
84  std::vector< std::string > cutNames = {
85  "NoCut" ,
86  "NoCutTrueNC",
87  "NoCutTrueVtxOutOfDet",
88  "Quality",
89  "QualityTrueNC",
90  "QualityTrueVtxOutOfDet",
91  "QualityVtx",
92  "QualityVtxTrueNC",
93  "QualityVtxTrueVtxOutOfDet",
94  "QualityVtxFiducial",
95  "QualityVtxFiducialTrueNC",
96  "QualityVtxFiducialTrueVtxOutOfDet",
97  "FullSelection",
98  "FullSelectionTrueNC",
99  "FullSelectionTrueVtxOutOfDet"
100  };
101 
102  std::vector<std::string> loaderNames = {
103  simName,
104  datName
105  };
106 
107  std::vector<std::string> loaderShortName = {
108  "sim",
109  "dat"
110  };
111 
112 
113  for (size_t iLoader = 0; iLoader < loaderNames.size(); iLoader++){
114 
115  for (size_t iCut = 0; iCut < cutLevels.size(); iCut++){
116 
117  if (loaderShortName.at(iLoader) == "dat" &&
118  (cutNames.at(iCut).find("True") != std::string::npos)){
119  std::cout
120  << "loaderShortName: "
121  << loaderShortName.at(iLoader)
122  << " and "
123  << "cutName: "
124  << cutNames.at(iCut)
125  << std::endl;
126  continue;
127  }
128 
129  if (loaderNames.at(iLoader) == "") continue;
130 
131  // Spectrum loader must be inside the loop, since
132  // thisLoader.Go() has to be run after the Spectrum
133  // have been defined
134  SpectrumLoader thisLoader(loaderNames.at(iLoader).c_str());
135 
136  // phase space cuts
137  Spectrum kHitsPerPlaneSpect("kHitsPerPlaneSpect",
139  thisLoader,
141  cutLevels.at(iCut));
142 
143 
144  Spectrum kNSliceHitsSpect("kNSliceHitsSpect",
145  kNSliceHitsBins,
146  thisLoader,
147  kNSliceHits,
148  cutLevels.at(iCut));
149 
150  Spectrum kPartPtpSpect("kPartPtpSpect",
151  kPartPtpBins,
152  thisLoader,
153  kPartPtp,
154  cutLevels.at(iCut));
155 
156  Spectrum kCVNncSpect("kCVNncSpect",
157  kCVNBins,
158  thisLoader,
159  kCVNnc,
160  cutLevels.at(iCut));
161 
162  Spectrum kCVNnc_looseptpSpect("kCVNnc_looseptpSpect",
163  kCVNBins,
164  thisLoader,
166  cutLevels.at(iCut));
167 
168  Spectrum kCVNnc_oldpreselSpect("kCVNnc_oldpreselSpect",
169  kCVNBins,
170  thisLoader,
172  cutLevels.at(iCut));
173 
174  Spectrum kCaloESpect("kCaloESpect",
175  kNus20EnergyBinning,
176  thisLoader,
177  kCaloE,
178  cutLevels.at(iCut));
179 
180  // for TH2s
181  Spectrum kHitsPerPlaneNSliceHitsSpect("kHitsPerPlaneNSliceHitsSpect",
182  thisLoader,
184  kNSliceHitsBins , kNSliceHits,
185  cutLevels.at(iCut));
186 
187  Spectrum kHitsPerPlanePartPtpSpect("kHitsPerPlanePartPtpSpect",
188  thisLoader,
191  cutLevels.at(iCut));
192 
193  Spectrum kNSliceHitsPartPtpSpect("kNSliceHitsPartPtpSpect",
194  thisLoader,
195  kNSliceHitsBins, kNSliceHits,
197  cutLevels.at(iCut));
198 
199  Spectrum kHitsPerPlaneCaloESpect("kHitsPerPlaneCaloESpect",
200  thisLoader,
201  kNus20EnergyBinning, kCaloE,
203  cutLevels.at(iCut));
204 
205  Spectrum kNSliceHitsCaloESpect("kNSliceHitsCaloESpect",
206  thisLoader,
207  kNus20EnergyBinning, kCaloE,
208  kNSliceHitsBins, kNSliceHits,
209  cutLevels.at(iCut));
210 
211  Spectrum kPartPtpCaloESpect("kPartPtpCaloESpect",
212  thisLoader,
213  kNus20EnergyBinning, kCaloE,
215  cutLevels.at(iCut));
216 
217  Spectrum kCVNnc_looseptpCaloESpect("kCVNnc_looseptpCaloESpec",
218  thisLoader,
219  kNus20EnergyBinning, kCaloE,
221  cutLevels.at(iCut));
222 
223  Spectrum kCVNnc_oldpreselCaloESpect("kCVNnc_oldpreselCaloESpect",
224  thisLoader,
225  kNus20EnergyBinning, kCaloE,
227  cutLevels.at(iCut));
228 
229  // quality cuts
230  Spectrum kIsElasticSpect("kIsElasticSpect",
232  thisLoader,
233  kIsElastic,
234  cutLevels.at(iCut));
235 
236  Spectrum kNFuzzyProngSpect("kNFuzzyProngSpect",
238  thisLoader,
239  kNFuzzyProng,
240  cutLevels.at(iCut));
241 
242  Spectrum kNContPlanesSpect("kNContPlanesSpect",
244  thisLoader,
245  kNContPlanes,
246  cutLevels.at(iCut));
247 
248  // containment cuts
249  Spectrum kVtxXSpect("kVtxXSpect",
250  kVtxXBins,
251  thisLoader,
252  kVtxX,
253  cutLevels.at(iCut));
254 
255  Spectrum kVtxYSpect("kVtxYSpect",
256  kVtxYBins,
257  thisLoader,
258  kVtxY,
259  cutLevels.at(iCut));
260 
261  Spectrum kVtxZSpect("kVtxZSpect",
262  kVtxZBins,
263  thisLoader,
264  kVtxZ,
265  cutLevels.at(iCut));
266 
267 
268  Spectrum kDistAllTopSpect("kDistAllTopSpect",
269  kDistAllTopBins,
270  thisLoader,
271  kDistAllTop,
272  cutLevels.at(iCut));
273 
274  Spectrum kDistAllBottomSpect("kDistAllBottomSpect",
275  kDistAllBottomBins,
276  thisLoader,
278  cutLevels.at(iCut));
279 
280  Spectrum kDistAllEastSpect("kDistAllEastSpect",
281  kDistAllEastBins,
282  thisLoader,
283  kDistAllEast,
284  cutLevels.at(iCut));
285 
286  Spectrum kDistAllWestSpect("kDistAllWestSpect",
287  kDistAllWestBins,
288  thisLoader,
289  kDistAllWest,
290  cutLevels.at(iCut));
291 
292  Spectrum kDistAllFrontSpect("kDistAllFrontSpect",
293  kDistAllFrontBins,
294  thisLoader,
296  cutLevels.at(iCut));
297 
298  Spectrum kDistAllBackSpect("kDistAllBackSpect",
299  kDistAllBackBins,
300  thisLoader,
301  kDistAllBack,
302  cutLevels.at(iCut));
303 
304 
305  thisLoader.Go();
306 
307  // now get histograms
308  float thisPOT = spectra.at(0) .POT();
309  TH1* kHitsPerPlaneHist = spectra.at(0) .ToTH1(thisPOT);
310  TH1* kNSliceHitsHist = kNSliceHitsSpect .ToTH1(thisPOT);
311  TH1* kPartPtpHist = kPartPtpSpect .ToTH1(thisPOT);
312  TH1* kCVNncHist = kCVNncSpect .ToTH1(thisPOT);
313  TH1* kCVNnc_looseptpHist = kCVNnc_looseptpSpect .ToTH1(thisPOT);
314  TH1* kCVNnc_oldpreselHist = kCVNnc_oldpreselSpect .ToTH1(thisPOT);
315  TH1* kIsElasticHist = kIsElasticSpect .ToTH1(thisPOT);
316  TH1* kNFuzzyProngHist = kNFuzzyProngSpect .ToTH1(thisPOT);
317  TH1* kNContPlanesHist = kNContPlanesSpect .ToTH1(thisPOT);
318  TH1* kVtxXHist = kVtxXSpect .ToTH1(thisPOT);
319  TH1* kVtxYHist = kVtxYSpect .ToTH1(thisPOT);
320  TH1* kVtxZHist = kVtxZSpect .ToTH1(thisPOT);
321  TH1* kDistAllTopHist = kDistAllTopSpect .ToTH1(thisPOT);
322  TH1* kDistAllBottomHist = kDistAllBottomSpect .ToTH1(thisPOT);
323  TH1* kDistAllEastHist = kDistAllEastSpect .ToTH1(thisPOT);
324  TH1* kDistAllWestHist = kDistAllWestSpect .ToTH1(thisPOT);
325  TH1* kDistAllFrontHist = kDistAllFrontSpect .ToTH1(thisPOT);
326  TH1* kDistAllBackHist = kDistAllBackSpect .ToTH1(thisPOT);
327  TH1* kCaloEHist = kCaloESpect .ToTH1(thisPOT);
328  TH2* kHitsPerPlaneNSliceHitsHist = kHitsPerPlaneNSliceHitsSpect.ToTH2(thisPOT);
329  TH2* kHitsPerPlanePartPtpHist = kHitsPerPlanePartPtpSpect .ToTH2(thisPOT);
330  TH2* kNSliceHitsPartPtpHist = kNSliceHitsPartPtpSpect .ToTH2(thisPOT);
331  TH2* kHitsPerPlaneCaloEHist = kHitsPerPlaneCaloESpect .ToTH2(thisPOT);
332  TH2* kNSliceHitsCaloEHist = kNSliceHitsCaloESpect .ToTH2(thisPOT);
333  TH2* kPartPtpCaloEHist = kPartPtpCaloESpect .ToTH2(thisPOT);
334  TH2* kCVNnc_looseptpCaloEHist = kCVNnc_looseptpCaloESpect .ToTH2(thisPOT);
335  TH2* kCVNnc_oldpreselCaloEHist = kCVNnc_oldpreselCaloESpect .ToTH2(thisPOT);
336 
337  fileOut.cd();
338 
339  kHitsPerPlaneHist -> Write((loaderShortName.at(iLoader) + "kHitsPerPlaneHist" + cutNames.at(iCut)).c_str());
340  kNSliceHitsHist -> Write((loaderShortName.at(iLoader) + "kNSliceHitsHist" + cutNames.at(iCut)).c_str());
341  kPartPtpHist -> Write((loaderShortName.at(iLoader) + "kPartPtpHist" + cutNames.at(iCut)).c_str());
342  if (loaderShortName.at(iLoader) != "dat" ||
343  (loaderShortName.at(iLoader)=="dat" && cutNames.at(iCut) == "FullSelection" )){
344  kCaloEHist -> Write((loaderShortName.at(iLoader) + "kCaloEHist" + cutNames.at(iCut)).c_str());
345  }
346  if (loaderShortName.at(iLoader) != "dat" ||
347  (loaderShortName.at(iLoader)=="dat" && cutNames.at(iCut) == "FullSelection" ) ||
348  (loaderShortName.at(iLoader)=="dat" && cutNames.at(iCut) == "QualityVtxFiducial")){
349  kCVNncHist -> Write((loaderShortName.at(iLoader) + "kCVNncHist" + cutNames.at(iCut)).c_str());
350  kCVNnc_looseptpHist -> Write((loaderShortName.at(iLoader) + "kCVNnc_looseptpHist" + cutNames.at(iCut)).c_str());
351  kCVNnc_oldpreselHist -> Write((loaderShortName.at(iLoader) + "kCVNnc_oldpreselHist" + cutNames.at(iCut)).c_str());
352  }
353  kIsElasticHist -> Write((loaderShortName.at(iLoader) + "kIsElasticHist" + cutNames.at(iCut)).c_str());
354  kNFuzzyProngHist -> Write((loaderShortName.at(iLoader) + "kNFuzzyProngHist" + cutNames.at(iCut)).c_str());
355  kNContPlanesHist -> Write((loaderShortName.at(iLoader) + "kNContPlanesHist" + cutNames.at(iCut)).c_str());
356  kVtxXHist -> Write((loaderShortName.at(iLoader) + "kVtxXHist" + cutNames.at(iCut)).c_str());
357  kVtxYHist -> Write((loaderShortName.at(iLoader) + "kVtxYHist" + cutNames.at(iCut)).c_str());
358  kVtxZHist -> Write((loaderShortName.at(iLoader) + "kVtxZHist" + cutNames.at(iCut)).c_str());
359  kDistAllTopHist -> Write((loaderShortName.at(iLoader) + "kDistAllTopHist" + cutNames.at(iCut)).c_str());
360  kDistAllBottomHist -> Write((loaderShortName.at(iLoader) + "kDistAllBottomHist" + cutNames.at(iCut)).c_str());
361  kDistAllEastHist -> Write((loaderShortName.at(iLoader) + "kDistAllEastHist" + cutNames.at(iCut)).c_str());
362  kDistAllWestHist -> Write((loaderShortName.at(iLoader) + "kDistAllWestHist" + cutNames.at(iCut)).c_str());
363  kDistAllFrontHist -> Write((loaderShortName.at(iLoader) + "kDistAllFrontHist" + cutNames.at(iCut)).c_str());
364  kDistAllBackHist -> Write((loaderShortName.at(iLoader) + "kDistAllBackHist" + cutNames.at(iCut)).c_str());
365  kHitsPerPlaneNSliceHitsHist -> Write((loaderShortName.at(iLoader)+"kHitsPerPlaneNSliceHitsHist" + cutNames.at(iCut)).c_str());
366  kHitsPerPlanePartPtpHist -> Write((loaderShortName.at(iLoader)+"kHitsPerPlanePartPtpHist" + cutNames.at(iCut)).c_str());
367  kNSliceHitsPartPtpHist -> Write((loaderShortName.at(iLoader)+"kNSliceHitsPartPtpHist" + cutNames.at(iCut)).c_str());
368  kHitsPerPlaneCaloEHist -> Write((loaderShortName.at(iLoader)+"kHitsPerPlaneCaloEHist" + cutNames.at(iCut)).c_str());
369  kNSliceHitsCaloEHist -> Write((loaderShortName.at(iLoader)+"kNSliceHitsCaloEHist" + cutNames.at(iCut)).c_str());
370  kPartPtpCaloEHist -> Write((loaderShortName.at(iLoader)+"kPartPtpCaloEHist" + cutNames.at(iCut)).c_str());
371  kCVNnc_looseptpCaloEHist -> Write((loaderShortName.at(iLoader)+"kCVNnc_looseptpCaloEHist" + cutNames.at(iCut)).c_str());
372  kCVNnc_oldpreselCaloEHist -> Write((loaderShortName.at(iLoader)+"kCVNnc_oldpreselCaloEHist" + cutNames.at(iCut)).c_str());
373 
374  fileOut.Write();
375  }
376  }
377 }
const ana::Binning kPartPtpBins
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
const ana::Binning kNContPlanesBins
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 Cut kNus20NDPreSelection
Definition: NusCuts20.h:61
const Var kCVNnc
PID
Definition: Vars.cxx:44
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 kIsElastic
Definition: NusVars.cxx:11
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
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
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
const Var kNContPlanes
Definition: NusVars.cxx:13
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
const Var kVtxX
const Cut kNus20NDTrueVtxOutOfDet
Definition: NusCuts20.h:112
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Cut kNus20NDCuts
Definition: NusCuts20.h:102
const Var kNFuzzyProng
Definition: NusVars.cxx:12
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
const Var kHitsPerPlane
const ana::Binning kNFuzzyProngBins
const Cut kNus20NDFiducialCuts
Definition: NusCuts20.h:88
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
const Cut kCCNC([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;return(sr->mc.nu[0].iscc)?false:true;})
caf::StandardRecord * sr
const Var kNSliceHits
Definition: NusVars.cxx:18
const ana::Binning kHitsPerPlaneBins
double POT() const
Definition: Spectrum.h:227
const ana::Binning kCVNBins
OStream cout
Definition: OStream.cxx:6
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
const Var kCVNnc_oldpresel
Definition: NusVars.cxx:15
const ana::Binning kIsElasticBins
const Var kCVNnc_looseptp
Definition: NusVars.cxx:14
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
const Var kVtxY
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Cut kNus20NDVtxPosCut
Definition: NusCuts20.h:71
void MakeSelectionPlots(const std::string simName="", const std::string datName="")
const Var kVtxZ
const Cut kTruthSelection
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
gm Write()
enum BeamMode string