NueEnergy2020.cxx
Go to the documentation of this file.
2 #include "CAFAna/Vars/Vars.h"
4 
6 
7 namespace ana
8 {
9 
10  double NueRecoE_2020RHCFit_2D3D(double rawEM, double rawHA)
11  {
12  if(rawEM < 0. || rawHA < 0.) return -5.;
13  const double p0 = 0.0;
14  const double p1 = 9.54388e-01;
15  const double p2 = 1.13311;
16  const double p3 = 0.0;
17  const double p4 = 3.39617e-06;
18  const double p5 = 3.49874e-06;
19  const double recoEn = (1./(1+0.0089272))*(rawHA*rawHA*p5 + rawEM*rawEM*p4 + rawEM*rawHA*p3 + rawHA*p2 + rawEM*p1 + p0);
20 
21  return recoEn;
22  }
23 
24  double NueRecoE_2020RHCFit(double rawEM, double rawHA)
25  {
26  if(rawEM < 0. || rawHA < 0.) return -5.;
27  const double p0 = 0.0;
28  const double p1 = 9.88258e-01;
29  const double p2 = 1.20084;
30  const double p3 = 0.0;
31  const double p4 = 1.92904e-07;
32  const double p5 = 1.20704e-07;
33  const double recoEn = (1./(1+0.0111873))*(rawHA*rawHA*p5 + rawEM*rawEM*p4 + rawEM*rawHA*p3 + rawHA*p2 + rawEM*p1 + p0);
34 
35  return recoEn;
36  }
37 
38  double NueRecoE_2020FHCFit_2D3D(double rawEM, double rawHA)
39  {
40  if(rawEM < 0. || rawHA < 0.) return -5.;
41  const double p0 = 0.0;
42  const double p1 = 9.24605e-01;
43  const double p2 = 1.15426;
44  const double p3 = 0.0;
45  const double p4 = 1.58171e-02;
46  const double p5 = 6.81363e-03;
47  const double recoEn = (1./(1+0.0217326))*(rawHA*rawHA*p5 + rawEM*rawEM*p4 + rawEM*rawHA*p3 + rawHA*p2 + rawEM*p1 + p0);
48  return recoEn;
49  }
50 
51  double NueRecoE_2020FHCFit(double rawEM, double rawHA)
52  {
53  if(rawEM < 0. || rawHA < 0.) return -5.;
54  const double p0 = 0.0;
55  const double p1 = 1.01777;
56  const double p2 = 1.10868;
57  const double p3 = 0.0;
58  const double p4 = 1.43541e-03;
59  const double p5 = 1.09628e-01;
60  const double recoEn = (1./(1+0.0355622))*(rawHA*rawHA*p5 + rawEM*rawEM*p4 + rawEM*rawHA*p3 + rawHA*p2 + rawEM*p1 + p0);
61  return recoEn;
62  }
63  const Var kEME2020_2D3D([](const caf::SRProxy* sr)
64  {
65  if(!sr->vtx.elastic.IsValid) return -1.0;
66  if(sr->vtx.elastic.fuzzyk.npng < 1) return -1.0;
67  const caf::SRFuzzyKProngProxy& prim_png = sr->vtx.elastic.fuzzyk.png[0];
68 
69  double CVNem_CalE = 0.0;
70  for (const caf::SRFuzzyKProngProxy& png: sr->vtx.elastic.fuzzyk.png ){
71  double png_CalE = png.shwlid.calE;
72  double emPID = ( (double)png.cvnpart.photonid +(double)png.cvnpart.electronid);
73 
74  if ( emPID <= 0 ) continue;
75  if ( emPID >= 0.5 ) CVNem_CalE += png_CalE;
76  }//ipng
77 
78  for (const caf::SRProngProxy& png: sr->vtx.elastic.fuzzyk.png2d ){
79  double png_CalE = png.calE;
80  double emPID = ( (double)png.cvnpart.photonid + (double)png.cvnpart.electronid);
81 
82  if ( emPID <= 0 ) continue;
83  if ( emPID >= 0.7 ) CVNem_CalE += png_CalE;
84  }
85  if(CVNem_CalE <=0 ||kLongestProng(sr) >= 500) CVNem_CalE = prim_png.shwlid.calE;
86  return CVNem_CalE;
87  });
88 
89  const Var kHADE2020_2D3D([](const caf::SRProxy* sr)
90  {
91  if(!sr->vtx.elastic.IsValid) return -1.0;
92  if(sr->vtx.elastic.fuzzyk.npng < 1) return -1.0;
93 
94  double CVNhad_CalE = 0.0;
95  double em = 0.0;
96  for (const caf::SRFuzzyKProngProxy& png: sr->vtx.elastic.fuzzyk.png ){
97  double png_CalE = png.shwlid.calE;
98  double emPID = ( (double)png.cvnpart.photonid +(double)png.cvnpart.electronid);
99 
100  if ( emPID < 0.5 ) CVNhad_CalE += png_CalE;
101  else em +=png_CalE;
102  }
103 
104  for (const caf::SRProngProxy& png: sr->vtx.elastic.fuzzyk.png2d ){
105  double png_CalE = png.calE;
106  double emPID = ( (double)png.cvnpart.photonid + (double)png.cvnpart.electronid);
107 
108  if ( emPID < 0.7 ) CVNhad_CalE += png_CalE;
109  else em += png_CalE;
110  }
111  CVNhad_CalE += sr->vtx.elastic.fuzzyk.orphCalE;
112  if (em<=0 || kLongestProng(sr) >= 500) CVNhad_CalE -= sr->vtx.elastic.fuzzyk.png[0].shwlid.calE;
113  if (CVNhad_CalE<0) abort();
114  return CVNhad_CalE;
115 
116  });
117 
118 
119 
120  const Var kEME_2020([](const caf::SRProxy* sr)
121  {
122  if(!sr->vtx.elastic.IsValid) return -1.0;
123  if(sr->vtx.elastic.fuzzyk.npng < 1) return -1.0;
124  double CVNem_CalE = 0.0;
125  const caf::SRFuzzyKProngProxy& prim_png = sr->vtx.elastic.fuzzyk.png[0];
126  for (const caf::SRFuzzyKProngProxy& png: sr->vtx.elastic.fuzzyk.png ){
127  double png_CalE = png.shwlid.calE;
128 
129  double emPID = ( (double)png.cvnpart.photonid +
130  (double)png.cvnpart.pizeroid +
131  (double)png.cvnpart.electronid);
132  double haPID = ( (double)png.cvnpart.protonid +
133  (double)png.cvnpart.pionid +
134  (double)png.cvnpart.neutronid +
135  (double)png.cvnpart.otherid +
136  (double)png.cvnpart.muonid);
137 
138  if ( emPID < 0 ) continue;
139  if ( emPID >= haPID ) CVNem_CalE += png_CalE;
140 
141  else continue;
142  }
143  if (kLongestProng(sr) >= 500 ) CVNem_CalE = prim_png.shwlid.calE;
144  return CVNem_CalE;
145  });
146 
147  const Var kHADE_2020([](const caf::SRProxy* sr)
148  {
149  if(!sr->vtx.elastic.IsValid) return -1.0;
150  if(sr->vtx.elastic.fuzzyk.npng < 1) return -1.0;
151 
152  const double CVNha_CalE = sr->slc.calE;
153  return std::max (CVNha_CalE-kEME_2020(sr), 0.);
154  });
155 
157  [](const caf::SRProxy* sr)
158  {
159  if(!sr->vtx.elastic.IsValid) return -1.0;
160  if(sr->vtx.elastic.fuzzyk.npng < 1) return -1.0;
161 
163  });
164 
165  const Var kNueEnergyFHC2020(
166  [](const caf::SRProxy* sr)
167  {
168  if(!sr->vtx.elastic.IsValid) return -1.0;
169  if(sr->vtx.elastic.fuzzyk.npng < 1) return -1.0;
170 
171  return NueRecoE_2020FHCFit(kEME_2020(sr), kHADE_2020(sr));
172  });
173 
174 
176  [](const caf::SRProxy* sr)
177  {
178  if(!sr->vtx.elastic.IsValid) return -1.0;
179  if(sr->vtx.elastic.fuzzyk.npng < 1) return -1.0;
181  });
182  const Var kNueEnergyRHC2020(
183  [](const caf::SRProxy* sr)
184  {
185  if(!sr->vtx.elastic.IsValid) return -1.0;
186  if(sr->vtx.elastic.fuzzyk.npng < 1) return -1.0;
187  return NueRecoE_2020RHCFit(kEME_2020(sr), kHADE_2020(sr));
188  });
189  /// Nue energy with 3d prong and 2d prong info. (new version of vars)
190  const Var kNueEnergy2020_2D3D(
191  [](const caf::SRProxy* sr)
192  {
193  if(!sr->vtx.elastic.IsValid) return -1.0;
194  if(sr->vtx.elastic.fuzzyk.npng < 1) return -1.0;
195  if(kIsRHC(sr)) return kNueEnergyRHC2020_2D3D(sr);
196  else return kNueEnergyFHC2020_2D3D(sr);
197  });
198  /// Nue energy with 3d prong info. only (old version of vars)
199  const Var kNueEnergy2020(
200  [](const caf::SRProxy* sr)
201  {
202  if(!sr->vtx.elastic.IsValid) return -1.0;
203  if(sr->vtx.elastic.fuzzyk.npng < 1) return -1.0;
204  if(kIsRHC(sr)) return kNueEnergyRHC2020(sr);
205  else return kNueEnergyFHC2020(sr);
206  });
207 
208  // Reconstructed nue electron energy Ee from estimator tuned for 2020
209  const Var kNueElectronE2020(
210  []( const caf::SRProxy* sr )
211  {
212  if ( !sr->vtx.elastic.IsValid ) return -1000.f;
213  if ( sr->vtx.elastic.fuzzyk.npng < 1 ) return -1000.f;
214 
215  float x = sr->vtx.elastic.fuzzyk.png[0].calE;
216  float energy = x;
217 
218  double a0, a1, a2, a3, a4, a5, a6;
219 
220  // To-do: provide a DocDb reference when available
221  if ( sr->hdr.det == caf::kFARDET )
222  {
223  if ( sr->spill.isRHC )
224  {
225  a0 = -0.727211;
226  a1 = 2.319804;
227  a2 = -2.697100;
228  a3 = 1.604802;
229  a4 = -0.493659;
230  a5 = 0.075183;
231  a6 = -0.004478;
232  }
233  else
234  {
235  a0 = -0.428873;
236  a1 = 1.486957;
237  a2 = -1.772140;
238  a3 = 1.114639;
239  a4 = -0.358208;
240  a5 = 0.056065;
241  a6 = -0.003376;
242  }
243  }
244  // Sorry, no tune for ND yet
245  else
246  {
247  return energy;
248  }
249 
250  if ( x > 0 ) energy = x - ( a0 + a1*x + a2*pow(x,2) +
251  a3*pow(x,3) + a4*pow(x,4) +
252  a5*pow(x,5) + a6*pow(x,6) );
253 
254  if ( energy < 0 ) energy = x;
255 
256  return energy;
257  });
258 
259  // Reconstructed nue electron transverse momentum Pt using kNueElectronE2020
260  const Var kNueElectronPt2020(
261  []( const caf::SRProxy* sr )
262  {
263  if ( !sr->vtx.elastic.IsValid ) return -1000.f;
264  if ( sr->vtx.elastic.fuzzyk.npng < 1 ) return -1000.f;
265 
266  TVector3 elecDir = (TVector3)sr->vtx.elastic.fuzzyk.png[0].dir;
267  TVector3 beamDir = ana::NuMIBeamDirection( sr->hdr.det );
268  float ct = elecDir.Dot( beamDir );
269  float p = kNueElectronE2020( sr ); // Ignore electron mass
270  if ( p < 0 || fabs( ct ) > 1 ) return -1000.f;
271  return float( p * sin( acos( ct ) ) );
272  });
273 
274 }
275 
276 
caf::Proxy< size_t > npng
Definition: SRProxy.h:2020
T max(const caf::Proxy< T > &a, T b)
TH1F * a3
Definition: f2_nu.C:606
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2041
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TH1F * a2
Definition: f2_nu.C:545
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const Var kNueEnergyRHC2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;return NueRecoE_2020RHCFit(kEME_2020(sr), kHADE_2020(sr));})
Definition: NueEnergy2020.h:33
caf::Proxy< std::vector< caf::SRProng > > png2d
Definition: SRProxy.h:2026
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2119
const char * p
Definition: xmltok.h:285
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2108
constexpr T pow(T x)
Definition: pow.h:75
const Var kNueEnergy2020_2D3D([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;if(kIsRHC(sr)) return kNueEnergyRHC2020_2D3D(sr);else return kNueEnergyFHC2020_2D3D(sr);})
Nue energy with 3d prong and 2d prong info. (new version of vars)
Definition: NueEnergy2020.h:36
const Var kEME2020_2D3D([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;const caf::SRFuzzyKProngProxy &prim_png=sr->vtx.elastic.fuzzyk.png[0];double CVNem_CalE=0.0;for(const caf::SRFuzzyKProngProxy &png:sr->vtx.elastic.fuzzyk.png){double png_CalE=png.shwlid.calE;double emPID=((double) png.cvnpart.photonid+(double) png.cvnpart.electronid);if(emPID<=0) continue;if(emPID >=0.5 ) CVNem_CalE+=png_CalE;}for(const caf::SRProngProxy &png:sr->vtx.elastic.fuzzyk.png2d){double png_CalE=png.calE;double emPID=((double) png.cvnpart.photonid+(double) png.cvnpart.electronid);if(emPID<=0) continue;if(emPID >=0.7 ) CVNem_CalE+=png_CalE;}if(CVNem_CalE<=0||kLongestProng(sr) >=500) CVNem_CalE=prim_png.shwlid.calE;return CVNem_CalE;})
Definition: NueEnergy2020.h:18
T acos(T number)
Definition: d0nt_math.hpp:54
double NueRecoE_2020RHCFit_2D3D(double rawEM, double rawHA)
double NueRecoE_2020FHCFit(double rawEM, double rawHA)
const Var kHADE2020_2D3D([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;double CVNhad_CalE=0.0;double em=0.0;for(const caf::SRFuzzyKProngProxy &png:sr->vtx.elastic.fuzzyk.png){double png_CalE=png.shwlid.calE;double emPID=((double) png.cvnpart.photonid+(double) png.cvnpart.electronid);if(emPID< 0.5 ) CVNhad_CalE+=png_CalE;else em+=png_CalE;}for(const caf::SRProngProxy &png:sr->vtx.elastic.fuzzyk.png2d){double png_CalE=png.calE;double emPID=((double) png.cvnpart.photonid+(double) png.cvnpart.electronid);if(emPID< 0.7 ) CVNhad_CalE+=png_CalE;else em+=png_CalE;}CVNhad_CalE+=sr->vtx.elastic.fuzzyk.orphCalE;if(em<=0||kLongestProng(sr) >=500) CVNhad_CalE-=sr->vtx.elastic.fuzzyk.png[0].shwlid.calE;if(CVNhad_CalE< 0) abort();return CVNhad_CalE;})
Definition: NueEnergy2020.h:19
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2100
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2025
const Var kLongestProng([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return 0.f;if(sr->vtx.elastic.fuzzyk.npng==0) return 0.f;auto idx=sr->vtx.elastic.fuzzyk.longestidx;return float(sr->vtx.elastic.fuzzyk.png[idx].len);})
Definition: Vars.h:89
if(dump)
const Var kNueEnergyRHC2020_2D3D([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;return NueRecoE_2020RHCFit_2D3D(kEME2020_2D3D(sr), kHADE2020_2D3D(sr));})
Definition: NueEnergy2020.h:30
TH1F * a1
Definition: f2_nu.C:476
double energy
Definition: plottest35.C:25
caf::StandardRecord * sr
double NueRecoE_2020RHCFit(double rawEM, double rawHA)
const Var kNueElectronE2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1000.f;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1000.f;float x=sr->vtx.elastic.fuzzyk.png[0].calE;float energy=x;double a0, a1, a2, a3, a4, a5, a6;if(sr->hdr.det==caf::kFARDET){if(sr->spill.isRHC){a0=-0.727211;a1=2.319804;a2=-2.697100;a3=1.604802;a4=-0.493659;a5=0.075183;a6=-0.004478;}else{a0=-0.428873;a1=1.486957;a2=-1.772140;a3=1.114639;a4=-0.358208;a5=0.056065;a6=-0.003376;}}else{return energy;}if(x > 0) energy=x-(a0+a1 *x+a2 *pow(x, 2)+ a3 *pow(x, 3)+a4 *pow(x, 4)+ a5 *pow(x, 5)+a6 *pow(x, 6));if(energy< 0) energy=x;return energy;})
Definition: NueEnergy2020.h:41
const Var kHADE_2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;const double CVNha_CalE=sr->slc.calE;return std::max(CVNha_CalE-kEME_2020(sr), 0.);})
Definition: NueEnergy2020.h:15
const Cut kIsRHC([](const caf::SRProxy *sr){return sr->spill.isRHC;})
Definition: Vars.h:16
caf::Proxy< float > orphCalE
Definition: SRProxy.h:2024
const Var kNueEnergy2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;if(kIsRHC(sr)) return kNueEnergyRHC2020(sr);else return kNueEnergyFHC2020(sr);})
Nue energy with 3d prong info. only (old version of vars)
Definition: NueEnergy2020.h:38
T sin(T number)
Definition: d0nt_math.hpp:132
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2040
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2124
caf::Proxy< float > calE
Definition: SRProxy.h:1271
const Var kNueEnergyFHC2020_2D3D([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;return NueRecoE_2020FHCFit_2D3D(kEME2020_2D3D(sr), kHADE2020_2D3D(sr));})
Definition: NueEnergy2020.h:22
const Var kNueElectronPt2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1000.f;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1000.f;TVector3 elecDir=(TVector3) sr->vtx.elastic.fuzzyk.png[0].dir;TVector3 beamDir=ana::NuMIBeamDirection(sr->hdr.det);float ct=elecDir.Dot(beamDir);float p=kNueElectronE2020(sr);if(p< 0||fabs(ct) > 1) return-1000.f;return float(p *sin(acos(ct)));})
Definition: NueEnergy2020.h:43
TH1F * a4
Definition: f2_nu.C:666
TVector3 NuMIBeamDirection(caf::Det_t det)
Average direction of NuMI neutrinos in a given detector This function is a copy of geo::GeometryBase:...
Definition: Utilities.cxx:333
double NueRecoE_2020FHCFit_2D3D(double rawEM, double rawHA)
const Var kNueEnergyFHC2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;return NueRecoE_2020FHCFit(kEME_2020(sr), kHADE_2020(sr));})
Definition: NueEnergy2020.h:25
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2128
#define for
Definition: msvc_pragmas.h:3
const Var kEME_2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;double CVNem_CalE=0.0;const caf::SRFuzzyKProngProxy &prim_png=sr->vtx.elastic.fuzzyk.png[0];for(const caf::SRFuzzyKProngProxy &png:sr->vtx.elastic.fuzzyk.png){double png_CalE=png.shwlid.calE;double emPID=((double) png.cvnpart.photonid+ (double) png.cvnpart.pizeroid+ (double) png.cvnpart.electronid);double haPID=((double) png.cvnpart.protonid+ (double) png.cvnpart.pionid+ (double) png.cvnpart.neutronid+ (double) png.cvnpart.otherid+ (double) png.cvnpart.muonid);if(emPID< 0) continue;if(emPID >=haPID ) CVNem_CalE+=png_CalE;else continue;}if(kLongestProng(sr) >=500) CVNem_CalE=prim_png.shwlid.calE;return CVNem_CalE;})
Definition: NueEnergy2020.h:14
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:231