nuisance is hosted by Hepforge, IPPP Durham
close Warning: Can't synchronize with repository "(default)" ("(default)" is not readable or not a Git repository.). Look in the Trac log for more information.

Changes between Version 44 and Version 45 of HowToAddSample


Ignore:
Timestamp:
Jan 6, 2017, 2:51:37 PM (8 years ago)
Author:
Clarence Wret
Comment:

--

Legend:

Unmodified
Added
Removed
Modified
  • HowToAddSample

    v44 v45  
    255255The `FitEvent` class provides numerous getter functions to make `FillEventVariables` as simple as possible.
    256256
    257 === For the p,,mu,, ===
     257[=#point_muon_dist]
     258=== For the p,,mu,, distribution ===
    258259
    259260To 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
     
    284285
    285286{{{
    286 fXVar = p_mu
     287fXVar = p_mu;
    287288
    288289return;
    289290}}}
    290 === For the E^rec^,,nu,, ===
    291 
    292 
     291
     292which concludes the `FillEventVariables(FitEvent *event)` function.
     293
     294=== For the E^rec^,,nu,, distribution ===
     295
     296E^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].
     297
     298To reconstructed E^rec^,,nu,, we need:
     299
     300 * Muon momentum
     301 * Pion momentum
     302 * Muon/neutrino angle
     303 * Pion/neutrino angle
     304 * Muon/pion angle
     305 * All the particle masses
     306
     307So 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:
     308
     309{{{
     310// Need to make sure there's a muon
     311if (event->NumFSParticle(13) == 0) return;
     312// Need to make sure there's a pion
     313if (event->NumFSParticle(211) == 0) return;
     314}}}
     315
     316Then get the relevant `TLorentzVector`s:
     317
     318{{{
     319// Get the incoming neutrino
     320TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
     321// Get the muon
     322TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
     323// Get the pion
     324TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
     325}}}
     326
     327We 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`.
     328
     329{{{
     330double FitUtils::EnuCC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi) {
     331
     332  // The muon energy, momentum and mass
     333  double E_mu = pmu.E() / 1000.;
     334  double p_mu = pmu.Vect().Mag() / 1000.;
     335  double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
     336
     337  // The pion energy, momentum and mass
     338  double E_pi = ppi.E() / 1000.;
     339  double p_pi = ppi.Vect().Mag() / 1000.;                                       
     340  double m_pi = sqrt(E_pi * E_pi - p_pi * p_pi);
     341
     342  // The relevant angles
     343  double th_nu_pi = pnu.Vect().Angle(ppi.Vect());
     344  double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
     345  double th_pi_mu = ppi.Vect().Angle(pmu.Vect());
     346
     347  // 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
     348  double rEnu = (m_mu * m_mu + m_pi * m_pi - 2 * m_n * (E_pi + E_mu) + 2 * E_pi * E_mu -
     349       2 * p_pi * p_mu * cos(th_pi_mu)) /
     350      (2 * (E_pi + E_mu - p_pi * cos(th_nu_pi) - p_mu * cos(th_nu_mu) - m_n));
     351
     352  return rEnu;
     353};
     354}}}
     355
     356Once we have the helper function implemented and tested we go back to the `FillEventVariables(FitEvent *event)` implementation and simply do
     357
     358{{{
     359double Enu = FitUtils::EnuCC1piprec(Pnu, Pmu, Ppip);
     360
     361fXVar = Enu;
     362
     363return;
     364}}}
     365
     366and `FillEventVariables(FitEvent *event)` is done!
    293367
    294368[=#point_signaldef]