CVNProngVars.cxx
Go to the documentation of this file.
2 
3 #include "CAFAna/Vars/Vars.h"
4 
6 
7 #include <cassert>
8 
9 namespace ana
10 {
11 
12  const Var kCVNMuonIdx([](const caf::SRProxy *sr)
13  {
14  float longest_idx = -5.0;
15  float longest_len = -5.0;
16 
17  // loop over all verticies, prongs, and tracks to find the longest muon track
18  if(sr->vtx.elastic.IsValid){
19 
20  for(unsigned int png_idx = 0; png_idx < sr->vtx.elastic.fuzzyk.png.size(); png_idx++) {
21 
22  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
23 
24  if(png.len > longest_len) {
25  longest_len = png.len;
26  longest_idx = png_idx;
27  }
28 
29  } // end loop over prongs
30  } // end loop over verticies
31 
32 
33 
34  // If there are ANY prongs longer than 500 cm, use the longest prong:
35  if(longest_len > 500.0) return longest_idx;
36 
37 
38 
39  // Otherwise find the highest scoring prong...
40  float best_idx = -5.0;
41  float best_score = -5.0;
42  if(sr->vtx.elastic.IsValid){
43 
44  for(unsigned int png_idx = 0; png_idx < sr->vtx.elastic.fuzzyk.png.size(); png_idx++) {
45 
46  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
47 
48  if(png.cvnpart.muonid > best_score) {
49  best_score = png.cvnpart.muonid;
50  best_idx = png_idx;
51  }
52 
53  } // end loop over prongs
54  } // end loop over verticies
55 
56  return best_idx;
57 
58  });
59 
60 
61 
62  const Var kCVNBestMuonScore([](const caf::SRProxy *sr)
63  {
64  float muonScore = -5.0;
65 
66  if(kCVNMuonIdx(sr) < 0.0) return muonScore;
67  if(!sr->vtx.elastic.IsValid) return muonScore;
68 
69  // HARDCODED ASSUMPTION: about only having one vertex... ...for now...
70  muonScore = sr->vtx.elastic.fuzzyk.png[(unsigned int)kCVNMuonIdx(sr)].cvnpart.muonid;
71 
72  // If the prong was longer than 500.0 cm, assign it a score of 1.0.
73  if(sr->vtx.elastic.fuzzyk.png[(unsigned int)kCVNMuonIdx(sr)].len > 500.0) muonScore = 1.0;
74 
75  return muonScore;
76 
77  });
78 
79  const Var kCVNBestMuonRawScore([](const caf::SRProxy *sr)
80  {
81  float muonScore = -5.0;
82 
83  if(kCVNMuonIdx(sr) < 0.0) return muonScore;
84  if(!sr->vtx.elastic.IsValid) return muonScore;
85 
86  // HARDCODED ASSUMPTION: about only having one vertex... ...for now...
87  muonScore = sr->vtx.elastic.fuzzyk.png[(unsigned int)kCVNMuonIdx(sr)].cvnpart.muonid;
88 
89  return muonScore;
90  });
91 
92 
93  /// Function to return "adjusted" CVN-prong muon score:
94  double GetCVNProngMuonScore(unsigned int ProngIdx, caf::SRProxy *sr) {
95 
96  // Sanity checks first to prevent segfaults...
97  if(!sr->vtx.elastic.IsValid) return -1.0;
98  if(sr->vtx.elastic.fuzzyk.npng == 0) return -1.0;
99  if(ProngIdx >= sr->vtx.elastic.fuzzyk.npng) return -1.0;
100 
101  double score = sr->vtx.elastic.fuzzyk.png[ProngIdx].cvnpart.muonid;
102  double len = sr->vtx.elastic.fuzzyk.png[ProngIdx].len;
103 
104  if(len < 500.0) return score;
105  else return 1.0;
106 
107  } // end GetCVNProngMuonScore
108 
109 
110 
111 
112 
113  // ==================== //
114  //
115  // Functions to return "combo PID" scores.
116  //
117  // ==================== //
118 
119  // EM score for this prong from summing all "shower like" PIDs.
120  double EMScore(unsigned int ProngIdx, const caf::SRProxy* sr) {
121 
122  // Sanity checks first to prevent segfaults...
123  if(!sr->vtx.elastic.IsValid) return -1.0;
124  if(sr->vtx.elastic.fuzzyk.npng == 0) return -1.0;
125  if(ProngIdx >= sr->vtx.elastic.fuzzyk.npng) return -1.0;
126 
127  double val = (sr->vtx.elastic.fuzzyk.png[ProngIdx].cvnpart.electronid +
128  sr->vtx.elastic.fuzzyk.png[ProngIdx].cvnpart.photonid +
129  sr->vtx.elastic.fuzzyk.png[ProngIdx].cvnpart.pizeroid);
130 
131  return val;
132 
133  } // end EMScore
134 
135  // Had score for this prong from summing all "hadronic like" PIDs.
136  //
137  // NOTE: This function literaly just returns the sum of the CVN-prong
138  // hadron scores. Some people think this is a bad function to use since
139  // CVN-prong can sometimes confuse pions and muons. If you think the muon
140  // score should be included, then use the HadScoreWithMuon function.
141  double HadScore(unsigned int ProngIdx, const caf::SRProxy* sr) {
142 
143  // Sanity checks first to prevent segfaults...
144  if(!sr->vtx.elastic.IsValid) return -1.0;
145  if(sr->vtx.elastic.fuzzyk.npng == 0) return -1.0;
146  if(ProngIdx >= sr->vtx.elastic.fuzzyk.npng) return -1.0;
147 
148  double val = (sr->vtx.elastic.fuzzyk.png[ProngIdx].cvnpart.protonid +
149  sr->vtx.elastic.fuzzyk.png[ProngIdx].cvnpart.pionid +
150  sr->vtx.elastic.fuzzyk.png[ProngIdx].cvnpart.neutronid);
151 
152  return val;
153 
154  } // end HadScore
155 
156  // BIG NOTE:
157  //
158  // For this, we are just taking the had score to be 1-EMScore. This will include the
159  // muon score in the Had score. This is based on the assumption that you have already
160  // selected the lepton prong and that pions might look a little bit like muons. Argue
161  // with me on slack if you disagree...
162  double HadScoreWithMuon(unsigned int ProngIdx, const caf::SRProxy* sr) {
163 
164  // Sanity checks first to prevent segfaults...
165  if(!sr->vtx.elastic.IsValid) return -1.0;
166  if(sr->vtx.elastic.fuzzyk.npng == 0) return -1.0;
167  if(ProngIdx >= sr->vtx.elastic.fuzzyk.npng) return -1.0;
168 
169  // just do the simple thing for now and take 1 - EMScore
170  double val = 1.0 - EMScore(ProngIdx,sr);
171 
172  return val;
173 
174  } // end HadScoreWithMuon
175 
176 
177 
178 
179 
180  // ==================== //
181  //
182  // Helper functions and variables for choosing prongs based on CVN score:
183  //
184  // ==================== //
185 
186 
187  /// IsNumuCCPionByCVN
188  /// \brief: Returns a bool if prong passes NumuCC event pion selection criteria
189  ///
190  /// Ignoring the prong selected to be a muon, return true if this prong passes pion selection.
191  ///
192  /// NOTE: Some assumptions are hardcoded in here:
193  /// 1) There is only one vertex.
194  /// 2) We are only asking this question of a 3D prong.
195  bool IsNumuCCPionByCVN(unsigned int ProngIdx, const caf::SRProxy* sr) {
196 
197  // Sanity checks first to prevent segfaults...
198  if(!sr->vtx.elastic.IsValid) return false;
199  if(sr->vtx.elastic.fuzzyk.npng == 0) return false;
200  if(ProngIdx >= sr->vtx.elastic.fuzzyk.npng) return false;
201 
202  // First check if this is the selected muon...
203  if(ProngIdx == (unsigned int)kCVNMuonIdx(sr)) return false;
204 
205  // Check the CVN muon and pion scores. If either score is high enough (because pions sort of look like muons) take it.
206  // To guarantee a prong is not selected more than once, we must also require that each id it selects on is also the maximum value for all prongs.
207  double muonid = sr->vtx.elastic.fuzzyk.png[ProngIdx].cvnpart.muonid;
208  double pionid = sr->vtx.elastic.fuzzyk.png[ProngIdx].cvnpart.pionid;
209  double maxval = sr->vtx.elastic.fuzzyk.png[ProngIdx].cvnpart.maxval;
210  if( (muonid > 0.7 && muonid == maxval) || (pionid > 0.7 && pionid == maxval) ) return true;
211 
212  // Otherwise, this is not a pion...
213  return false;
214 
215  } // end IsNumuCCPionByCVN
216 
217 
218 
219  /// IsNumuCCProtonByCVN
220  /// \brief: Returns a bool if prong passes NumuCC event proton selection criteria
221  ///
222  /// Ignoring the prong selected to be a muon, return true if this prong passes proton selection.
223  ///
224  /// NOTE: Some assumptions are hardcoded in here:
225  /// 1) There is only one vertex.
226  /// 2) We are only asking this question of a 3D prong.
227  bool IsNumuCCProtonByCVN(unsigned int ProngIdx, const caf::SRProxy* sr) {
228 
229  // Sanity checks first to prevent segfaults...
230  if(!sr->vtx.elastic.IsValid) return false;
231  if(sr->vtx.elastic.fuzzyk.npng == 0) return false;
232  if(ProngIdx >= sr->vtx.elastic.fuzzyk.npng) return false;
233 
234  // First check if this is the selected muon...
235  if(ProngIdx == (unsigned int)kCVNMuonIdx(sr)) return false;
236 
237  // Check the CVN proton score...
238  // To guarantee a prong is not selected more than once, we must also require that each id it selects on is also the maximum value for all prongs.
239  double protonid = sr->vtx.elastic.fuzzyk.png[ProngIdx].cvnpart.protonid;
240  double maxval = sr->vtx.elastic.fuzzyk.png[ProngIdx].cvnpart.maxval;
241  if( protonid > 0.7 && protonid == maxval ) return true;
242 
243  // Otherwise, this is not a proton...
244  return false;
245 
246  } // end IsNumuCCProtonByCVN
247 
248 
249  /// IsNumuCCHadronShowerByCVN
250  /// \brief: Returns a bool if prong passes NumuCC event hadron selection criteria
251  ///
252  /// Ignoring the prong selected to be a muon, return true if this prong passes hadron selection
253  ///
254  /// NOTE: Some assumptions are hardcoded in here:
255  /// 1) There is only one vertex.
256  /// 2) We are only asking this question of a 3D prong.
257  bool IsNumuCCHadronByCVN(unsigned int ProngIdx, const caf::SRProxy* sr) {
258 
259  // Sanity checks first to prevent segfaults...
260  if(!sr->vtx.elastic.IsValid) return false;
261  if(sr->vtx.elastic.fuzzyk.npng == 0) return false;
262  if(ProngIdx >= sr->vtx.elastic.fuzzyk.npng) return false;
263 
264  // First check if this is the selected muon...
265  if(ProngIdx == (unsigned int)kCVNMuonIdx(sr)) return false;
266 
267  if ( HadScoreWithMuon(ProngIdx,sr) >= 0.5 ) return true;
268 
269  // Otherwise, not a hadron
270  return false;
271  } // end IsNumuCCHadronByCVN
272 
273 
274  /// IsNumuCCEMShowerByCVN
275  /// \brief: Returns a bool if prong passes NumuCC event EM selection criteria
276  ///
277  /// Ignoring the prong selected to be a muon, return true if this prong passes EM shower selection.
278  ///
279  /// NOTE: Some assumptions are hardcoded in here:
280  /// 1) There is only one vertex.
281  /// 2) We are only asking this question of a 3D prong.
282  bool IsNumuCCEMShowerByCVN(unsigned int ProngIdx, const caf::SRProxy* sr) {
283 
284  // Sanity checks first to prevent segfaults...
285  if(!sr->vtx.elastic.IsValid) return false;
286  if(sr->vtx.elastic.fuzzyk.npng == 0) return false;
287  if(ProngIdx >= sr->vtx.elastic.fuzzyk.npng) return false;
288 
289  // First check if this is the selected muon...
290  if(ProngIdx == (unsigned int)kCVNMuonIdx(sr)) return false;
291 
292  if( EMScore(ProngIdx,sr) > 0.5 ) return true;
293 
294  // Otherwise, this is not an EM shower...
295  return false;
296 
297  } // end IsNumuCCEMShowerByCVN
298 
299 
300  /// kNumuCCPionIndices
301  /// \brief: Returns a vector of indices for prongs selected as pions
302  ///
303  const MultiVar kNumuCCPionIndices([](const caf::SRProxy *sr) {
304 
305  std::vector<double> indices;
306  if(sr->vtx.elastic.IsValid){
307  for(unsigned int png_idx = 0; png_idx < sr->vtx.elastic.fuzzyk.png.size(); png_idx++) {
308 
309  // Call helper function defined above to make a decision about this prong.
310  if( IsNumuCCPionByCVN(png_idx,sr) ) indices.push_back((double)png_idx);
311 
312  } // end loop over prongs
313  } // end loop over verticies
314 
315  return indices;
316 
317  }); // end kNumuCCPionIndices
318 
319 
320  const MultiVar kNumuCCProtonIndices([](const caf::SRProxy *sr) {
321 
322  std::vector<double> indices;
323  if(sr->vtx.elastic.IsValid){
324  for(unsigned int png_idx = 0; png_idx < sr->vtx.elastic.fuzzyk.png.size(); png_idx++) {
325 
326  // Call helper function defined above to make a decision about this prong.
327  if( IsNumuCCProtonByCVN(png_idx,sr) ) indices.push_back((double)png_idx);
328 
329  } // end loop over prongs
330  } // end loop over verticies
331 
332  return indices;
333 
334  }); // end kNumuCCProtonIndices
335 
336 
337  const MultiVar kNumuCCEMShowerIndices([](const caf::SRProxy *sr) {
338 
339  std::vector<double> indices;
340 
341  if(sr->vtx.elastic.IsValid){
342  for(unsigned int png_idx = 0; png_idx < sr->vtx.elastic.fuzzyk.png.size(); png_idx++) {
343 
344  // Call helper function defined above to make a decision about this prong.
345  if( IsNumuCCEMShowerByCVN(png_idx,sr) ) indices.push_back((double)png_idx);
346 
347  } // end loop over prongs
348  } // end loop over verticies
349 
350  return indices;
351 
352  }); // end kNumuCCEMShowerIndices
353 
354 }
const XML_Char int len
Definition: expat.h:262
caf::Proxy< size_t > npng
Definition: SRProxy.h:2038
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double EMScore(unsigned int ProngIdx, const caf::SRProxy *sr)
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const MultiVar kNumuCCPionIndices([](const caf::SRProxy *sr){std::vector< double > indices;if(sr->vtx.elastic.IsValid){for(unsigned int png_idx=0;png_idx< sr->vtx.elastic.fuzzyk.png.size();png_idx++){if(IsNumuCCPionByCVN(png_idx, sr)) indices.push_back((double) png_idx);}}return indices;})
: Returns a vector of indices for prongs selected as pions
Definition: CVNProngVars.h:69
double HadScoreWithMuon(unsigned int ProngIdx, const caf::SRProxy *sr)
const Var kCVNMuonIdx([](const caf::SRProxy *sr){float longest_idx=-5.0;float longest_len=-5.0;if(sr->vtx.elastic.IsValid){for(unsigned int png_idx=0;png_idx< sr->vtx.elastic.fuzzyk.png.size();png_idx++){auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];if(png.len > longest_len){longest_len=png.len;longest_idx=png_idx;}}} if(longest_len > 500.0) return longest_idx;float best_idx=-5.0;float best_score=-5.0;if(sr->vtx.elastic.IsValid){for(unsigned int png_idx=0;png_idx< sr->vtx.elastic.fuzzyk.png.size();png_idx++){auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];if(png.cvnpart.muonid > best_score){best_score=png.cvnpart.muonid;best_idx=png_idx;}}}return best_idx;})
: Prong index of best muon prong by CVN score & length
Definition: CVNProngVars.h:19
A PID for muons.
Definition: FillPIDs.h:11
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
bool IsNumuCCPionByCVN(unsigned int ProngIdx, const caf::SRProxy *sr)
: Returns a bool if prong passes NumuCC event pion selection criteria
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2043
const MultiVar kNumuCCEMShowerIndices([](const caf::SRProxy *sr){std::vector< double > indices;if(sr->vtx.elastic.IsValid){for(unsigned int png_idx=0;png_idx< sr->vtx.elastic.fuzzyk.png.size();png_idx++){if(IsNumuCCEMShowerByCVN(png_idx, sr)) indices.push_back((double) png_idx);}}return indices;})
: Returns a vector of indices for prongs selected as EM showers
Definition: CVNProngVars.h:79
bool IsNumuCCEMShowerByCVN(unsigned int ProngIdx, const caf::SRProxy *sr)
: Returns a bool if prong passes NumuCC event EM selection criteria
const MultiVar kNumuCCProtonIndices([](const caf::SRProxy *sr){std::vector< double > indices;if(sr->vtx.elastic.IsValid){for(unsigned int png_idx=0;png_idx< sr->vtx.elastic.fuzzyk.png.size();png_idx++){if(IsNumuCCProtonByCVN(png_idx, sr)) indices.push_back((double) png_idx);}}return indices;})
: Returns a vector of indices for prongs selected as protons
Definition: CVNProngVars.h:74
caf::StandardRecord * sr
double HadScore(unsigned int ProngIdx, const caf::SRProxy *sr)
double GetCVNProngMuonScore(unsigned int ProngIdx, caf::SRProxy *sr)
Function to return "adjusted" CVN-prong muon score:
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
const Var kCVNBestMuonRawScore([](const caf::SRProxy *sr){float muonScore=-5.0;if(kCVNMuonIdx(sr)< 0.0) return muonScore;if(!sr->vtx.elastic.IsValid) return muonScore;muonScore=sr->vtx.elastic.fuzzyk.png[(unsigned int) kCVNMuonIdx(sr)].cvnpart.muonid;return muonScore;})
: Raw muon score for best muon prong by CVN score & length
Definition: CVNProngVars.h:43
bool IsNumuCCHadronByCVN(unsigned int ProngIdx, const caf::SRProxy *sr)
: Returns a bool if prong passes NumuCC event hadron selection criteria
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
const Var kCVNBestMuonScore([](const caf::SRProxy *sr){float muonScore=-5.0;if(kCVNMuonIdx(sr)< 0.0) return muonScore;if(!sr->vtx.elastic.IsValid) return muonScore;muonScore=sr->vtx.elastic.fuzzyk.png[(unsigned int) kCVNMuonIdx(sr)].cvnpart.muonid;if(sr->vtx.elastic.fuzzyk.png[(unsigned int) kCVNMuonIdx(sr)].len > 500.0) muonScore=1.0;return muonScore;})
: Muon score for best muon prong by CVN score & length
Definition: CVNProngVars.h:27
bool IsNumuCCProtonByCVN(unsigned int ProngIdx, const caf::SRProxy *sr)
: Returns a bool if prong passes NumuCC event proton selection criteria