NumuCC2p2hVars.cxx
Go to the documentation of this file.
3 
4 
5 namespace ana
6 {
7  //gotta make the muon mass a var for reasons
8  inline double MuonMass()
9  {
10  static TDatabasePDG* pdgdb = TDatabasePDG::Instance();
11  static double mass = pdgdb->GetParticle(13)->Mass();
12  return mass;
13  }
14 
15  const Var kMuMass([](const caf::SRProxy* sr)
16  {
17  return MuonMass();
18  });
19 
20  //functions for energy in different areas of the detector
21  //---------------------------------
22  // Muon Energy for Active region
23  //---------------------------------
24  float MuonEActive(double trklenact){
25  double p0 = 1.67012e-01;
26  double p1 = 1.79305e-01;
27  double p2 = 3.74708e-03;
28  double p3 = -1.54232e-04;
29  float MuonE = 0.0;
30  if (trklenact <= 0.0) return 0.0;
31  MuonE = p0
32  + p1 * trklenact
33  + p2 * std::pow(trklenact, 2)
34  + p3 * std::pow(trklenact, 3);
35  return MuonE;
36  }
37  //---------------------------------------
38  // Muon Energy for Active+Catcher region
39  //(for Active region Track length)
40  //---------------------------------------
41  float MuonEActiveandCatcher(double trklenactandcat){
42  double p0 = 1.21130e-02;
43  double p1 = 1.97903e-01;
44  double p2 = 7.82459e-04;
45  float MuonE = 0.0;
46  if (trklenactandcat <= 0.0) return 0.0;
47  MuonE = p0
48  + p1 * trklenactandcat
49  + p2 * std::pow(trklenactandcat, 2);
50  return MuonE;
51  }
52  //---------------------------------------
53  // Muon Energy for Active+Catcher region
54  //(for catcher region Track length)
55  //---------------------------------------
56  float MuonECatcher(double trklencat){
57  double offset = 1.31325e-01;
58  double slope = 5.35146e-01;
59  float MuonE = 0.0;
60  if (trklencat <= 0.0) return 0.0;
61  MuonE = slope*trklencat + offset;
62  return MuonE;
63  }
64  //---------------------------------------
65  // Haronic Energy
66  //----------------------------------------
67  float VisibleHadE(float vishadE) {
68  double p0 = 5.85254e-02;
69  double p1 = 1.27796e+00;
70  double p2 = 3.75457e-01;
71  double p3 = -5.45618e-01;
72  double p4 = 1.65975e-01;
73  float HadE = 0.0;
74  HadE = p0
75  + p1 * vishadE
76  + p2 * std::pow(vishadE, 2)
77  + p3 * std::pow(vishadE, 3)
78  + p4 * std::pow(vishadE, 4);
79  return HadE;
80  }
81 
82  //Nu truth vars
83  const NuTruthVar kTrueENu_NT([](const caf::SRNeutrinoProxy* nu)
84  {
85  return nu->E;
86  });
87 
88  const NuTruthVar kTrueEMu_NT([](const caf::SRNeutrinoProxy* nu) //True Muon energy
89  {
90  float MuE = -5.;
91  if(abs(nu->pdg) != 14 || !nu->iscc)
92  return MuE;
93  int nprims = nu->prim.size();
94  for(int iprim = 0; iprim < nprims; iprim++){
95  if(abs(nu->prim[iprim].pdg) == 13){
96  MuE = nu->prim[iprim].p.T();
97  }
98  }//end loop over primaries
99  return MuE;
100  });
101 
102  const NuTruthVar kTrueKEMu_NT([](const caf::SRNeutrinoProxy* nu) //true Muon kinetic energy
103  {
104  float ke = -5;
105  if(abs(nu->pdg) != 14 || !nu->iscc)
106  return ke;
107  int nprims = nu->prim.size();
108  for(int iprim = 0; iprim < nprims; iprim++){
109  if(abs(nu->prim[iprim].pdg) == 13){
110  double E= nu->prim[iprim].p.T();
111  ke = E-MuonMass();
112  }
113  }//end loop over primaries
114  return ke;
115  });
116 
117  //true cosine of muon angle with respect to the beam
119  {
120  if(abs(nu->pdg) != 14 || !nu->iscc)
121  return -5.0;
122  int nprims = nu->prim.size();
123  for(int iprim = 0; iprim < nprims; iprim++){
124  if(abs(nu->prim[iprim].pdg) == 13){
125  TVector3 mudir = nu->prim[iprim].p.Vect();
126  TVector3 beamdir = NuMIBeamDirection(caf::kNEARDET);
127  return mudir.Unit().Dot(beamdir.Unit());
128  }
129  }//end loop over primaries
130  return -5.0;
131  });
132  // True hadronic y
133  const NuTruthVar kTrueY_NT([](const caf::SRNeutrinoProxy* sr)
134  {
135  return float(sr->y);
136  });
137  // True hadronic energy
138  const NuTruthVar kTrueEHad_NT([](const caf::SRNeutrinoProxy* sr)
139  {
140  return float(sr->y * sr->E);
141  });
142 
143  //truth Vars from nu truth Vars
150 
151  //reco cuts
152 
153  //reconstructed available energy
154  //mostly a scale factor on visibile hadronic energy
155  const Var kRecoEAvail([] (const caf::SRProxy *sr)
156  {
157  return 1.68*kNumuHadVisE(sr)+0.02*pow(kNumuHadVisE(sr),2);
158  });
159 
160  const Var kTrkLenActive([](const caf::SRProxy* sr)
161  {
162  int ibesttrk = kBestTrack(sr);
163  if(sr->trk.kalman.ntracks < 1)
164  return -1000.f;
165  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
166  return -1000.f;
167  // If leninact is positive and lenincat is negative,
168  // the transition plane is to the right of the track.
169  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
170  sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
171  return float((sr->trk.kalman.tracks[ibesttrk].leninact / 100.)
172  +(sr->trk.kalman.tracks[ibesttrk].lenincat / 100.));
173  // If leninact is positive and lenincat is positive,
174  // the transition plane is in the middle of the track.
175  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
176  sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
177  return float(sr->trk.kalman.tracks[ibesttrk].leninact / 100.);
178  // If leninact is negative and lenincat is positive,
179  // the transition plane is to the left of the track.
180  if(sr->trk.kalman.tracks[ibesttrk].leninact < 0 &&
181  sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
182  return 0.f;
183  return -1000.f;
184  });
185 
186  const Var kTrkLenCatcher([](const caf::SRProxy* sr)
187  {
188  int ibesttrk = kBestTrack(sr);
189  if(sr->trk.kalman.ntracks < 1)
190  return -1000.f;
191  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
192  return -1000.f;
193  // If leninact is positive and lenincat is negative,
194  // the transition plane is to the right of the track.
195  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
196  sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
197  return 0.f;
198  // If leninact is positive and lenincat is positive,
199  // the transition plane is in the middle of the track.
200  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
201  sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
202  return float(sr->trk.kalman.tracks[ibesttrk].lenincat / 100.);
203  // If leninact is negative and lenincat is positive,
204  // the transition plane is to the left of the track.
205  if(sr->trk.kalman.tracks[ibesttrk].leninact < 0 &&
206  sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
207  return float((sr->trk.kalman.tracks[ibesttrk].leninact / 100.)
208  +(sr->trk.kalman.tracks[ibesttrk].lenincat / 100.));
209  return -1000.f;
210  });
211 
212  // Reconstructed candidate muon angle wrt beam
213  const Var kRecoCosthetaMu([](const caf::SRProxy* sr)
214  {
215  int ibesttrk = kBestTrack(sr);
216  if(sr->trk.kalman.ntracks < 1)
217  return -1000.f;
218  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
219  return -1000.f;
220  TVector3 dir = sr->trk.kalman.tracks[ibesttrk].dir;
221  TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
222  return (float)dir.Dot(beamdir);
223  });
224 
225  // Muon Energy for numuCC 2p2h
226  //------------------------------------
227  const Var k2p2hXSecMuonE([](const caf::SRProxy *sr)
228  {
229  float muonE = 0.0;
230  float muonEact = 0.0;
231  float muonEcat = 0.0;
232  float muonEactandcat = 0.0;
233  float trkLenAct = 0.f;
234  float trkLenCat = 0.f;
235  trkLenAct = kTrkLenActive(sr);
236  trkLenCat = kTrkLenCatcher(sr);
237  int ibesttrk = kBestTrack(sr);
238  if(sr->trk.kalman.ntracks < 1)
239  return -1000.f;
240  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
241  return -1000.f;
242  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
243  sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
244  muonEact = MuonEActive(trkLenAct);
245  else if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
246  sr->trk.kalman.tracks[ibesttrk].lenincat > 0){
247  muonEcat = MuonECatcher(trkLenCat);
248  muonEactandcat = MuonEActiveandCatcher(trkLenAct);
249  muonE = muonEactandcat + muonEcat;
250  }
251  return muonE + muonEact;
252  });
253 
254  // Hadronic Energy for numuCC 2p2h
255  //---------------------------------------
256  const Var k2p2hXSecHadE([](const caf::SRProxy* sr)
257  {
258  int ibesttrk = kBestTrack(sr);
259  if(sr->trk.kalman.ntracks < 1)
260  return -1000.f;
261  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
262  return -1000.f;
263  // Extra energy (hadronic contamination) associated with muon track
264  float ExtraHadE = sr->trk.kalman.tracks[ibesttrk].overlapE;
265  // calorimetric slice energy - Calorimetric Energy of Muon Track
266  float CalhadE = sr->slc.calE - sr->trk.kalman.tracks[ibesttrk].calE;
267  // Add calorimetric hadE and Overlap energy in the muon track
268  float vishadE = CalhadE + ExtraHadE;
269  float hadE = 0.0;
270  hadE = VisibleHadE(vishadE);
271  return hadE;
272  });
273 
274 
276  const Var k2p2hXSecMuKE = k2p2hXSecMuonE - kMuMass;
277 
278  //2d Variables for unfolding
279  //reconstructed available energy vs three momentum transfer
281  q3binsCustom,
282  kRecoEAvail,
284  //true available energy vs three momentum transfer
286  q3binsCustom,
287  kTrueEavail,
289 
290 
291 }
const Var kRecoEAvail([](const caf::SRProxy *sr){return 1.68 *kNumuHadVisE(sr)+0.02 *pow(kNumuHadVisE(sr), 2);})
const Var k2p2hXSecMuKE
Near Detector underground.
Definition: SREnums.h:10
const NuTruthVar kTrueY_NT([](const caf::SRNeutrinoProxy *sr){return float(sr->y);})
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kTrueEMu
const Var kTrueEavail
Definition: TruthVars.h:43
caf::Proxy< size_t > ntracks
Definition: SRProxy.h:1756
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2119
const Var kRecoQmag([](const caf::SRProxy *sr){return sqrt(kRecoQ2(sr)+util::sqr(kHadE(sr)));})
Reconstructed three-momentum transfer.
Definition: NumuVars.h:149
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2108
_Var< T > Var2D(const _Var< T > &a, const Binning &binsa, const _Var< T > &b, const Binning &binsb)
Variable formed from two input variables.
Definition: Var.cxx:247
constexpr T pow(T x)
Definition: pow.h:75
const NuTruthVar kTrueCosthetaMu_NT([](const caf::SRNeutrinoProxy *nu){if(abs(nu->pdg)!=14||!nu->iscc) return-5.0;int nprims=nu->prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(nu->prim[iprim].pdg)==13){TVector3 mudir=nu->prim[iprim].p.Vect();TVector3 beamdir=NuMIBeamDirection(caf::kNEARDET);return mudir.Unit().Dot(beamdir.Unit());}}return-5.0;})
const Binning q3binsCustom
void abs(TH1 *hist)
const Var kRecoEAvailVSq3
const NuTruthVar kTrueKEMu_NT([](const caf::SRNeutrinoProxy *nu){float ke=-5;if(abs(nu->pdg)!=14||!nu->iscc) return ke;int nprims=nu->prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(nu->prim[iprim].pdg)==13){double E=nu->prim[iprim].p.T();ke=E-MuonMass();}}return ke;})
const NuTruthVar kTrueEMu_NT([](const caf::SRNeutrinoProxy *nu){float MuE=-5.;if(abs(nu->pdg)!=14||!nu->iscc) return MuE;int nprims=nu->prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(nu->prim[iprim].pdg)==13){MuE=nu->prim[iprim].p.T();}}return MuE;})
const Var kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
Definition: NumuVars.h:124
const Var kTrkLenCatcher([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f; if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&sr->trk.kalman.tracks[ibesttrk].lenincat< 0) return 0.f; if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&sr->trk.kalman.tracks[ibesttrk].lenincat > 0) return float(sr->trk.kalman.tracks[ibesttrk].lenincat/100.); if(sr->trk.kalman.tracks[ibesttrk].leninact< 0 &&sr->trk.kalman.tracks[ibesttrk].lenincat > 0) return float((sr->trk.kalman.tracks[ibesttrk].leninact/100.) +(sr->trk.kalman.tracks[ibesttrk].lenincat/100.));return-1000.f;})
float MuonEActive(double trklenact)
const Var kTrueEAvailVSq3
const Var k2p2hXSecNuE
const Var kTrueQ3
Definition: TruthVars.h:38
float VisibleHadE(float vishadE)
caf::Proxy< caf::SRTrackBranch > trk
Definition: SRProxy.h:2127
Float_t E
Definition: plot.C:20
const Var kTrueCosthetaMu
double MuonMass()
caf::StandardRecord * sr
const Var kTrueENu
const Var k2p2hXSecHadE([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f;float ExtraHadE=sr->trk.kalman.tracks[ibesttrk].overlapE;float CalhadE=sr->slc.calE-sr->trk.kalman.tracks[ibesttrk].calE;float vishadE=CalhadE+ExtraHadE;float hadE=0.0;hadE=VisibleHadE(vishadE);return hadE;})
const Binning EAvailbinsCustom
float MuonEActiveandCatcher(double trklenactandcat)
const Var kTrueY
caf::Proxy< caf::SRKalman > kalman
Definition: SRProxy.h:1775
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:7
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2124
const Var kRecoCosthetaMu([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f;TVector3 dir=sr->trk.kalman.tracks[ibesttrk].dir;TVector3 beamdir=NuMIBeamDirection(sr->hdr.det);return(float) dir.Dot(beamdir);})
TDirectory * dir
Definition: macro.C:5
caf::Proxy< float > calE
Definition: SRProxy.h:1271
const NuTruthVar kTrueEHad_NT([](const caf::SRNeutrinoProxy *sr){return float(sr->y *sr->E);})
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
const NuTruthVar kTrueENu_NT([](const caf::SRNeutrinoProxy *nu){return nu->E;})
const Var k2p2hXSecMuonE([](const caf::SRProxy *sr){float muonE=0.0;float muonEact=0.0;float muonEcat=0.0;float muonEactandcat=0.0;float trkLenAct=0.f;float trkLenCat=0.f;trkLenAct=kTrkLenActive(sr);trkLenCat=kTrkLenCatcher(sr);int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f;if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&sr->trk.kalman.tracks[ibesttrk].lenincat< 0) muonEact=MuonEActive(trkLenAct);else if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat > 0){muonEcat=MuonECatcher(trkLenCat);muonEactandcat=MuonEActiveandCatcher(trkLenAct);muonE=muonEactandcat+muonEcat;}return muonE+muonEact;})
caf::Proxy< std::vector< caf::SRKalmanTrack > > tracks
Definition: SRProxy.h:1758
const Var kTrkLenActive([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f; if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat< 0) return float((sr->trk.kalman.tracks[ibesttrk].leninact/100.) +(sr->trk.kalman.tracks[ibesttrk].lenincat/100.)); if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat > 0) return float(sr->trk.kalman.tracks[ibesttrk].leninact/100.); if(sr->trk.kalman.tracks[ibesttrk].leninact< 0 && sr->trk.kalman.tracks[ibesttrk].lenincat > 0) return 0.f;return-1000.f;})
float MuonECatcher(double trklencat)
const Var kMuMass([](const caf::SRProxy *sr){return MuonMass();})
const Var kTrueEHad
const Var kBestTrack
Definition: NDXSecMuonPID.h:33
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:231
const Var kTrueKEMu