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