OutgoingDarkGenerator.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: Joshua Berger <jberger \at physics.wisc.edu>
8  University of Wisconsin-Madison
9  Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
10  University of Liverpool & STFC Rutherford Appleton Lab
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 //___________________________________________________________________________
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 //___________________________________________________________________________
272 
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)
void Configure(const Registry &config)
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
TLorentzVector * GetP4(void) const
virtual void SetPolarization(GHepRecord *ev) 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
virtual void AddToEventRecord(GHepRecord *ev, int pdgc, const TLorentzVector &p4) const
bool IsTau(int pdgc)
Definition: PDGUtils.cxx:200
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
virtual void ProcessEventRecord(GHepRecord *evrec) const
#define pINFO
Definition: Messenger.h:63
bool IsMuon(int pdgc)
Definition: PDGUtils.cxx:190
bool PolzIsSet(void) const
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
static const double A
Definition: Units.h:82
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 TVector3 NucRestFrame2Lab(GHepRecord *ev) const
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
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