GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
NuEPrimaryLeptonGenerator.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <TMath.h>
12 #include <TVector3.h>
13 
23 
24 using namespace genie;
25 using namespace genie::constants;
26 
27 //___________________________________________________________________________
29 PrimaryLeptonGenerator("genie::NuEPrimaryLeptonGenerator")
30 {
31 
32 }
33 //___________________________________________________________________________
35 PrimaryLeptonGenerator("genie::NuEPrimaryLeptonGenerator", config)
36 {
37 
38 }
39 //___________________________________________________________________________
41 {
42 
43 }
44 //___________________________________________________________________________
46 {
47 // This method generates the final state primary lepton for NuE events
48 
49  Interaction * interaction = evrec->Summary();
50  const InitialState & init_state = interaction->InitState();
51 
52  // Get selected kinematics
53  double y = interaction->Kine().y(true);
54  assert(y>0 && y<1);
55 
56  // Final state primary lepton PDG code
57  int pdgc = interaction->FSPrimLeptonPdg();
58  assert(pdgc!=0);
59 
60  // Compute the neutrino and muon energy
61  double Ev = init_state.ProbeE(kRfLab);
62  double El = (1-y)*Ev;
63 
64  LOG("LeptonicVertex", pINFO)
65  << "Ev = " << Ev << ", y = " << y << ", -> El = " << El;
66 
67  // Compute the momentum transfer and scattering angle
68  double El2 = TMath::Power(El,2);
69  double me = kElectronMass;
70  double ml = interaction->FSPrimLepton()->Mass();
71  double ml2 = TMath::Power(ml,2);
72  double pl = TMath::Sqrt(El2-ml2);
73 
74  assert(El2>=ml2);
75 
76  double Q2 = 2*(Ev-El)*me + me*me;
77  double costh = (El-0.5*(Q2+ml2)/Ev)/pl;
78  double sinth = TMath::Sqrt( TMath::Max(0., 1-TMath::Power(costh,2.)) );
79 
80  LOG("LeptonicVertex", pNOTICE)
81  << "Q2 = " << Q2 << ", cos(theta) = " << costh;
82 
83  //warn about overflow in costheta and ignore it if it is small or abort
84  if( TMath::Abs(costh)>1 ) {
85  LOG("LeptonicVertex", pWARN)
86  << "El = " << El << ", Ev = " << Ev << ", cos(theta) = " << costh;
87  //if(TMath::Abs(costh)-1<0.3) costh = 1.0; //why?
88  }
89  assert(TMath::Abs(costh)<=1);
90 
91  // Compute the p components along and perpendicular the v direction
92  double plp = pl * costh; // p(//)
93  double plt = pl * sinth; // p(-|)
94 
95  LOG("LeptonicVertex", pNOTICE)
96  << "fsl: E = " << El << ", |p//| = " << plp << "[pT] = " << plt;
97 
98  // Randomize transverse components
100  double phi = 2*kPi * rnd->RndLep().Rndm();
101  double pltx = plt * TMath::Cos(phi);
102  double plty = plt * TMath::Sin(phi);
103 
104  // Take a unit vector along the neutrino direction
105  TVector3 unit_nudir = evrec->Probe()->P4()->Vect().Unit();
106 
107  // Rotate lepton momentum vector from the reference frame (x'y'z') where
108  // {z':(neutrino direction), z'x':(theta plane)} to the LAB
109  TVector3 p3l(pltx,plty,plp);
110  p3l.RotateUz(unit_nudir);
111 
112  // Lepton 4-momentum in the LAB
113  TLorentzVector p4l(p3l,El);
114 
115  // Create a GHepParticle and add it to the event record
116  this->AddToEventRecord(evrec, pdgc, p4l);
117 
118  // Set final state lepton polarization
119  this->SetPolarization(evrec);
120 }
121 //___________________________________________________________________________
TRandom3 & RndLep(void) const
rnd number generator used by final state primary lepton generators
Definition: RandomGen.h:62
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:91
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
static const double kElectronMass
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
int FSPrimLeptonPdg(void) const
final state primary lepton pdg
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:284
double y(bool selected=false) const
Definition: Kinematics.cxx:112
Summary information for an interaction.
Definition: Interaction.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
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:62
#define pWARN
Definition: Messenger.h:60
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
void ProcessEventRecord(GHepRecord *event_rec) const
virtual void SetPolarization(GHepRecord *ev) const
const InitialState & InitState(void) const
Definition: Interaction.h:69
#define pNOTICE
Definition: Messenger.h:61
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:45
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Initial State information.
Definition: InitialState.h:48