PrimaryLeptonGenerator.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 12, 2013 - CA (code from Rosen Matev)
14  Handle the IMD annihilation channel.
15 
16 */
17 //____________________________________________________________________________
18 
19 #include <TVector3.h>
20 #include <TLorentzVector.h>
21 
36 
37 using namespace genie;
38 using namespace genie::constants;
39 
40 //___________________________________________________________________________
43 {
44 
45 }
46 //___________________________________________________________________________
49 {
50 
51 }
52 //___________________________________________________________________________
54 EventRecordVisitorI(name, config)
55 {
56 
57 }
58 //___________________________________________________________________________
60 {
61 
62 }
63 //___________________________________________________________________________
65 {
66 // This method generates the final state primary lepton
67 
68  Interaction * interaction = evrec->Summary();
69 
70  // Boost vector for [LAB] <-> [Nucleon Rest Frame] transforms
71  TVector3 beta = this->NucRestFrame2Lab(evrec);
72 
73  // Neutrino 4p
74  TLorentzVector * p4v = evrec->Probe()->GetP4(); // v 4p @ LAB
75  p4v->Boost(-1.*beta); // v 4p @ Nucleon rest frame
76 
77  // Look-up selected kinematics & other needed kinematical params
78  double Q2 = interaction->Kine().Q2(true);
79  double y = interaction->Kine().y(true);
80  double Ev = p4v->E();
81  double ml = interaction->FSPrimLepton()->Mass();
82  double ml2 = TMath::Power(ml,2);
83 
84  LOG("LeptonicVertex", pNOTICE)
85  << "Ev = " << Ev << ", Q2 = " << Q2 << ", y = " << y;
86 
87  // Compute the final state primary lepton energy and momentum components
88  // along and perpendicular the neutrino direction
89  double El = (1-y)*Ev;
90  double plp = El - 0.5*(Q2+ml2)/Ev; // p(//)
91  double plt = TMath::Sqrt(TMath::Max(0.,El*El-plp*plp-ml2)); // p(-|)
92 
93  LOG("LeptonicVertex", pNOTICE)
94  << "fsl: E = " << El << ", |p//| = " << plp << ", [pT] = " << plt;
95 
96  // Randomize transverse components
98  double phi = 2*kPi * rnd->RndLep().Rndm();
99  double pltx = plt * TMath::Cos(phi);
100  double plty = plt * TMath::Sin(phi);
101 
102  // Take a unit vector along the neutrino direction @ the nucleon rest frame
103  TVector3 unit_nudir = p4v->Vect().Unit();
104 
105  // Rotate lepton momentum vector from the reference frame (x'y'z') where
106  // {z':(neutrino direction), z'x':(theta plane)} to the nucleon rest frame
107  TVector3 p3l(pltx,plty,plp);
108  p3l.RotateUz(unit_nudir);
109 
110  // Lepton 4-momentum in the nucleon rest frame
111  TLorentzVector p4l(p3l,El);
112 
113  LOG("LeptonicVertex", pNOTICE)
114  << "fsl @ NRF: " << utils::print::P4AsString(&p4l);
115 
116  // Boost final state primary lepton to the lab frame
117  p4l.Boost(beta); // active Lorentz transform
118 
119  LOG("LeptonicVertex", pNOTICE)
120  << "fsl @ LAB: " << utils::print::P4AsString(&p4l);
121 
122  // Figure out the Final State Lepton PDG Code
123  int pdgc = interaction->FSPrimLepton()->PdgCode();
124 
125  // Create a GHepParticle and add it to the event record
126  this->AddToEventRecord(evrec, pdgc, p4l);
127 
128  // Set final state lepton polarization
129  this->SetPolarization(evrec);
130 
131  delete p4v;
132 }
133 //___________________________________________________________________________
135 {
136 // Velocity for an active Lorentz transform taking the final state primary
137 // lepton from the [nucleon rest frame] --> [LAB]
138 
139  Interaction * interaction = evrec->Summary();
140  const InitialState & init_state = interaction->InitState();
141 
142  const TLorentzVector & pnuc4 = init_state.Tgt().HitNucP4(); //[@LAB]
143  TVector3 beta = pnuc4.BoostVector();
144 
145  return beta;
146 }
147 //___________________________________________________________________________
149  GHepRecord * evrec, int pdgc, const TLorentzVector & p4) const
150 {
151 // Adds the final state primary lepton GHepParticle to the event record.
152 // To be called by all concrete PrimaryLeptonGenerators before exiting.
153 
154  Interaction * interaction = evrec->Summary();
155 
156  GHepParticle * mom = evrec->Probe();
157  int imom = evrec->ProbePosition();
158 
159  const TLorentzVector & vtx = *(mom->X4());
160 
161  TLorentzVector x4l(vtx); // position 4-vector
162  TLorentzVector p4l(p4); // momentum 4-vector
163 
164  GHepParticle * nucltgt = evrec->TargetNucleus();
165 
166  bool is_ve = interaction->ProcInfo().IsInverseMuDecay() ||
167  interaction->ProcInfo().IsIMDAnnihilation() ||
168  interaction->ProcInfo().IsNuElectronElastic();
169 
170  bool can_correct = fApplyCoulombCorrection && nucltgt!=0 && !is_ve;
171  if(can_correct) {
172  LOG("LeptonicVertex", pINFO)
173  << "Correcting f/s lepton energy for Coulomb effects";
174 
175  double m = interaction->FSPrimLepton()->Mass();
176  double Z = nucltgt->Z();
177  double A = nucltgt->A();
178 
179  // charge radius of nucleus in GeV^-1
180  double Rc = (1.1*TMath::Power(A,1./3.) + 0.86*TMath::Power(A,-1./3.))/0.197;
181 
182  // shift of lepton energy in homogenous sphere with radius Rc
183  double Vo = 3*kAem*Z/(2*Rc);
184  Vo *= 0.75; // as suggested in R.Gran's note
185 
186  double Elo = p4l.Energy();
187  double e = TMath::Min(Vo, Elo-m);
188  double El = TMath::Max(0., Elo-e);
189 
190  LOG("LeptonicVertex", pINFO)
191  << "Lepton energy subtraction: E = " << Elo << " --> " << El;
192 
193  double pmag_old = p4l.P();
194  double pmag_new = TMath::Sqrt(utils::math::NonNegative(El*El-m*m));
195  double scale = pmag_new / pmag_old;
196  LOG("LeptonicVertex", pDEBUG)
197  << "|pnew| = " << pmag_new << ", |pold| = " << pmag_old
198  << ", scale = " << scale;
199 
200  double pxl = scale * p4l.Px();
201  double pyl = scale * p4l.Py();
202  double pzl = scale * p4l.Pz();
203 
204  p4l.SetPxPyPzE(pxl,pyl,pzl,El);
205 
206  TLorentzVector p4c = p4 - p4l;
207  TLorentzVector x4dummy(0,0,0,0);;
208 
209  evrec->AddParticle(
210  kPdgCoulobtron, kIStStableFinalState, -1,-1,-1,-1, p4c, x4dummy);
211  }
212 
213  evrec->AddParticle(pdgc, kIStStableFinalState, imom,-1,-1,-1, p4l, x4l);
214 
215  // update the interaction summary
216  evrec->Summary()->KinePtr()->SetFSLeptonP4(p4l);
217 }
218 //___________________________________________________________________________
220 {
221 // Set the final state lepton polarization. A mass-less lepton would be fully
222 // polarized. This would be exact for neutrinos and a very good approximation
223 // for electrons for the energies this generator is going to be used. This is
224 // not the case for muons and, mainly, for taus. I need to refine this later.
225 // How? See Kuzmin, Lyubushkin and Naumov, hep-ph/0312107
226 
227  // get the final state primary lepton
229  if(!fsl) {
230  LOG("LeptonicVertex", pERROR)
231  << "Final state lepton not set yet! \n" << *ev;
232  return;
233  }
234 
235  // get (px,py,pz) @ LAB
236  TVector3 plab(fsl->Px(), fsl->Py(), fsl->Pz());
237 
238  // in the limit m/E->0: leptons are left-handed and their anti-particles
239  // are right-handed
240  int pdgc = fsl->Pdg();
241  if(pdg::IsNeutrino(pdgc) || pdg::IsElectron(pdgc) ||
242  pdg::IsMuon(pdgc) || pdg::IsTau(pdgc) ) {
243  plab *= -1; // left-handed
244  }
245 
246  LOG("LeptonicVertex", pINFO)
247  << "Setting polarization angles for particle: " << fsl->Name();
248 
249  fsl->SetPolarization(plab);
250 
251  if(fsl->PolzIsSet()) {
252  LOG("LeptonicVertex", pINFO)
253  << "Polarization (rad): Polar = " << fsl->PolzPolarAngle()
254  << ", Azimuthal = " << fsl->PolzAzimuthAngle();
255  }
256 }
257 //___________________________________________________________________________
259 {
260  Algorithm::Configure(config);
261  this->LoadConfig();
262 }
263 //____________________________________________________________________________
265 {
266  Algorithm::Configure(config);
267  this->LoadConfig();
268 }
269 //____________________________________________________________________________
271 {
272  GetParam( "ApplyCoulombCorrection", fApplyCoulombCorrection ) ;
273 
274 }
275 //___________________________________________________________________________
276 
const double kPi
const XML_Char * name
Definition: expat.h:151
int Z(void) const
Basic constants.
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
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
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
double PolzPolarAngle(void) const
Definition: GHepParticle.h:120
bool IsInverseMuDecay(void) const
const int kPdgCoulobtron
Definition: PDGCodes.h:190
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:34
Definition: config.py:1
bool IsIMDAnnihilation(void) const
Double_t beta
void SetPolarization(double theta, double phi)
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:91
double Px(void) const
Get Px.
Definition: GHepParticle.h:89
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
double y(bool selected=false) const
Definition: Kinematics.cxx:122
static const double kAem
Definition: Constants.h:57
Float_t Z
Definition: plot.C:38
int Pdg(void) const
Definition: GHepParticle.h:64
string Name(void) const
Name that corresponds to the PDG code.
Summary information for an interaction.
Definition: Interaction.h:56
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:92
Double_t scale
Definition: plot.C:25
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void SetFSLeptonP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:307
bool IsNuElectronElastic(void) const
virtual void ProcessEventRecord(GHepRecord *evrec) const
TLorentzVector * GetP4(void) const
virtual GHepParticle * FinalStatePrimaryLepton(void) const
Definition: GHepRecord.cxx:370
const Kinematics & Kine(void) const
Definition: Interaction.h:71
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
bool IsTau(int pdgc)
Definition: PDGUtils.cxx:200
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
#define pINFO
Definition: Messenger.h:63
bool IsMuon(int pdgc)
Definition: PDGUtils.cxx:190
virtual TVector3 NucRestFrame2Lab(GHepRecord *ev) const
bool PolzIsSet(void) const
void Configure(const Registry &config)
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
static const double A
Definition: Units.h:82
virtual void SetPolarization(GHepRecord *ev) const
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:535
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:280
int A(void) const
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
#define pNOTICE
Definition: Messenger.h:62
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:67
virtual int ProbePosition(void) const
Definition: GHepRecord.cxx:389
double PolzAzimuthAngle(void) const
Definition: GHepParticle.h:121
Float_t e
Definition: plot.C:35
virtual void AddToEventRecord(GHepRecord *ev, int pdgc, const TLorentzVector &p4) 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...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
bool IsElectron(int pdgc)
Definition: PDGUtils.cxx:180
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
double Py(void) const
Get Py.
Definition: GHepParticle.h:90