DrawPreSelectionPlots.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/Vars/Vars.h"
33 
38 
39 #include "NuXAna/Vars/NusVars.h"
40 #include "NuXAna/Cuts/NusCuts20.h"
41 #include "NuXAna/Cuts/NusCuts.h"
42 
43 // selection includes
45 
46 // root includes
47 #include "TFile.h"
48 #include "TCanvas.h"
49 #include "TH1.h"
50 #include "TH2.h"
51 
52 using namespace ana;
53 
54 const Cut kCCNC([](const caf::SRProxy* sr) {
55  if (sr->mc.nnu == 0) return false;
56  return (sr->mc.nu[0].iscc) ? false : true;
57  });
58 
60 
61 void DrawPreSelectionPlots(const std::string simName = "",
62  const std::string datName = ""){
63 
64  // setup output file
65  std::string fileOutName = "PreSelectionPlots.root";
66  TFile fileOut(fileOutName.c_str(), "RECREATE");
67 
68  std::vector< ana::Cut > cutLevels = {
69  kPreSelection && kNoCut,
70  kPreSelection && kNoCut && kTruthSelection,
71  kPreSelection && kNoCut && kNus20NDTrueVtxOutOfDet,
72  kPreSelection && kNus20NDQualityCuts,
73  kPreSelection && kNus20NDQualityCuts && kTruthSelection,
74  kPreSelection && kNus20NDQualityCuts && kNus20NDTrueVtxOutOfDet,
75  kPreSelection && kNus20NDQualityCuts && kNus20NDVtxPosCut,
76  kPreSelection && kNus20NDQualityCuts && kNus20NDVtxPosCut && kTruthSelection,
77  kPreSelection && kNus20NDQualityCuts && kNus20NDVtxPosCut && kNus20NDTrueVtxOutOfDet,
78  kPreSelection && kNus20NDQualityCuts && kNus20NDFiducialCuts,
79  kPreSelection && kNus20NDQualityCuts && kNus20NDFiducialCuts && kTruthSelection,
80  kPreSelection && kNus20NDQualityCuts && kNus20NDFiducialCuts &&kNus20NDTrueVtxOutOfDet,
81  kPreSelection && kNus20NDCuts,
82  kPreSelection && kNus20NDCuts && kTruthSelection,
83  kPreSelection && kNus20NDCuts && kNus20NDTrueVtxOutOfDet
84  };
85 
86  std::vector< std::string > cutNames = {
87  "NoCut" ,
88  "NoCutTrueNC",
89  "NoCutTrueVtxOutOfDet",
90  "Quality",
91  "QualityTrueNC",
92  "QualityTrueVtxOutOfDet",
93  "QualityVtx",
94  "QualityVtxTrueNC",
95  "QualityVtxTrueVtxOutOfDet",
96  "QualityVtxFiducial",
97  "QualityVtxFiducialTrueNC",
98  "QualityVtxFiducialTrueVtxOutOfDet",
99  "FullSelection",
100  "FullSelectionTrueNC",
101  "FullSelectionTrueVtxOutOfDet"
102  };
103 
104  std::vector<std::string> loaderNames = {
105  simName,
106  datName
107  };
108 
109  std::vector<std::string> loaderShortName = {
110  "sim",
111  "dat"
112  };
113 
114 
115  for (size_t iLoader = 0; iLoader < loaderNames.size(); iLoader++){
116 
117  for (size_t iCut = 0; iCut < cutLevels.size(); iCut++){
118 
119  if (loaderShortName.at(iLoader) == "dat" &&
120  (cutNames.at(iCut).find("True") != std::string::npos)){
121  std::cout
122  << "loaderShortName: "
123  << loaderShortName.at(iLoader)
124  << " and "
125  << "cutName: "
126  << cutNames.at(iCut)
127  << std::endl;
128  continue;
129  }
130 
131  if (loaderNames.at(iLoader) == "") continue;
132 
133  // Spectrum loader must be inside the loop, since
134  // thisLoader.Go() has to be run after the Spectrum
135  // have been defined
136  SpectrumLoader thisLoader(loaderNames.at(iLoader).c_str());
137 
138  // phase space cuts
139  Spectrum kHitsPerPlaneSpect("kHitsPerPlaneSpect",
141  thisLoader,
143  cutLevels.at(iCut));
144 
145 
146  Spectrum kNSliceHitsSpect("kNSliceHitsSpect",
147  kNSliceHitsBins,
148  thisLoader,
149  kNSliceHits,
150  cutLevels.at(iCut));
151 
152  Spectrum kPartPtpSpect("kPartPtpSpect",
153  kPartPtpBins,
154  thisLoader,
155  kPartPtp,
156  cutLevels.at(iCut));
157 
158  Spectrum kCVNncSpect("kCVNncSpect",
159  kCVNBins,
160  thisLoader,
161  kCVNnc,
162  cutLevels.at(iCut));
163 
164  Spectrum kCVNnc_looseptpSpect("kCVNnc_looseptpSpect",
165  kCVNBins,
166  thisLoader,
168  cutLevels.at(iCut));
169 
170  Spectrum kCVNnc_oldpreselSpect("kCVNnc_oldpreselSpect",
171  kCVNBins,
172  thisLoader,
174  cutLevels.at(iCut));
175 
176  Spectrum kCaloESpect("kCaloESpect",
177  kNus20EnergyBinning,
178  thisLoader,
179  kCaloE,
180  cutLevels.at(iCut));
181 
182  // for TH2s
183  Spectrum kHitsPerPlaneNSliceHitsSpect("kHitsPerPlaneNSliceHitsSpect",
184  thisLoader,
186  kNSliceHitsBins , kNSliceHits,
187  cutLevels.at(iCut));
188 
189  Spectrum kHitsPerPlanePartPtpSpect("kHitsPerPlanePartPtpSpect",
190  thisLoader,
193  cutLevels.at(iCut));
194 
195  Spectrum kNSliceHitsPartPtpSpect("kNSliceHitsPartPtpSpect",
196  thisLoader,
197  kNSliceHitsBins, kNSliceHits,
199  cutLevels.at(iCut));
200 
201  Spectrum kHitsPerPlaneCaloESpect("kHitsPerPlaneCaloESpect",
202  thisLoader,
203  kNus20EnergyBinning, kCaloE,
205  cutLevels.at(iCut));
206 
207  Spectrum kNSliceHitsCaloESpect("kNSliceHitsCaloESpect",
208  thisLoader,
209  kNus20EnergyBinning, kCaloE,
210  kNSliceHitsBins, kNSliceHits,
211  cutLevels.at(iCut));
212 
213  Spectrum kPartPtpCaloESpect("kPartPtpCaloESpect",
214  thisLoader,
215  kNus20EnergyBinning, kCaloE,
217  cutLevels.at(iCut));
218 
219  Spectrum kCVNnc_looseptpCaloESpect("kCVNnc_looseptpCaloESpec",
220  thisLoader,
221  kNus20EnergyBinning, kCaloE,
223  cutLevels.at(iCut));
224 
225  Spectrum kCVNnc_oldpreselCaloESpect("kCVNnc_oldpreselCaloESpect",
226  thisLoader,
227  kNus20EnergyBinning, kCaloE,
229  cutLevels.at(iCut));
230 
231  // quality cuts
232  Spectrum kIsElasticSpect("kIsElasticSpect",
234  thisLoader,
235  kIsElastic,
236  cutLevels.at(iCut));
237 
238  Spectrum kNFuzzyProngSpect("kNFuzzyProngSpect",
240  thisLoader,
241  kNFuzzyProng,
242  cutLevels.at(iCut));
243 
244  Spectrum kNContPlanesSpect("kNContPlanesSpect",
246  thisLoader,
247  kNContPlanes,
248  cutLevels.at(iCut));
249 
250  // containment cuts
251  Spectrum kVtxXSpect("kVtxXSpect",
252  kVtxXBins,
253  thisLoader,
254  kVtxX,
255  cutLevels.at(iCut));
256 
257  Spectrum kVtxYSpect("kVtxYSpect",
258  kVtxYBins,
259  thisLoader,
260  kVtxY,
261  cutLevels.at(iCut));
262 
263  Spectrum kVtxZSpect("kVtxZSpect",
264  kVtxZBins,
265  thisLoader,
266  kVtxZ,
267  cutLevels.at(iCut));
268 
269 
270  Spectrum kDistAllTopSpect("kDistAllTopSpect",
271  kDistAllTopBins,
272  thisLoader,
273  kDistAllTop,
274  cutLevels.at(iCut));
275 
276  Spectrum kDistAllBottomSpect("kDistAllBottomSpect",
277  kDistAllBottomBins,
278  thisLoader,
280  cutLevels.at(iCut));
281 
282  Spectrum kDistAllEastSpect("kDistAllEastSpect",
283  kDistAllEastBins,
284  thisLoader,
285  kDistAllEast,
286  cutLevels.at(iCut));
287 
288  Spectrum kDistAllWestSpect("kDistAllWestSpect",
289  kDistAllWestBins,
290  thisLoader,
291  kDistAllWest,
292  cutLevels.at(iCut));
293 
294  Spectrum kDistAllFrontSpect("kDistAllFrontSpect",
295  kDistAllFrontBins,
296  thisLoader,
298  cutLevels.at(iCut));
299 
300  Spectrum kDistAllBackSpect("kDistAllBackSpect",
301  kDistAllBackBins,
302  thisLoader,
303  kDistAllBack,
304  cutLevels.at(iCut));
305 
306 
307  thisLoader.Go();
308 
309  // now get histograms
310  float thisPOT = kIsElasticSpect.POT();
311  TH1* kHitsPerPlaneHist = kHitsPerPlaneSpect .ToTH1(thisPOT);
312  TH1* kNSliceHitsHist = kNSliceHitsSpect .ToTH1(thisPOT);
313  TH1* kPartPtpHist = kPartPtpSpect .ToTH1(thisPOT);
314  TH1* kCVNncHist = kCVNncSpect .ToTH1(thisPOT);
315  TH1* kCVNnc_looseptpHist = kCVNnc_looseptpSpect .ToTH1(thisPOT);
316  TH1* kCVNnc_oldpreselHist = kCVNnc_oldpreselSpect .ToTH1(thisPOT);
317  TH1* kIsElasticHist = kIsElasticSpect .ToTH1(thisPOT);
318  TH1* kNFuzzyProngHist = kNFuzzyProngSpect .ToTH1(thisPOT);
319  TH1* kNContPlanesHist = kNContPlanesSpect .ToTH1(thisPOT);
320  TH1* kVtxXHist = kVtxXSpect .ToTH1(thisPOT);
321  TH1* kVtxYHist = kVtxYSpect .ToTH1(thisPOT);
322  TH1* kVtxZHist = kVtxZSpect .ToTH1(thisPOT);
323  TH1* kDistAllTopHist = kDistAllTopSpect .ToTH1(thisPOT);
324  TH1* kDistAllBottomHist = kDistAllBottomSpect .ToTH1(thisPOT);
325  TH1* kDistAllEastHist = kDistAllEastSpect .ToTH1(thisPOT);
326  TH1* kDistAllWestHist = kDistAllWestSpect .ToTH1(thisPOT);
327  TH1* kDistAllFrontHist = kDistAllFrontSpect .ToTH1(thisPOT);
328  TH1* kDistAllBackHist = kDistAllBackSpect .ToTH1(thisPOT);
329  TH1* kCaloEHist = kCaloESpect .ToTH1(thisPOT);
330  TH2* kHitsPerPlaneNSliceHitsHist = kHitsPerPlaneNSliceHitsSpect.ToTH2(thisPOT);
331  TH2* kHitsPerPlanePartPtpHist = kHitsPerPlanePartPtpSpect .ToTH2(thisPOT);
332  TH2* kNSliceHitsPartPtpHist = kNSliceHitsPartPtpSpect .ToTH2(thisPOT);
333  TH2* kHitsPerPlaneCaloEHist = kHitsPerPlaneCaloESpect .ToTH2(thisPOT);
334  TH2* kNSliceHitsCaloEHist = kNSliceHitsCaloESpect .ToTH2(thisPOT);
335  TH2* kPartPtpCaloEHist = kPartPtpCaloESpect .ToTH2(thisPOT);
336  TH2* kCVNnc_looseptpCaloEHist = kCVNnc_looseptpCaloESpect .ToTH2(thisPOT);
337  TH2* kCVNnc_oldpreselCaloEHist = kCVNnc_oldpreselCaloESpect .ToTH2(thisPOT);
338 
339  fileOut.cd();
340 
341  kHitsPerPlaneHist -> Write((loaderShortName.at(iLoader) + "kHitsPerPlaneHist" + cutNames.at(iCut)).c_str());
342  kNSliceHitsHist -> Write((loaderShortName.at(iLoader) + "kNSliceHitsHist" + cutNames.at(iCut)).c_str());
343  kPartPtpHist -> Write((loaderShortName.at(iLoader) + "kPartPtpHist" + cutNames.at(iCut)).c_str());
344  if (loaderShortName.at(iLoader) != "dat" ||
345  (loaderShortName.at(iLoader)=="dat" && cutNames.at(iCut) == "FullSelection" )){
346  kCaloEHist -> Write((loaderShortName.at(iLoader) + "kCaloEHist" + cutNames.at(iCut)).c_str());
347  }
348  if (loaderShortName.at(iLoader) != "dat" ||
349  (loaderShortName.at(iLoader)=="dat" && cutNames.at(iCut) == "FullSelection" ) ||
350  (loaderShortName.at(iLoader)=="dat" && cutNames.at(iCut) == "QualityVtxFiducial")){
351  kCVNncHist -> Write((loaderShortName.at(iLoader) + "kCVNncHist" + cutNames.at(iCut)).c_str());
352  kCVNnc_looseptpHist -> Write((loaderShortName.at(iLoader) + "kCVNnc_looseptpHist" + cutNames.at(iCut)).c_str());
353  kCVNnc_oldpreselHist -> Write((loaderShortName.at(iLoader) + "kCVNnc_oldpreselHist" + cutNames.at(iCut)).c_str());
354  }
355  kIsElasticHist -> Write((loaderShortName.at(iLoader) + "kIsElasticHist" + cutNames.at(iCut)).c_str());
356  kNFuzzyProngHist -> Write((loaderShortName.at(iLoader) + "kNFuzzyProngHist" + cutNames.at(iCut)).c_str());
357  kNContPlanesHist -> Write((loaderShortName.at(iLoader) + "kNContPlanesHist" + cutNames.at(iCut)).c_str());
358  kVtxXHist -> Write((loaderShortName.at(iLoader) + "kVtxXHist" + cutNames.at(iCut)).c_str());
359  kVtxYHist -> Write((loaderShortName.at(iLoader) + "kVtxYHist" + cutNames.at(iCut)).c_str());
360  kVtxZHist -> Write((loaderShortName.at(iLoader) + "kVtxZHist" + cutNames.at(iCut)).c_str());
361  kDistAllTopHist -> Write((loaderShortName.at(iLoader) + "kDistAllTopHist" + cutNames.at(iCut)).c_str());
362  kDistAllBottomHist -> Write((loaderShortName.at(iLoader) + "kDistAllBottomHist" + cutNames.at(iCut)).c_str());
363  kDistAllEastHist -> Write((loaderShortName.at(iLoader) + "kDistAllEastHist" + cutNames.at(iCut)).c_str());
364  kDistAllWestHist -> Write((loaderShortName.at(iLoader) + "kDistAllWestHist" + cutNames.at(iCut)).c_str());
365  kDistAllFrontHist -> Write((loaderShortName.at(iLoader) + "kDistAllFrontHist" + cutNames.at(iCut)).c_str());
366  kDistAllBackHist -> Write((loaderShortName.at(iLoader) + "kDistAllBackHist" + cutNames.at(iCut)).c_str());
367  kHitsPerPlaneNSliceHitsHist -> Write((loaderShortName.at(iLoader)+"kHitsPerPlaneNSliceHitsHist" + cutNames.at(iCut)).c_str());
368  kHitsPerPlanePartPtpHist -> Write((loaderShortName.at(iLoader)+"kHitsPerPlanePartPtpHist" + cutNames.at(iCut)).c_str());
369  kNSliceHitsPartPtpHist -> Write((loaderShortName.at(iLoader)+"kNSliceHitsPartPtpHist" + cutNames.at(iCut)).c_str());
370  kHitsPerPlaneCaloEHist -> Write((loaderShortName.at(iLoader)+"kHitsPerPlaneCaloEHist" + cutNames.at(iCut)).c_str());
371  kNSliceHitsCaloEHist -> Write((loaderShortName.at(iLoader)+"kNSliceHitsCaloEHist" + cutNames.at(iCut)).c_str());
372  kPartPtpCaloEHist -> Write((loaderShortName.at(iLoader)+"kPartPtpCaloEHist" + cutNames.at(iCut)).c_str());
373  kCVNnc_looseptpCaloEHist -> Write((loaderShortName.at(iLoader)+"kCVNnc_looseptpCaloEHist" + cutNames.at(iCut)).c_str());
374  kCVNnc_oldpreselCaloEHist -> Write((loaderShortName.at(iLoader)+"kCVNnc_oldpreselCaloEHist" + cutNames.at(iCut)).c_str());
375 
376  fileOut.Write();
377  }
378  }
379 }
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 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 Cut kCCNC([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;return(sr->mc.nu[0].iscc)?false:true;})
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 Cut kTruthSelection
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
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
const Var kVtxZ
void DrawPreSelectionPlots(const std::string simName="", const std::string datName="")
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
gm Write()
enum BeamMode string