Changes between Version 44 and Version 45 of HowToAddSample

Timestamp:

Jan 6, 2017, 2:51:37 PM (7 years ago)

Author:

Clarence Wret

Comment:

--

HowToAddSample

@@ -255 +255 @@
 The `FitEvent` class provides numerous getter functions to make `FillEventVariables` as simple as possible.
 
-=== For the p,,mu,, ===
+[=#point_muon_dist]
+=== For the p,,mu,, distribution ===
 
 To get the muon momentum, we first need a muon in the event. The function `FitEvent::NumFSParticle(int pdg)` counts the number of final-state particles with pdg-code `pdg` and returns that number. So our first check is 
@@ -284 +285 @@
 
 {{{
-fXVar = p_mu
+fXVar = p_mu;
 
 return;
 }}}
-=== For the E^rec^,,nu,, ===
-
-
+
+which concludes the `FillEventVariables(FitEvent *event)` function.
+
+=== For the E^rec^,,nu,, distribution ===
+
+E^rec^,,nu,, needs to be reconstructed from the outgoing particles, so requires a little more work. However it follows very similar steps to the [#point_muon_dist muon case].
+
+To reconstructed E^rec^,,nu,, we need:
+
+ * Muon momentum
+ * Pion momentum
+ * Muon/neutrino angle
+ * Pion/neutrino angle
+ * Muon/pion angle
+ * All the particle masses
+
+So we proceed to find the `TLorentzVector`s in the event as before. We first need to make sure we have a positive pion and a negative muon in the event:
+
+{{{
+// Need to make sure there's a muon
+if (event->NumFSParticle(13) == 0) return;
+// Need to make sure there's a pion
+if (event->NumFSParticle(211) == 0) return;
+}}}
+
+Then get the relevant `TLorentzVector`s:
+
+{{{
+// Get the incoming neutrino
+TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
+// Get the muon
+TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
+// Get the pion
+TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
+}}}
+
+We should now write a helper function in the `FitUtils` namespace to get `E^rec^,,nu,,` from the muon, pion and neutrino `TLorentzVector`s. The `FitUtils` implementation lives in `src/Utils/FitUtils.cxx`.
+
+{{{
+double FitUtils::EnuCC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi) {
+
+  // The muon energy, momentum and mass
+  double E_mu = pmu.E() / 1000.;
+  double p_mu = pmu.Vect().Mag() / 1000.;
+  double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
+
+  // The pion energy, momentum and mass
+  double E_pi = ppi.E() / 1000.;
+  double p_pi = ppi.Vect().Mag() / 1000.;                                        
+  double m_pi = sqrt(E_pi * E_pi - p_pi * p_pi);
+
+  // The relevant angles
+  double th_nu_pi = pnu.Vect().Angle(ppi.Vect());
+  double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
+  double th_pi_mu = ppi.Vect().Angle(pmu.Vect());
+
+  // Now write down the derived neutrino energy, assuming the initial state nucleon is at rest and we only have a pion, muon and nucleon in the final state
+  double rEnu = (m_mu * m_mu + m_pi * m_pi - 2 * m_n * (E_pi + E_mu) + 2 * E_pi * E_mu -
+       2 * p_pi * p_mu * cos(th_pi_mu)) /
+      (2 * (E_pi + E_mu - p_pi * cos(th_nu_pi) - p_mu * cos(th_nu_mu) - m_n));
+
+  return rEnu;
+};
+}}}
+
+Once we have the helper function implemented and tested we go back to the `FillEventVariables(FitEvent *event)` implementation and simply do
+
+{{{
+double Enu = FitUtils::EnuCC1piprec(Pnu, Pmu, Ppip);
+
+fXVar = Enu;
+
+return;
+}}}
+
+and `FillEventVariables(FitEvent *event)` is done!
 
 [=#point_signaldef]