GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PrimaryLeptonGenerator.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 <TVector3.h>
12 #include <TLorentzVector.h>
13 
29 
30 using namespace genie;
31 using namespace genie::constants;
32 
33 //___________________________________________________________________________
36 {
37 
38 }
39 //___________________________________________________________________________
42 {
43 
44 }
45 //___________________________________________________________________________
47 EventRecordVisitorI(name, config)
48 {
49 
50 }
51 //___________________________________________________________________________
53 {
54 
55 }
56 //___________________________________________________________________________
58 {
59 // This method generates the final state primary lepton
60 
61  Interaction * interaction = evrec->Summary();
62 
63  // Boost vector for [LAB] <-> [Nucleon Rest Frame] transforms
64  TVector3 beta = this->NucRestFrame2Lab(evrec);
65 
66  // Neutrino 4p
67  TLorentzVector * p4v = evrec->Probe()->GetP4(); // v 4p @ LAB
68  p4v->Boost(-1.*beta); // v 4p @ Nucleon rest frame
69 
70  // Look-up selected kinematics & other needed kinematical params
71  double Q2 = interaction->Kine().Q2(true);
72  double y = interaction->Kine().y(true);
73  double Ev = p4v->E();
74  double ml = interaction->FSPrimLepton()->Mass();
75  double ml2 = TMath::Power(ml,2);
76 
77  LOG("LeptonicVertex", pNOTICE)
78  << "Ev = " << Ev << ", Q2 = " << Q2 << ", y = " << y;
79 
80  // Compute the final state primary lepton energy and momentum components
81  // along and perpendicular the neutrino direction
82  double El = (1-y)*Ev;
83  double plp = El - 0.5*(Q2+ml2)/Ev; // p(//)
84  double plt = TMath::Sqrt(TMath::Max(0.,El*El-plp*plp-ml2)); // p(-|)
85 
86  LOG("LeptonicVertex", pNOTICE)
87  << "fsl: E = " << El << ", |p//| = " << plp << ", [pT] = " << plt;
88 
89  // Randomize transverse components
91  double phi = 2*kPi * rnd->RndLep().Rndm();
92  double pltx = plt * TMath::Cos(phi);
93  double plty = plt * TMath::Sin(phi);
94 
95  // Take a unit vector along the neutrino direction @ the nucleon rest frame
96  TVector3 unit_nudir = p4v->Vect().Unit();
97 
98  // Rotate lepton momentum vector from the reference frame (x'y'z') where
99  // {z':(neutrino direction), z'x':(theta plane)} to the nucleon rest frame
100  TVector3 p3l(pltx,plty,plp);
101  p3l.RotateUz(unit_nudir);
102 
103  // Lepton 4-momentum in the nucleon rest frame
104  TLorentzVector p4l(p3l,El);
105 
106  LOG("LeptonicVertex", pNOTICE)
107  << "fsl @ NRF: " << utils::print::P4AsString(&p4l);
108 
109  // Boost final state primary lepton to the lab frame
110  p4l.Boost(beta); // active Lorentz transform
111 
112  LOG("LeptonicVertex", pNOTICE)
113  << "fsl @ LAB: " << utils::print::P4AsString(&p4l);
114 
115  // Figure out the Final State Lepton PDG Code
116  int pdgc = interaction->FSPrimLepton()->PdgCode();
117 
118  // Create a GHepParticle and add it to the event record
119  this->AddToEventRecord(evrec, pdgc, p4l);
120 
121  // Set final state lepton polarization
122  this->SetPolarization(evrec);
123 
124  delete p4v;
125 }
126 //___________________________________________________________________________
128 {
129 // Velocity for an active Lorentz transform taking the final state primary
130 // lepton from the [nucleon rest frame] --> [LAB]
131 
132  Interaction * interaction = evrec->Summary();
133  const InitialState & init_state = interaction->InitState();
134 
135  const TLorentzVector & pnuc4 = init_state.Tgt().HitNucP4(); //[@LAB]
136  TVector3 beta = pnuc4.BoostVector();
137 
138  return beta;
139 }
140 //___________________________________________________________________________
142  GHepRecord * evrec, int pdgc, const TLorentzVector & p4) const
143 {
144 // Adds the final state primary lepton GHepParticle to the event record.
145 // To be called by all concrete PrimaryLeptonGenerators before exiting.
146 
147  Interaction * interaction = evrec->Summary();
148 
149  GHepParticle * mom = evrec->Probe();
150  int imom = evrec->ProbePosition();
151 
152  const TLorentzVector & vtx = *(mom->X4());
153 
154  TLorentzVector x4l(vtx); // position 4-vector
155  TLorentzVector p4l(p4); // momentum 4-vector
156 
157  GHepParticle * nucltgt = evrec->TargetNucleus();
158 
159  bool is_ve = interaction->ProcInfo().IsInverseMuDecay() ||
160  interaction->ProcInfo().IsIMDAnnihilation() ||
161  interaction->ProcInfo().IsNuElectronElastic();
162 
163  bool can_correct = fApplyCoulombCorrection && nucltgt!=0 && !is_ve;
164  if(can_correct) {
165  LOG("LeptonicVertex", pINFO)
166  << "Correcting f/s lepton energy for Coulomb effects";
167 
168  double m = interaction->FSPrimLepton()->Mass();
169  double Z = nucltgt->Z();
170  double A = nucltgt->A();
171 
172  // charge radius of nucleus in GeV^-1
173  double Rc = (1.1*TMath::Power(A,1./3.) + 0.86*TMath::Power(A,-1./3.))/0.197;
174 
175  // shift of lepton energy in homogenous sphere with radius Rc
176  double Vo = 3*kAem*Z/(2*Rc);
177  Vo *= 0.75; // as suggested in R.Gran's note
178 
179  double Elo = p4l.Energy();
180  double e = TMath::Min(Vo, Elo-m);
181  double El = TMath::Max(0., Elo-e);
182 
183  LOG("LeptonicVertex", pINFO)
184  << "Lepton energy subtraction: E = " << Elo << " --> " << El;
185 
186  double pmag_old = p4l.P();
187  double pmag_new = TMath::Sqrt(utils::math::NonNegative(El*El-m*m));
188  double scale = pmag_new / pmag_old;
189  LOG("LeptonicVertex", pDEBUG)
190  << "|pnew| = " << pmag_new << ", |pold| = " << pmag_old
191  << ", scale = " << scale;
192 
193  double pxl = scale * p4l.Px();
194  double pyl = scale * p4l.Py();
195  double pzl = scale * p4l.Pz();
196 
197  p4l.SetPxPyPzE(pxl,pyl,pzl,El);
198 
199  TLorentzVector p4c = p4 - p4l;
200  TLorentzVector x4dummy(0,0,0,0);;
201 
202  evrec->AddParticle(
203  kPdgCoulobtron, kIStStableFinalState, -1,-1,-1,-1, p4c, x4dummy);
204  }
205 
206  evrec->AddParticle(pdgc, kIStStableFinalState, imom,-1,-1,-1, p4l, x4l);
207 
208  // update the interaction summary
209  evrec->Summary()->KinePtr()->SetFSLeptonP4(p4l);
210 }
211 //___________________________________________________________________________
213 {
214  // The default treatment is just to delegate to the utility function.
215  // Derived classes that need a different approach may override this
216  // member function.
218 }
219 //___________________________________________________________________________
221 {
222  Algorithm::Configure(config);
223  this->LoadConfig();
224 }
225 //____________________________________________________________________________
227 {
228  Algorithm::Configure(config);
229  this->LoadConfig();
230 }
231 //____________________________________________________________________________
233 {
234  GetParam( "ApplyCoulombCorrection", fApplyCoulombCorrection ) ;
235 
236 }
237 //___________________________________________________________________________
int Z(void) const
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
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
bool IsInverseMuDecay(void) const
const int kPdgCoulobtron
Definition: PDGCodes.h:213
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:27
bool IsIMDAnnihilation(void) const
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
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
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:91
const double e
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void SetFSLeptonP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:297
bool IsNuElectronElastic(void) const
virtual void ProcessEventRecord(GHepRecord *evrec) const
static constexpr double A
Definition: Units.h:74
TLorentzVector * GetP4(void) const
const Kinematics & Kine(void) const
Definition: Interaction.h:71
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
#define pINFO
Definition: Messenger.h:62
virtual TVector3 NucRestFrame2Lab(GHepRecord *ev) const
void Configure(const Registry &config)
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
virtual void SetPolarization(GHepRecord *ev) const
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double NonNegative(double x)
Definition: MathUtils.cxx:273
int A(void) const
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
#define pNOTICE
Definition: Messenger.h:61
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
virtual int ProbePosition(void) const
Definition: GHepRecord.cxx:352
void SetPrimaryLeptonPolarization(GHepRecord *ev)
virtual void AddToEventRecord(GHepRecord *ev, int pdgc, const TLorentzVector &p4) const
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
static constexpr double m
Definition: Units.h:71
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63