NuEPrimaryLeptonGenerator.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Feb 09, 2009 - CA
14  Moved into the NuE package from its previous location (EVGModules package)
15  @ Feb 12, 2013 - CA (code from Rosen Matev)
16  Add mass_electron^2 term in kinematical calculation.
17 */
18 //____________________________________________________________________________
19 
20 #include <TMath.h>
21 #include <TVector3.h>
22 
32 
33 using namespace genie;
34 using namespace genie::constants;
35 
36 //___________________________________________________________________________
38 PrimaryLeptonGenerator("genie::NuEPrimaryLeptonGenerator")
39 {
40 
41 }
42 //___________________________________________________________________________
44 PrimaryLeptonGenerator("genie::NuEPrimaryLeptonGenerator", config)
45 {
46 
47 }
48 //___________________________________________________________________________
50 {
51 
52 }
53 //___________________________________________________________________________
55 {
56 // This method generates the final state primary lepton for NuE events
57 
58  Interaction * interaction = evrec->Summary();
59  const InitialState & init_state = interaction->InitState();
60 
61  // Get selected kinematics
62  double y = interaction->Kine().y(true);
63  assert(y>0 && y<1);
64 
65  // Final state primary lepton PDG code
66  int pdgc = interaction->FSPrimLeptonPdg();
67  assert(pdgc!=0);
68 
69  // Compute the neutrino and muon energy
70  double Ev = init_state.ProbeE(kRfLab);
71  double El = (1-y)*Ev;
72 
73  LOG("LeptonicVertex", pINFO)
74  << "Ev = " << Ev << ", y = " << y << ", -> El = " << El;
75 
76  // Compute the momentum transfer and scattering angle
77  double El2 = TMath::Power(El,2);
78  double me = kElectronMass;
79  double ml = interaction->FSPrimLepton()->Mass();
80  double ml2 = TMath::Power(ml,2);
81  double pl = TMath::Sqrt(El2-ml2);
82 
83  assert(El2>=ml2);
84 
85  double Q2 = 2*(Ev-El)*me + me*me;
86  double costh = (El-0.5*(Q2+ml2)/Ev)/pl;
87  double sinth = TMath::Sqrt( TMath::Max(0., 1-TMath::Power(costh,2.)) );
88 
89  LOG("LeptonicVertex", pNOTICE)
90  << "Q2 = " << Q2 << ", cos(theta) = " << costh;
91 
92  //warn about overflow in costheta and ignore it if it is small or abort
93  if( TMath::Abs(costh)>1 ) {
94  LOG("LeptonicVertex", pWARN)
95  << "El = " << El << ", Ev = " << Ev << ", cos(theta) = " << costh;
96  //if(TMath::Abs(costh)-1<0.3) costh = 1.0; //why?
97  }
98  assert(TMath::Abs(costh)<=1);
99 
100  // Compute the p components along and perpendicular the v direction
101  double plp = pl * costh; // p(//)
102  double plt = pl * sinth; // p(-|)
103 
104  LOG("LeptonicVertex", pNOTICE)
105  << "fsl: E = " << El << ", |p//| = " << plp << "[pT] = " << plt;
106 
107  // Randomize transverse components
109  double phi = 2*kPi * rnd->RndLep().Rndm();
110  double pltx = plt * TMath::Cos(phi);
111  double plty = plt * TMath::Sin(phi);
112 
113  // Take a unit vector along the neutrino direction
114  TVector3 unit_nudir = evrec->Probe()->P4()->Vect().Unit();
115 
116  // Rotate lepton momentum vector from the reference frame (x'y'z') where
117  // {z':(neutrino direction), z'x':(theta plane)} to the LAB
118  TVector3 p3l(pltx,plty,plp);
119  p3l.RotateUz(unit_nudir);
120 
121  // Lepton 4-momentum in the LAB
122  TLorentzVector p4l(p3l,El);
123 
124  // Create a GHepParticle and add it to the event record
125  this->AddToEventRecord(evrec, pdgc, p4l);
126 
127  // Set final state lepton polarization
128  this->SetPolarization(evrec);
129 }
130 //___________________________________________________________________________
const double kPi
Basic constants.
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
TRandom3 & RndLep(void) const
rnd number generator used by final state primary lepton generators
Definition: RandomGen.h:63
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:135
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
Definition: config.py:1
static const double kElectronMass
Definition: Constants.h:71
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
int FSPrimLeptonPdg(void) const
final state primary lepton pdg
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
double y(bool selected=false) const
Definition: Kinematics.cxx:122
Summary information for an interaction.
Definition: Interaction.h:56
Definition: NueSkimmer.h:24
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const Kinematics & Kine(void) const
Definition: Interaction.h:71
Abstract class. Is used to pass common implementation to concrete implementations of the EventRecordV...
#define pINFO
Definition: Messenger.h:63
#define pWARN
Definition: Messenger.h:61
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
void ProcessEventRecord(GHepRecord *event_rec) const
virtual void SetPolarization(GHepRecord *ev) const
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
assert(nhit_max >=nhit_nbins)
const InitialState & InitState(void) const
Definition: Interaction.h:69
#define pNOTICE
Definition: Messenger.h:62
virtual void AddToEventRecord(GHepRecord *ev, int pdgc, const TLorentzVector &p4) const
double ProbeE(RefFrame_t rf) const
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Initial State information.
Definition: InitialState.h:49