NCPi0SingleProngCuts.h
Go to the documentation of this file.
1 //This header is to define any cuts that may be necessary when performing the single prong neutral current pi0 analysis.
2 
3 #pragma once
4 #include "CAFAna/Core/Cut.h"
6 #include "CAFAna/Vars/Vars.h"
9 #include "CAFAna/Cuts/Cuts.h"
10 #include "TMath.h"
11 #include "TVector3.h"
12 #include "TParticlePDG.h"
13 
14 namespace ana
15 {
16  /****************************************************************************/
17  //Make the individual cuts to be used
18 
19 
20  //Create the fiducial region in which vertex hits must be located.
21  const Cut kFiducial([](const caf::SRProxy* sr)
22  {
23  if(sr->vtx.nelastic == 0 && sr ->mc.nnu == 0) return false;
24  if((sr->vtx.elastic[0].vtx.x) < -170.0) return false;
25  if((sr->vtx.elastic[0].vtx.x) > 170.0) return false;
26  if((sr->vtx.elastic[0].vtx.y) < -160.0) return false;
27  if((sr->vtx.elastic[0].vtx.y) > 170.0) return false;
28  if((sr->vtx.elastic[0].vtx.z) < 50.0) return false;
29  if((sr->vtx.elastic[0].vtx.z) > 1100.0) return false;
30  return true;
31  });
32 
33  //Make sure the true vertices are inside the detector
34  const Cut kTrueFiducial([](const caf::SRProxy* sr)
35  {
36  if (sr->mc.nnu == 0) return false;
37  if (sr->vtx.nelastic == 0) return false;
38  if(fabs(sr->mc.nu[0].vtx.x) > 200.0) return false;
39  if(fabs(sr->mc.nu[0].vtx.y) > 200.0) return false;
40  if((sr->mc.nu[0].vtx.z) < 0.0) return false;
41  if((sr->mc.nu[0].vtx.z) > 1600.0) return false;
42  return true;
43  });
44 
45 
46  //Create the containtment region in which shower hits must be located.
47  const Cut kContainment([](const caf::SRProxy* sr)
48  { //Check that there is a vertex, shower, and prongs.
49  if(sr->mc.nnu==0) return false;
50  if(sr->vtx.nelastic==0) return false;
51  if(sr->vtx.elastic[0].fuzzyk.nshwlid==0) return false;
52  if(sr->vtx.elastic[0].fuzzyk.npng==0) return false;
53  for (unsigned int i=0; i<sr->vtx.elastic[0].fuzzyk.nshwlid; i++)
54  { //Define the area of containment.
55  if ((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.x) < -100.0) return false;
56  if ((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.x) > 125.0) return false;
57  if ((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.y) < -120.0) return false;
58  if ((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.y) > 100.0) return false;
59  if ((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.z) < 300.0) return false;
60  if ((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.z) > 1200.0) return false;
61  }
62  return true;
63  });
64 
65 
66  //Choose the slices that have a single vertex in them.
67  const Cut kVtx([](const caf::SRProxy* sr)
68  {
69  if(sr->vtx.nelastic == 0 && sr->mc.nnu == 0)return false;
70  return (sr->vtx.nelastic==1);
71  });
72 
73  const Cut kVertexDistance([](const caf::SRProxy* sr)
74  {
75  if(sr->vtx.nelastic == 0)return false;
76  float Distance;
77  Distance = sqrt(
78  (sr->mc.nu[0].vtx.x - sr->vtx.elastic[0].vtx.x) *
79  (sr->mc.nu[0].vtx.x - sr->vtx.elastic[0].vtx.x) +
80  (sr->mc.nu[0].vtx.y - sr->vtx.elastic[0].vtx.y) *
81  (sr->mc.nu[0].vtx.y - sr->vtx.elastic[0].vtx.y) +
82  (sr->mc.nu[0].vtx.z - sr->vtx.elastic[0].vtx.z) *
83  (sr->mc.nu[0].vtx.z - sr->vtx.elastic[0].vtx.z)
84  );
85  return Distance > 10.0;
86  });
87 
88 
89 
90  //Select only single prong events.
91  const Cut k1Prong([](const caf::SRProxy* sr)
92  {
93  if(sr->vtx.nelastic==0 && sr->mc.nnu == 0) return false;
94  return (sr->vtx.elastic[0].fuzzyk.npng == 1 );
95  });
96 
97  //Identify particles as muon type neutrinos.
98  const Cut kNumu([](const caf::SRProxy* sr)
99  {
100  if (sr->mc.nnu == 0) return false;
101  if (sr->vtx.nelastic == 0) return false;
102  return(sr->mc.nu[0].iscc &&
103  sr->mc.nu[0].pdg == 14);
104  });
105 
106  //Identify particles as muon type neutrinos.
107  const Cut kNue([](const caf::SRProxy* sr)
108  {
109  if (sr->mc.nnu == 0) return false;
110  if (sr->vtx.nelastic == 0) return false;
111  return(sr->mc.nu[0].iscc &&
112  sr->mc.nu[0].pdg == 12);
113  });
114 
115 
116  //Pick out the neutral current interactions.
117  const Cut kNCurrent([](const caf::SRProxy* sr)
118  {
119  if (sr->mc.nnu == 0) return false;
120  if (sr->vtx.nelastic == 0) return false;
121  return (sr->mc.nu[0].iscc == 0);
122  });
123 
124  //Pick out the charged current interactions.
125  const Cut kCCurrent([](const caf::SRProxy* sr)
126  {
127  if (sr->mc.nnu == 0) return false;
128  if (sr->vtx.nelastic == 0) return false;
129  return (sr->mc.nu[0].iscc==1);
130  });
131 
132 
133  //Pick out events that have a neutral pion within them.
134  const Cut kPi0([](const caf::SRProxy* sr)
135  {
136  float MassOfPi0=0.135;
137  float kinetic=-10;
138  float en = -99;
139  if(sr->mc.nnu == 0) return false;
140  if(sr->vtx.nelastic == 0) return false;
141  if(sr->mc.nu[0].prim.size() == 0) return false;
142  int nbofprim=sr->mc.nu[0].prim.size();
143  int countpi=0;
144  for(int i = 0; i < nbofprim; i++)
145  {
146  if(sr->mc.nu[0].prim[i].pdg == 111)
147  {
148  en= sr->mc.nu[0].prim[i].p.E;
149  kinetic = en - MassOfPi0;
150  //if(kinetic > 0.5) countpi++;
151  if(kinetic > 0.0) countpi++;
152  }
153  }
154  if(countpi > 0) return true;
155  return false;
156  });
157 
158 
159  //Define neutral current background events.
160  const Cut NCBkg([](const caf::SRProxy* sr)
161  {
162  float MassOfPi0=0.135;
163  float kinetic=-10;
164  float en = -99;
165  if(sr->mc.nnu == 0) return false;
166  if(sr->vtx.nelastic == 0) return false;
167  if(sr->mc.nu[0].prim.size() == 0) return false;
168  int nbofprim=sr->mc.nu[0].prim.size();
169  int countpi=0;
170  for(int i = 0; i < nbofprim; i++)
171  {
172  if(sr->mc.nu[0].prim[i].pdg == 111)
173  { //May need to change kinetic energy.
174  en= sr->mc.nu[0].prim[i].p.E;
175  kinetic = en - MassOfPi0;
176  countpi++;
177  //if(kinetic > 0.0) countpi++;
178  }
179  }
180  if (countpi == 0 && !sr->mc.nu[0].iscc) return true;
181  return false;
182  });
183 
184 
185  //Define charged current pi0 background events.
186  const Cut CCPi0Bkg([](const caf::SRProxy* sr)
187  {
188  if(sr->mc.nnu == 0) return false;
189  if(sr->vtx.nelastic == 0) return false;
190  if(sr->mc.nu[0].prim.size() == 0) return false;
191  int nbofprim=sr->mc.nu[0].prim.size();
192  int countpi=0;
193  // bool ispi0=false;
194  for(int i=0;i<nbofprim;i++)
195  {
196  if(sr->mc.nu[0].prim[i].pdg==111) countpi++;
197  }
198  if (countpi > 0 && sr->mc.nu[0].iscc) return true;
199  return false;
200  });
201 
202 
203  //Define charged current no-pi0 background events.
204  const Cut CCNonPi0Bkg([](const caf::SRProxy* sr)
205  {
206  if(sr->mc.nnu == 0) return false;
207  if(sr->vtx.nelastic == 0) return false;
208  if(sr->mc.nu[0].prim.size() == 0) return false;
209  int nbofprim=sr->mc.nu[0].prim.size();
210  int countpi=0;
211  // bool ispi0=false;
212  for(int i=0;i<nbofprim;i++)
213  {
214  if(sr->mc.nu[0].prim[i].pdg==111) countpi++;
215  }
216  if (countpi == 0 && sr->mc.nu[0].iscc) return true;
217  return false;
218  });
219 
220 
221  //Create a cut based on the ReMId of the particle.
222  const Cut kremid([](const caf::SRProxy* sr)
223  { //May need to change the pid value.
224  //return (sr->sel.remid.pid < 0.35);
225  return (sr->sel.remid.pid < 0.4125);
226  });
227 
228 
229  //Ensure that the total energy of the pi0 is greater than .5MeV.
230  //May need to change the energy cut.
231  const Cut krecoprongenergy([](const caf::SRProxy* sr)
232  {
233  float MassOfPi0=0.135;
234  float kinetic=-10;
235  float en = -99;
236  float totEn;
237  float reco;
238  if(sr->vtx.nelastic == 0 && sr->mc.nnu==0) return false;
239  totEn=(sr->vtx.elastic[0].fuzzyk.png[0].shwlid.calE);
240  reco= (totEn-MassOfPi0);
241  // if (reco > 0.5) return true;
242  if(reco > .2) return true;
243  return false;
244  });
245 
246  //Make a cut to specify the total energy of the slice.
247  const Cut ksliceenergy([](const caf::SRProxy* sr)
248  {
249  if(sr->vtx.nelastic == 0 && sr->mc.nnu==0) return false;
250  if(sr->slc.calE > .2) return true;
251  return false;
252  });
253 
254  //Ensure that prongs have at least 5 hits in both the x and y views
255  const Cut kprongnhits([](const caf::SRProxy* sr)
256  {
257  if(sr->vtx.nelastic == 0 && sr->mc.nnu==0) return false;
258  if(sr->vtx.elastic[0].fuzzyk.png[0].nhitx<5) return false;
259  if(sr->vtx.elastic[0].fuzzyk.png[0].nhity<5) return false;
260  return true;
261  });
262 
263 
264  const Cut kElectron([](const caf::SRProxy* sr)
265  {
266  if(sr->vtx.nelastic == 0) return false;
267  if(sr->mc.nnu==0) return false;
268  if(sr->vtx.elastic[0].fuzzyk.npng == 0) return false;
269  if(sr->vtx.elastic[0].fuzzyk.png[0].truth.pdg==11) return true;
270  return false;
271  });
272 
273  const Cut kMuon([](const caf::SRProxy* sr)
274  {
275  if(sr->vtx.nelastic == 0) return false;
276  if(sr->mc.nnu==0) return false;
277  if(sr->vtx.elastic[0].fuzzyk.npng == 0) return false;
278  if(sr->vtx.elastic[0].fuzzyk.png[0].truth.pdg==13) return true;
279  return false;
280  });
281 
282  const Cut kPhoton([](const caf::SRProxy* sr)
283  {
284  if(sr->vtx.nelastic == 0) return false;
285  if(sr->mc.nnu==0) return false;
286  if(sr->vtx.elastic[0].fuzzyk.npng == 0) return false;
287  if(sr->vtx.elastic[0].fuzzyk.png[0].truth.pdg==22) return true;
288  return false;
289  });
290 
291  const Cut kPiPlus([](const caf::SRProxy* sr)
292  {
293  if(sr->vtx.nelastic == 0) return false;
294  if(sr->mc.nnu==0) return false;
295  if(sr->vtx.elastic[0].fuzzyk.npng == 0) return false;
296  if(sr->vtx.elastic[0].fuzzyk.png[0].truth.pdg==211)return true;
297  return false;
298  });
299 
300  const Cut kCVNPhoton([](const caf::SRProxy* sr)
301  {
302  return (sr->vtx.elastic[0].fuzzyk.png[0].cvnpart.photonid>0.7);
303  });
304 
305 
306 
307  /*****************************************************************************/
308  //Combine cuts to define signal and background
309 
310  const Cut Signal = kNCurrent && kPi0;
311  const Cut Bkg = !Signal;
314 
315 
316 
317 
318  /*****************************************************************************/
319 
320 }
caf::Proxy< size_t > npng
Definition: SRProxy.h:2038
caf::Proxy< unsigned int > nshwlid
Definition: SRProxy.h:2040
const Cut kCVNPhoton([](const caf::SRProxy *sr){return(sr->vtx.elastic[0].fuzzyk.png[0].cvnpart.photonid >0.7);})
const Cut kCCurrent([](const caf::SRProxy *sr){if(sr->vtx.nelastic==0) return false;if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);if(sr->vtx.elastic[0].fuzzyk.nshwlid==0) return false;if(sr->vtx.elastic[0].fuzzyk.npng==0) return false;return(sr->mc.nu[0].iscc==1);})
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Cut kElectron([](const caf::SRProxy *sr){if(sr->vtx.nelastic==0) return false;if(sr->mc.nnu==0) return false;if(sr->vtx.elastic[0].fuzzyk.npng==0) return false;if(sr->vtx.elastic[0].fuzzyk.png[0].truth.pdg==11) return true;return false;})
const Cut Signal
Definition: ncpi0Cuts.h:811
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const Cut k1Prong([](const caf::SRProxy *sr){if(sr->vtx.nelastic==0 &&sr->mc.nnu==0) return false;return(sr->vtx.elastic[0].fuzzyk.npng==1);})
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
T sqrt(T number)
Definition: d0nt_math.hpp:156
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
caf::Proxy< float > pid
Definition: SRProxy.h:1136
const Cut kNCurrent([](const caf::SRProxy *sr){if(sr->vtx.nelastic==0) return false;if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);if(sr->vtx.elastic[0].fuzzyk.nshwlid==0) return false;if(sr->vtx.elastic[0].fuzzyk.npng==0) return false;return(sr->mc.nu[0].iscc==0);})
CC and NC current cuts///.
double Distance(double x1, double y1, double x2, double y2)
const Cut kVertexDistance([](const caf::SRProxy *sr){if(sr->vtx.nelastic==0) return false;float Distance;Distance=sqrt( (sr->mc.nu[0].vtx.x-sr->vtx.elastic[0].vtx.x)* (sr->mc.nu[0].vtx.x-sr->vtx.elastic[0].vtx.x)+ (sr->mc.nu[0].vtx.y-sr->vtx.elastic[0].vtx.y)* (sr->mc.nu[0].vtx.y-sr->vtx.elastic[0].vtx.y)+ (sr->mc.nu[0].vtx.z-sr->vtx.elastic[0].vtx.z)* (sr->mc.nu[0].vtx.z-sr->vtx.elastic[0].vtx.z) );return Distance > 10.0;})
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
const Cut kNumu([](const caf::SRProxy *sr){if(sr->vtx.nelastic==0) return false;if(sr->mc.nnu==0) return false;return(sr->mc.nu[0].pdg==14);})
const Cut Preselection
const Cut kPhoton([](const caf::SRProxy *sr){if(sr->vtx.nelastic==0) return false;if(sr->mc.nnu==0) return false;if(sr->vtx.elastic[0].fuzzyk.npng==0) return false;if(sr->vtx.elastic[0].fuzzyk.png[0].truth.pdg==22) return true;return false;})
const Cut NCBkg([](const caf::SRProxy *sr){float MassOfPi0=0.135;float kinetic=-10;float en=-99;if(sr->mc.nnu==0) return false;if(sr->vtx.nelastic==0) return false;if(sr->mc.nu[0].prim.size()==0) return false;int nbofprim=sr->mc.nu[0].prim.size();int countpi=0;for(int i=0;i< nbofprim;i++){if(sr->mc.nu[0].prim[i].pdg==111){en=sr->mc.nu[0].prim[i].p.E;kinetic=en-MassOfPi0;countpi++;}}if(countpi==0 &&!sr->mc.nu[0].iscc) return true;return false;})
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
const Cut kPiPlus([](const caf::SRProxy *sr){if(sr->vtx.nelastic==0) return false;if(sr->mc.nnu==0) return false;if(sr->vtx.elastic[0].fuzzyk.npng==0) return false;if(sr->vtx.elastic[0].fuzzyk.png[0].truth.pdg==211) return true;return false;})
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2043
const Cut CCPi0Bkg([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;if(sr->vtx.nelastic==0) return false;if(sr->mc.nu[0].prim.size()==0) return false;int nbofprim=sr->mc.nu[0].prim.size();int countpi=0;for(int i=0;i< nbofprim;i++){if(sr->mc.nu[0].prim[i].pdg==111) countpi++;}if(countpi > 0 &&sr->mc.nu[0].iscc) return true;return false;})
caf::Proxy< float > z
Definition: SRProxy.h:108
if(dump)
caf::Proxy< float > x
Definition: SRProxy.h:106
const Cut Energy
const Cut kTrueFiducial
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
const Cut Bkg
caf::StandardRecord * sr
const Cut kFiducial([](const caf::SRProxy *sr){if(sr->vtx.nelastic==0 &&sr->mc.nnu==0) return false;if((sr->vtx.elastic[0].vtx.x)< -170.0) return false;if((sr->vtx.elastic[0].vtx.x) > 170.0) return false;if((sr->vtx.elastic[0].vtx.y)< -160.0) return false;if((sr->vtx.elastic[0].vtx.y) > 170.0) return false;if((sr->vtx.elastic[0].vtx.z)< 50.0) return false;if((sr->vtx.elastic[0].vtx.z) > 1100.0) return false;return true;})
caf::Proxy< caf::SRRemid > remid
Definition: SRProxy.h:1269
const Cut kContainment([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;if(sr->vtx.nelastic==0) return false;if(sr->vtx.elastic[0].fuzzyk.nshwlid==0) return false;if(sr->vtx.elastic[0].fuzzyk.npng==0) return false;for(unsigned int i=0;i< sr->vtx.elastic[0].fuzzyk.nshwlid;i++){if((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.x)< -160.0) return false;if((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.x) > 120.0) return false;if((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.y)< -180.0) return false;if((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.y) > 155.0) return false;if((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.z)< 200.0) return false;if((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.z) > 1200.0) return false;}return true;})
const Cut kNue([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;if(sr->vtx.nelastic==0) return false;return(sr->mc.nu[0].iscc && sr->mc.nu[0].pdg==12);})
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
const Cut kremid([](const caf::SRProxy *sr){ return(sr->sel.remid.pid< 0.4125);})
const Cut CCNonPi0Bkg([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;if(sr->vtx.nelastic==0) return false;if(sr->mc.nu[0].prim.size()==0) return false;int nbofprim=sr->mc.nu[0].prim.size();int countpi=0;for(int i=0;i< nbofprim;i++){if(sr->mc.nu[0].prim[i].pdg==111) countpi++;}if(countpi==0 &&sr->mc.nu[0].iscc) return true;return false;})
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2142
caf::Proxy< float > y
Definition: SRProxy.h:107
caf::Proxy< float > calE
Definition: SRProxy.h:1292
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2073
const Cut kVtx([](const caf::SRProxy *sr){if(sr->vtx.nelastic==0 &&sr->mc.nnu==0) return false;return(sr->vtx.nelastic==0);})
const Cut kprongnhits([](const caf::SRProxy *sr){if(sr->vtx.nelastic==0 &&sr->mc.nnu==0) return false;if(sr->vtx.elastic[0].fuzzyk.png[0].nhitx< 5) return false;if(sr->vtx.elastic[0].fuzzyk.png[0].nhity< 5) return false;return true;})
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
Template for Cut and SpillCut.
Definition: Cut.h:15
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
const Cut kMuon([](const caf::SRProxy *sr){if(sr->vtx.nelastic==0) return false;if(sr->mc.nnu==0) return false;if(sr->vtx.elastic[0].fuzzyk.npng==0) return false;if(sr->vtx.elastic[0].fuzzyk.png[0].truth.pdg==13) return true;return false;})
const Cut ksliceenergy([](const caf::SRProxy *sr){if(sr->vtx.nelastic==0 &&sr->mc.nnu==0) return false;if(sr->slc.calE >.2) return true;return false;})
#define for
Definition: msvc_pragmas.h:3
const Cut kPi0([](const caf::SRProxy *sr){float MassOfPi0=0.135;float kinetic=-10;float en=-99;if(sr->vtx.nelastic==0) return false;if(sr->vtx.elastic[0].fuzzyk.nshwlid==0) return false;if(sr->vtx.elastic[0].fuzzyk.npng==0) return false;if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);if(sr->mc.nu[0].prim.size()==0) return false;int nbofprim=sr->mc.nu[0].prim.size();int countpi0=0;for(int i=0;i< nbofprim;i++){if(sr->mc.nu[0].prim[i].pdg==111){en=sr->mc.nu[0].prim[i].p.E; if(en > 0.3){countpi0++;}}}if(countpi0 > 0) return true;return false;})
const Cut krecoprongenergy([](const caf::SRProxy *sr){float MassOfPi0=0.135;float kinetic=-10;float en=-99;float totEn;float reco;if(sr->vtx.nelastic==0 &&sr->mc.nnu==0) return false;totEn=(sr->vtx.elastic[0].fuzzyk.png[0].shwlid.calE);reco=(totEn-MassOfPi0);if(reco >.2) return true;return false;})