NusVars.cxx
Go to the documentation of this file.
2 #include "NuXAna/Vars/NusVars.h"
3 
5 
6 
7 namespace ana
8 {
9 
10  // Nus20 base variables for selection
11  const Var kIsElastic = SIMPLEVAR(vtx.elastic.IsValid);
12  const Var kNFuzzyProng = SIMPLEVAR(vtx.elastic.fuzzyk.npng);
13  const Var kNContPlanes = SIMPLEVAR(slc.ncontplanes);
14  const Var kCVNnc_looseptp = SIMPLEVAR(sel.cvnloosepreselptp.ncid);
15  const Var kCVNnc_oldpresel = SIMPLEVAR(sel.cvnoldpresel.ncid);
16 
17  // Nus20 BDT variables
18  const Var kNSliceHits = SIMPLEVAR(slc.nhit);
19  const Var kShwNHitX = SIMPLEVAR(vtx.elastic.fuzzyk.png[0].shwlid.nhitx);
20  const Var kShwNHitY = SIMPLEVAR(vtx.elastic.fuzzyk.png[0].shwlid.nhity);
21  const Var kShwDirY = SIMPLEVAR(vtx.elastic.fuzzyk.png[0].shwlid.dir.y);
22  const Var kNShwLid = SIMPLEVAR(vtx.elastic.fuzzyk.nshwlid);
23  const Var kNMipHit = SIMPLEVAR(slc.nmiphit);
24 
25  const Var kShwNHitXPlusY([](const caf::SRProxy* sr)
26  {
27  double nHits = sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx +
28  sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;
29 
30  return nHits;
31  });
32 
33  const Var kShwNHitXMinY([](const caf::SRProxy* sr)
34  {
35  double nHits = sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx -
36  sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;
37 
38  return nHits;
39  });
40 
41  const Var kShwNHitXOverY([](const caf::SRProxy* sr)
42  {
43  double hitRatio = (double)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx /
44  (double)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;
45 
46  return hitRatio;
47  });
48 
49  // Nus17 reconstructed (visible) energy - see docdb 17797 and 17985
50  const double FDscaleCalE17 = 0.6377;
51  const double NDscaleCalE17 = 0.6484;
52 
53  const Var kNus17Energy([](const caf::SRProxy* sr)
54  {
55  double cale = sr->slc.calE;
56  double recoE = 0.;
57  if(sr->hdr.det == caf::kFARDET) recoE = FDscaleCalE17*cale;
58  if(sr->hdr.det == caf::kNEARDET) recoE = NDscaleCalE17*cale;
59  return recoE;
60  });
61 
62 
63  const double NDscaleCalE18RHC = 1.15;
64  const double NDscaleCalE18 = 1.11;
65  const double FDscaleCalE18RHC = 1.18;
66  const double FDscaleCalE18 = 1.20;
67 
68  const Var kNus18Energy([](const caf::SRProxy* sr) {
69  bool h_FHC = sr->spill.isFHC;
70  bool h_RHC = sr->spill.isRHC;
71  double cale = sr->slc.calE;
72  double recoE = 0.;
73  if (h_FHC && sr->hdr.det == caf::kFARDET) recoE = FDscaleCalE18*cale;
74  else if (h_FHC && sr->hdr.det == caf::kNEARDET) recoE = NDscaleCalE18*cale;
75  else if (h_RHC && sr->hdr.det == caf::kFARDET) recoE = FDscaleCalE18RHC*cale;
76  else if (h_RHC && sr->hdr.det == caf::kNEARDET) recoE = NDscaleCalE18RHC*cale;
77  return recoE ; // default case
78  });
79 
80 
81  const Var kNus20Energy([](const caf::SRProxy* sr) {
82  if (sr->hdr.det != caf::kNEARDET && sr->hdr.det != caf::kFARDET) return (double) sr->slc.calE; // Only defined for ND & FD (eg, not test beam)
83  double pars[4][6] = { {1.049, 0.795, 0.8409, 0.17, 0.82, -1.00}, // ND FHC (or 0HC)
84  {1.025, 0.797, 0.9162, 0.53, -0.26, -1.48}, // FD FHC (or 0HC)
85  {1.000, 1.000, 1.0000, 0.00, 0.00, 0.00}, // ND RHC (just return calE = e_EM + e_Had)
86  {1.000, 1.000, 1.0000, 0.00, 0.00, 0.00} };// FD RHC (just return calE = e_EM + e_Had)
87  int detCur = (sr->hdr.det == caf::kFARDET) + ((sr->spill.isRHC)<<1); // bitshifting is magic
88  double e_EM = ana:: kEME_2020(sr);
89  double e_Had = ana::kHADE_2020(sr);
90  double e_Cal = sr->slc.calE;
91  return (e_EM / pars[detCur][0] + e_Had / pars[detCur][1]) / (pars[detCur][2] + pars[detCur][3] * std::pow(e_Cal+pars[detCur][4],2) * std::exp(pars[detCur][5] * e_Cal));
92  });
93 
94  /// NC Cosmic Rejection (G)BDT variables
95  const Var kNCCosRej = SIMPLEVAR(sel.nccosrej.cospidbdt);
96  const Var kNCCosRejG = SIMPLEVAR(sel.nccosrej.cospidbdtg);
97 
98  const Var kPartPtp([](const caf::SRProxy* sr) {
99  if(std::isnan(1.*sr->sel.nuecosrej.partptp) ||
100  std::isinf(1.*sr->sel.nuecosrej.partptp)) return -5.f;
101  return float(sr->sel.nuecosrej.partptp);
102  });
103 
104  const Var kDistTop([](const caf::SRProxy* sr) {
105  return std::min(sr->sel.nuecosrej.starttop,
106  sr->sel.nuecosrej.stoptop);
107  });
108 
109  const Var kMinDistFace([](const caf::SRProxy* sr) {
110  const caf::SRNueCosRejProxy& cr = sr->sel.nuecosrej;
111  float mindist = 10000.;
112  mindist = std::min(mindist, std::min(cr.starteast, cr.stopeast));
113  mindist = std::min(mindist, std::min(cr.startwest, cr.stopwest));
114  mindist = std::min(mindist, std::min(cr.starttop, cr.stoptop));
115  mindist = std::min(mindist, std::min(cr.startbottom, cr.stopbottom));
116  mindist = std::min(mindist, std::min(cr.startfront, cr.stopfront));
117  mindist = std::min(mindist, std::min(cr.startback, cr.stopback));
118  return (mindist > 9999. ? (float)-5. : mindist);
119  });
120 
121  const Var kClosestSlcTime([](const caf::SRProxy* sr) {
122  return (float)sr->slc.closestslicetime;
123  });
124 
125  const Var kClosestSlcMinDist([](const caf::SRProxy* sr) {
126  return (float)sr->slc.closestslicemindist;
127  });
128 
129  const Var kClosestSlcMinTop([](const caf::SRProxy* sr) {
130  return (float)sr->slc.closestsliceminfromtop;
131  });
132 
133  //---------------------------------------------------------------------------------
134  // Nus18 kaon
136  if (abs(nu->beam.ptype) == 321 || nu->beam.ptype == 310 || nu->beam.ptype == 311 || nu->beam.ptype == 130)
137  return 1.3;
138  return 1.;
139  });
140 
142  if (abs(nu->beam.ptype) == 321 || nu->beam.ptype == 310 || nu->beam.ptype == 311 || nu->beam.ptype == 130)
143  return 0.7;
144  return 1.;
145  });
146 
147  //---------------------------------------------------------------------------------
148  // Nus18 tau
149  const NuTruthVar kFDTauWgtM_NT([](const caf::SRNeutrinoProxy* nu) {
150  if (abs(nu->pdg) == 16 && nu->iscc)
151  return 0.4;
152  return 1.;
153  });
154 
155  const NuTruthVar kFDTauWgtP_NT([](const caf::SRNeutrinoProxy* nu) {
156  if (abs(nu->pdg) == 16 && nu->iscc)
157  return 1.6;
158  return 1.;
159  });
160 }
caf::Proxy< caf::SRSpill > spill
Definition: SRProxy.h:2143
Near Detector underground.
Definition: SREnums.h:10
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
const double FDscaleCalE18
Definition: NusVars.cxx:66
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 Var kShwDirY
Definition: NusVars.cxx:21
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
Far Detector at Ash River.
Definition: SREnums.h:11
const NuTruthVar kNus18BeamWeightP_NT([](const caf::SRNeutrinoProxy *nu){if(abs(nu->beam.ptype)==321||nu->beam.ptype==310||nu->beam.ptype==311||nu->beam.ptype==130) return 1.3;return 1.;})
Definition: NusVars.h:86
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kNus18Energy([](const caf::SRProxy *sr){bool h_FHC=sr->spill.isFHC;bool h_RHC=sr->spill.isRHC;double cale=sr->slc.calE;double recoE=0.;if(h_FHC &&sr->hdr.det==caf::kFARDET) recoE=FDscaleCalE18 *cale;else if(h_FHC &&sr->hdr.det==caf::kNEARDET) recoE=NDscaleCalE18 *cale;else if(h_RHC &&sr->hdr.det==caf::kFARDET) recoE=FDscaleCalE18RHC *cale;else if(h_RHC &&sr->hdr.det==caf::kNEARDET) recoE=NDscaleCalE18RHC *cale;return recoE;})
Definition: NusVars.h:63
const Var kShwNHitX
Definition: NusVars.cxx:19
caf::Proxy< float > stoptop
Definition: SRProxy.h:1080
caf::Proxy< float > stopfront
Definition: SRProxy.h:1079
const Var kDistTop([](const caf::SRProxy *sr){return std::min(sr->sel.nuecosrej.starttop, sr->sel.nuecosrej.stoptop);})
Definition: NusVars.h:76
const Var kIsElastic
Definition: NusVars.cxx:11
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
const Var kShwNHitXOverY([](const caf::SRProxy *sr){double hitRatio=(double) sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx/(double) sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;return hitRatio;})
Definition: NusVars.h:50
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const Var kNus20Energy([](const caf::SRProxy *sr){if(sr->hdr.det!=caf::kNEARDET &&sr->hdr.det!=caf::kFARDET) return(double) sr->slc.calE;double pars[4][6]={{1.049, 0.795, 0.8409, 0.17, 0.82,-1.00},{1.025, 0.797, 0.9162, 0.53,-0.26,-1.48},{1.000, 1.000, 1.0000, 0.00, 0.00, 0.00},{1.000, 1.000, 1.0000, 0.00, 0.00, 0.00}};int detCur=(sr->hdr.det==caf::kFARDET)+((sr->spill.isRHC)<< 1);double e_EM=ana::kEME_2020(sr);double e_Had=ana::kHADE_2020(sr);double e_Cal=sr->slc.calE;return(e_EM/pars[detCur][0]+e_Had/pars[detCur][1])/(pars[detCur][2]+pars[detCur][3]*std::pow(e_Cal+pars[detCur][4], 2)*std::exp(pars[detCur][5]*e_Cal));})
Definition: NusVars.h:64
constexpr T pow(T x)
Definition: pow.h:75
caf::Proxy< float > startbottom
Definition: SRProxy.h:1071
void abs(TH1 *hist)
const Var kNCCosRejG
Definition: NusVars.cxx:96
const Var kNContPlanes
Definition: NusVars.cxx:13
const Var kMinDistFace([](const caf::SRProxy *sr){const caf::SRNueCosRejProxy &cr=sr->sel.nuecosrej;float mindist=10000.;mindist=std::min(mindist, std::min(cr.starteast, cr.stopeast));mindist=std::min(mindist, std::min(cr.startwest, cr.stopwest));mindist=std::min(mindist, std::min(cr.starttop, cr.stoptop));mindist=std::min(mindist, std::min(cr.startbottom, cr.stopbottom));mindist=std::min(mindist, std::min(cr.startfront, cr.stopfront));mindist=std::min(mindist, std::min(cr.startback, cr.stopback));return(mindist > 9999.?(float)-5.:mindist);})
Definition: NusVars.h:78
caf::Proxy< float > starttop
Definition: SRProxy.h:1074
caf::Proxy< float > stopback
Definition: SRProxy.h:1076
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
const Var kNFuzzyProng
Definition: NusVars.cxx:12
const Var kNus17Energy([](const caf::SRProxy *sr){double cale=sr->slc.calE;double recoE=0.;if(sr->hdr.det==caf::kFARDET) recoE=FDscaleCalE17 *cale;if(sr->hdr.det==caf::kNEARDET) recoE=NDscaleCalE17 *cale;return recoE;})
Definition: NusVars.h:62
caf::Proxy< float > closestslicetime
Definition: SRProxy.h:1302
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
caf::Proxy< float > startback
Definition: SRProxy.h:1070
caf::Proxy< caf::SRNueCosRej > nuecosrej
Definition: SRProxy.h:1265
const NuTruthVar kFDTauWgtM_NT([](const caf::SRNeutrinoProxy *nu){if(abs(nu->pdg)==16 &&nu->iscc) return 0.4;return 1.;})
Definition: NusVars.h:95
caf::Proxy< bool > isFHC
Definition: SRProxy.h:1375
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2043
const NuTruthVar kNus18BeamWeightM_NT([](const caf::SRNeutrinoProxy *nu){if(abs(nu->beam.ptype)==321||nu->beam.ptype==310||nu->beam.ptype==311||nu->beam.ptype==130) return 0.7;return 1.;})
Definition: NusVars.h:89
caf::Proxy< float > stopeast
Definition: SRProxy.h:1078
if(dump)
caf::Proxy< float > stopbottom
Definition: SRProxy.h:1077
const Var kNMipHit
Definition: NusVars.cxx:23
caf::Proxy< float > startfront
Definition: SRProxy.h:1073
caf::Proxy< bool > iscc
Definition: SRProxy.h:538
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
_Var< caf::SRNeutrinoProxy > NuTruthVar
Var designed to be used over the nuTree, ie all neutrinos, not just those that got slices...
Definition: Var.h:82
caf::StandardRecord * sr
caf::Proxy< float > partptp
Definition: SRProxy.h:1057
caf::Proxy< float > stopwest
Definition: SRProxy.h:1081
caf::Proxy< float > closestsliceminfromtop
Definition: SRProxy.h:1299
const Var kNSliceHits
Definition: NusVars.cxx:18
caf::Proxy< caf::SRBeam > beam
Definition: SRProxy.h:523
const NuTruthVar kFDTauWgtP_NT([](const caf::SRNeutrinoProxy *nu){if(abs(nu->pdg)==16 &&nu->iscc) return 1.6;return 1.;})
Definition: NusVars.h:98
const Var kClosestSlcMinTop([](const caf::SRProxy *sr){return(float) sr->slc.closestsliceminfromtop;})
Definition: NusVars.h:82
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
caf::Proxy< float > starteast
Definition: SRProxy.h:1072
const double NDscaleCalE17
Definition: NusVars.cxx:51
const Var kCVNnc_oldpresel
Definition: NusVars.cxx:15
const Var kClosestSlcMinDist([](const caf::SRProxy *sr){return(float) sr->slc.closestslicemindist;})
Definition: NusVars.h:81
const Var kShwNHitXPlusY([](const caf::SRProxy *sr){double nHits=sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx+sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;return nHits;})
Definition: NusVars.h:49
const Var kNShwLid
Definition: NusVars.cxx:22
const Var kCVNnc_looseptp
Definition: NusVars.cxx:14
caf::Proxy< int > ptype
Definition: SRProxy.h:366
std::string pars("Th23Dmsq32")
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2142
caf::Proxy< float > calE
Definition: SRProxy.h:1292
caf::Proxy< float > closestslicemindist
Definition: SRProxy.h:1294
Proxy for caf::SRNueCosRej.
Definition: SRProxy.h:1024
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
const Var kShwNHitXMinY([](const caf::SRProxy *sr){double nHits=sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx-sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;return nHits;})
Definition: NusVars.h:43
const Var kShwNHitY
Definition: NusVars.cxx:20
caf::Proxy< short int > pdg
Definition: SRProxy.h:552
caf::Proxy< bool > isRHC
Definition: SRProxy.h:1376
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
T min(const caf::Proxy< T > &a, T b)
const double FDscaleCalE17
Definition: NusVars.cxx:50
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
caf::Proxy< float > startwest
Definition: SRProxy.h:1075
const double NDscaleCalE18RHC
Definition: NusVars.cxx:63
const Var kNCCosRej
NC Cosmic Rejection (G)BDT variables.
Definition: NusVars.cxx:95
const Var kClosestSlcTime([](const caf::SRProxy *sr){return(float) sr->slc.closestslicetime;})
Definition: NusVars.h:80
const double NDscaleCalE18
Definition: NusVars.cxx:64
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
const double FDscaleCalE18RHC
Definition: NusVars.cxx:65
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:232