Geant4 11.3.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4RiGeAngularGenerator Class Reference

#include <G4RiGeAngularGenerator.hh>

+ Inheritance diagram for G4RiGeAngularGenerator:

Public Member Functions

 G4RiGeAngularGenerator ()
 
 ~G4RiGeAngularGenerator () override=default
 
G4ThreeVectorSampleDirection (const G4DynamicParticle *dp, G4double gEnergy, G4int Z, const G4Material *mat=nullptr) override
 
G4LorentzVector Sample5DPairDirections (const G4DynamicParticle *dp, G4ThreeVector &dirElectron, G4ThreeVector &dirPositron, const G4double gEnergy, const G4double q2, const G4double gMomentum, G4double muFinalMomentum, G4double muFinalEnergy, const G4double *randNumbs, const G4double *W)
 
void PhiRotation (G4ThreeVector &dir, G4double phi)
 
G4LorentzVector eDP2 (G4double x1, G4double x2, G4double x3, G4double x4, G4double x5)
 
G4LorentzVector pDP2 (G4double x3, const G4LorentzVector &x6)
 
void PrintGeneratorInformation () const override
 
G4RiGeAngularGeneratoroperator= (const G4RiGeAngularGenerator &right)=delete
 
 G4RiGeAngularGenerator (const G4RiGeAngularGenerator &)=delete
 
- Public Member Functions inherited from G4VEmAngularDistribution
 G4VEmAngularDistribution (const G4String &name)
 
virtual ~G4VEmAngularDistribution ()
 
virtual G4ThreeVectorSampleDirectionForShell (const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
 
virtual void SamplePairDirections (const G4DynamicParticle *dp, G4double elecKinEnergy, G4double posiKinEnergy, G4ThreeVector &dirElectron, G4ThreeVector &dirPositron, G4int Z=0, const G4Material *mat=nullptr)
 
const G4StringGetName () const
 
G4VEmAngularDistributionoperator= (const G4VEmAngularDistribution &right)=delete
 
 G4VEmAngularDistribution (const G4VEmAngularDistribution &)=delete
 

Additional Inherited Members

- Protected Attributes inherited from G4VEmAngularDistribution
G4ThreeVector fLocalDirection
 
G4bool fPolarisation
 

Detailed Description

Definition at line 49 of file G4RiGeAngularGenerator.hh.

Constructor & Destructor Documentation

◆ G4RiGeAngularGenerator() [1/2]

G4RiGeAngularGenerator::G4RiGeAngularGenerator ( )

Definition at line 49 of file G4RiGeAngularGenerator.cc.

50 : G4VEmAngularDistribution("RiGeAngularGen")
51{}
G4VEmAngularDistribution(const G4String &name)

Referenced by G4RiGeAngularGenerator(), and operator=().

◆ ~G4RiGeAngularGenerator()

G4RiGeAngularGenerator::~G4RiGeAngularGenerator ( )
overridedefault

◆ G4RiGeAngularGenerator() [2/2]

G4RiGeAngularGenerator::G4RiGeAngularGenerator ( const G4RiGeAngularGenerator & )
delete

Member Function Documentation

◆ eDP2()

G4LorentzVector G4RiGeAngularGenerator::eDP2 ( G4double x1,
G4double x2,
G4double x3,
G4double x4,
G4double x5 )

Definition at line 383 of file G4RiGeAngularGenerator.cc.

386{
387 G4double sint = std::sqrt((1.0 - x4)*(1.0 + x4));
388 G4double cosp = std::cos(x5);
389 G4double sinp = std::sin(x5);
390
391 G4double QJM2 = (x1 + x3 - x2)*(x1 + x3 - x2)/(4.*x1) - x3;
392
393 if (QJM2 < 0.) {
394 QJM2 = 1.e-13;
395 }
396
397 G4double QJM = std::sqrt(QJM2);
398 G4LorentzVector x6(std::sqrt(x2 + QJM2), QJM*sint*cosp, QJM*sint*sinp, QJM*x4);
399
400 return x6;
401}
CLHEP::HepLorentzVector G4LorentzVector
double G4double
Definition G4Types.hh:83

Referenced by Sample5DPairDirections().

◆ operator=()

G4RiGeAngularGenerator & G4RiGeAngularGenerator::operator= ( const G4RiGeAngularGenerator & right)
delete

◆ pDP2()

G4LorentzVector G4RiGeAngularGenerator::pDP2 ( G4double x3,
const G4LorentzVector & x6 )

Definition at line 405 of file G4RiGeAngularGenerator.cc.

406{
407 G4LorentzVector x7(x3 + x6.vect().dot(x6.vect()), -x6.vect());
408 return x7;
409}
double dot(const Hep3Vector &) const
Hep3Vector vect() const

Referenced by Sample5DPairDirections().

◆ PhiRotation()

void G4RiGeAngularGenerator::PhiRotation ( G4ThreeVector & dir,
G4double phi )

Definition at line 369 of file G4RiGeAngularGenerator.cc.

370{
371 G4double sinp = std::sin(phi);
372 G4double cosp = std::cos(phi);
373
374 G4double newX = dir.x()*cosp + dir.y()*sinp;
375 G4double newY = -dir.x()*sinp + dir.y()*cosp;
376
377 dir.setX(newX);
378 dir.setY(newY);
379}
double x() const
void setY(double)
double y() const
void setX(double)

Referenced by Sample5DPairDirections().

◆ PrintGeneratorInformation()

void G4RiGeAngularGenerator::PrintGeneratorInformation ( ) const
overridevirtual

Reimplemented from G4VEmAngularDistribution.

Definition at line 413 of file G4RiGeAngularGenerator.cc.

414{
415 G4cout << "\n" << G4endl;
416 G4cout << "Angular Generator by RiGe algorithm" << G4endl;
417}
#define G4endl
Definition G4ios.hh:67
G4GLOB_DLL std::ostream G4cout

◆ Sample5DPairDirections()

G4LorentzVector G4RiGeAngularGenerator::Sample5DPairDirections ( const G4DynamicParticle * dp,
G4ThreeVector & dirElectron,
G4ThreeVector & dirPositron,
const G4double gEnergy,
const G4double q2,
const G4double gMomentum,
G4double muFinalMomentum,
G4double muFinalEnergy,
const G4double * randNumbs,
const G4double * W )

Definition at line 89 of file G4RiGeAngularGenerator.cc.

98{
99 G4double muEnergy = dp->GetKineticEnergy();
100 G4ThreeVector muMomentumVector = dp->GetMomentum();
101 G4double muMomentum = muMomentumVector.mag();
102 G4LorentzVector muFinalFourMomentum(muEnergy, muMomentumVector);
103
104 // Electron mass
105 G4double eMass = CLHEP::electron_mass_c2;
106 G4double eMass2 = eMass*eMass;
107
108 // Muon mass
109 G4double muMass = dp->GetDefinition()->GetPDGMass();
110
111 G4double mint3 = 0.;
112 G4double maxt3 = CLHEP::pi;
113 G4double Cmin = std::cos(maxt3);
114 G4double Cmax = std::cos(mint3);
115
116 if (randNumbs[7] < W[0]) {
117 G4double A1 = -(q2 - 2.*muEnergy*gEnergy);
118 G4double B1 = -(2.*gMomentum*muMomentum);
119 G4double tginterval = G4Log((A1 + B1)/(A1 - B1))/B1;
120
121 G4double costg = (-A1 + (A1 - B1)*G4Exp(B1*tginterval*randNumbs[1]))/B1;
122 G4double sintg = std::sqrt((1.0 - costg)*(1.0 + costg));
123 G4double phig = CLHEP::twopi*randNumbs[2];
124 G4double sinpg = std::sin(phig);
125 G4double cospg = std::cos(phig);
126
127 G4ThreeVector dirGamma;
128 dirGamma.set(sintg*cospg, sintg*sinpg, costg);
129 G4LorentzVector gFourMomentum(gEnergy, dirGamma*gMomentum);
130
131 G4double Ap = muMomentum*muMomentum + muFinalMomentum*muFinalMomentum + gMomentum*gMomentum;
132 G4double A = Ap - 2.*muMomentum*gMomentum*costg;
133 G4double B = 2.*muFinalMomentum*gMomentum*sintg*cospg;
134 G4double C = 2.*muFinalMomentum*gMomentum*costg - 2.*muMomentum*muFinalMomentum;
135 G4double absB = std::abs(B);
136 G4double t1interval = (1./(A + C + absB*mint3) - 1./(A + C + absB*maxt3))/absB;
137 G4double t1 = (-(A + C) + 1./(1./(A + C + absB*mint3) - absB*t1interval*randNumbs[0]))/absB;
138 G4double sint1 = std::sin(t1);
139 G4double cost1 = std::cos(t1);
140
141 G4ThreeVector dirMuon;
142 dirMuon.set(sint1, 0., cost1);
143
144 G4double cost5 = -1. + 2.*randNumbs[6];
145 G4double phi5 = CLHEP::twopi*randNumbs[8];
146
147 G4LorentzVector eFourMomentumMQ = eDP2(q2, eMass2, eMass2, cost5, phi5);
148 G4LorentzVector pFourMomentumMQ = pDP2(eMass2, eFourMomentumMQ);
149
150 G4LorentzVector eFourMomentum = eFourMomentumMQ.boost(gFourMomentum.boostVector());
151 G4LorentzVector pFourMomentum = pFourMomentumMQ.boost(gFourMomentum.boostVector());
152
153 dirElectron = eFourMomentum.vect().unit();
154 dirPositron = pFourMomentum.vect().unit();
155
156 G4double phi = CLHEP::twopi*randNumbs[3];
157 PhiRotation(dirElectron, phi);
158 PhiRotation(dirPositron, phi);
159
160 } else if (randNumbs[7] >= W[0] && randNumbs[7] < W[1]) {
161 G4double A3 = q2 + 2.*gEnergy*muFinalEnergy;
162 G4double B3 = -2.*gMomentum*muFinalMomentum;
163
164 G4double tQ3interval = G4Log((A3 + B3)/(A3 - B3))/B3;
165 G4double tQMG = (-A3 + (A3 - B3)*G4Exp(B3*tQ3interval*randNumbs[0]))/B3;
166 G4double phiQP = CLHEP::twopi*randNumbs[2];
167
168 G4double sintQ3 = std::sqrt(1. - tQMG*tQMG);
169 G4double cospQP = std::cos(phiQP);
170 G4double sinpQP = std::sin(phiQP);
171
172 G4double Ap = muMomentum*muMomentum + muFinalMomentum*muFinalMomentum + gMomentum*gMomentum;
173 G4double A = Ap + 2.*muFinalMomentum*gMomentum*tQMG;
174 G4double B = -2.*muMomentum*gMomentum*sintQ3*cospQP;
175 G4double C = -2.*muMomentum*gMomentum*tQMG - 2.*muMomentum*muFinalMomentum;
176
177 G4double absB = std::abs(B);
178 G4double t3interval = (1./(A + C + absB*mint3) - 1./(A + C + absB*maxt3))/absB;
179 G4double t3 = (-(A + C) + 1./(1./(A + C + absB*mint3) - absB*t3interval*randNumbs[0]))/absB;
180 G4double sint3 = std::sin(t3);
181 G4double cost3 = std::cos(t3);
182
183 G4double cost = -sint3*sintQ3*cospQP + cost3*tQMG;
184 G4double sint = std::sqrt((1. + cost)*(1. - cost));
185 G4double cosp = (sintQ3*cospQP*cost3 + sint3*tQMG)/sint;
186 G4double sinp = sintQ3*sinpQP/sint;
187
188 G4ThreeVector dirGamma;
189 dirGamma.set(sint*cosp, sint*sinp, cost);
190 G4LorentzVector gFourMomentum(gEnergy, dirGamma*gMomentum);
191
192 G4ThreeVector dirMuon;
193 dirMuon.set(sint3, 0., cost3);
194
195 G4double cost5 = -1. + 2.*randNumbs[6];
196 G4double phi5 = CLHEP::twopi*randNumbs[8];
197
198 G4LorentzVector eFourMomentumMQ = eDP2(q2, eMass2, eMass2, cost5, phi5);
199 G4LorentzVector pFourMomentumMQ = pDP2(eMass2, eFourMomentumMQ);
200
201 G4LorentzVector eFourMomentum = eFourMomentumMQ.boost(gFourMomentum.boostVector());
202 G4LorentzVector pFourMomentum = pFourMomentumMQ.boost(gFourMomentum.boostVector());
203
204 dirElectron = eFourMomentum.vect().unit();
205 dirPositron = pFourMomentum.vect().unit();
206
207 G4double phi = CLHEP::twopi*randNumbs[3];
208 PhiRotation(dirElectron, phi);
209 PhiRotation(dirPositron, phi);
210
211 } else if (randNumbs[7] >= W[1] && randNumbs[7] < W[2]) {
212 G4double phi3 = CLHEP::twopi*randNumbs[0];
213 G4double phi5 = CLHEP::twopi*randNumbs[1];
214 G4double phi6 = CLHEP::twopi*randNumbs[2];
215 G4double minmuFinalEnergy = muMass;
216 G4double muEnergyInterval = muEnergy - 2.*eMass - minmuFinalEnergy;
217 G4double muFEnergy = minmuFinalEnergy + muEnergyInterval*randNumbs[3];
218
219 G4double mineEnergy = eMass;
220 G4double maxeEnergy = muEnergy - muFEnergy - eMass;
221 G4double eEnergyInterval = maxeEnergy - mineEnergy;
222 G4double eEnergy = mineEnergy + eEnergyInterval*randNumbs[4];
223
224 G4double cosp3 = 1.;
225 G4double sinp3 = 0.;
226 G4double cosp5 = std::cos(phi5);
227 G4double sinp5 = std::sin(phi5);
228 G4double cosp6 = std::cos(phi6);
229 G4double sinp6 = std::sin(phi6);
230
231 G4double muFMomentum = std::sqrt(muFEnergy*muFEnergy - muMass*muMass);
232 G4double eMomentum = std::sqrt(eEnergy*eEnergy - eMass*eMass);
233 G4double pEnergy = muEnergy - muFEnergy - eEnergy;
234 G4double pMomentum = std::sqrt(pEnergy*pEnergy - eMass*eMass);
235
236 G4double A3 = -2.*muMass*muMass + 2.*muEnergy*muFinalEnergy;
237 G4double B3 = -2.*muMomentum*muFinalMomentum;
238 G4double cost3interval = G4Log((A3 + B3*Cmax)/(A3 + B3*Cmin))/B3;
239 G4double expanCost3r6 = G4Exp(B3*cost3interval*randNumbs[5]);
240 G4double cost3 = A3*(expanCost3r6 - 1.)/B3 + Cmin*expanCost3r6;
241 G4double sint3 = std::sqrt((1. - cost3)*(1. + cost3));
242
243 G4ThreeVector muFinalMomentumVector(muFMomentum*sint3, 0., muFMomentum*cost3);
244
245 G4LorentzVector muFourMomentum(muMomentum, muMomentumVector);
246 muFinalFourMomentum.set(muFEnergy, muFinalMomentumVector);
247 G4LorentzVector auxVec1 = muFourMomentum - muFinalFourMomentum;
248 G4double A5 = auxVec1.mag2() - 2.*eEnergy*(muEnergy - muFEnergy) +
249 2.*muMomentumVector[2]*eMomentum - 2.*muFMomentum*eMomentum*cost3;
250 G4double B5 = -2.*muFMomentum*eMomentum*(sint3*cosp3*cosp5 + sint3*sinp3*sinp5);
251 G4double absA5 = std::abs(A5);
252 G4double absB5 = std::abs(B5);
253 G4double mint5 = 0.;
254 G4double maxt5 = CLHEP::pi;
255 G4double t5interval = G4Log((absA5 + absB5*maxt5)/(absA5 + absB5*mint5))/absB5;
256 G4double argexp = absB5*t5interval*randNumbs[6] + G4Log(absA5 + absB5*mint5);
257 G4double t5 = -absA5/absB5 + G4Exp(argexp)/absB5;
258 G4double sint5 = std::sin(t5);
259 G4double cost5 = std::cos(t5);
260
261 dirElectron.set(sint5*cosp5, sint5*sinp5, cost5);
262 G4ThreeVector eMomentumVector = eMomentum*dirElectron;
263
264 G4ThreeVector auxVec2 = muMomentumVector - muFinalMomentumVector - eMomentumVector;
265 G4double p1mp3mp52 = auxVec2.dot(auxVec2);
266 G4double Bp = muFinalMomentum*(sint3*cosp3*cosp6 + sint3*sinp3*sinp6) +
267 eMomentum*(sint5*cosp5*cosp6 + sint5*sinp5*sinp6);
268 G4double Cp = -muMomentum + muFMomentum*cost3 + eMomentum*cost5;
269 G4double A6 = p1mp3mp52 + pMomentum*pMomentum;
270 G4double B6 = 2.*pMomentum*Bp;
271 G4double C6 = 2.*pMomentum*Cp;
272 G4double mint6 = 0.;
273 G4double maxt6 = CLHEP::pi;
274 G4double absA6C6 = std::abs(A6 + C6);
275 G4double absB6 = std::abs(B6);
276 G4double t6interval = (1./(absA6C6 + absB6*mint6) - 1./(absA6C6 + absB6*maxt6))/absB6;
277 G4double t6 = (-absA6C6 + 1./(1./(absA6C6 + absB6*mint6) - absB6*t6interval*randNumbs[8]))/absB6;
278 G4double sint6 = std::sin(t6);
279 G4double cost6 = std::cos(t6);
280
281 dirPositron.set(sint6*cosp6, sint6*sinp6, cost6);
282
283 PhiRotation(dirElectron, phi3);
284 PhiRotation(dirPositron, phi3);
285
286 } else if (randNumbs[7] >= W[2]) {
287 G4double phi3 = CLHEP::twopi*randNumbs[0];
288 G4double phi6 = CLHEP::twopi*randNumbs[1];
289 G4double phi5 = CLHEP::twopi*randNumbs[2];
290 G4double minmuFinalEnergy = muMass;
291 G4double muFinalEnergyinterval = muEnergy - 2.*eMass - minmuFinalEnergy;
292 G4double muFEnergy = minmuFinalEnergy + muFinalEnergyinterval*randNumbs[3];
293
294 G4double minpEnergy = eMass;
295 G4double maxpEnergy = muEnergy - muFEnergy - eMass;
296 G4double pEnergyinterval = maxpEnergy - minpEnergy;
297 G4double pEnergy = minpEnergy + pEnergyinterval*randNumbs[4];
298
299 G4double cosp3 = 1.;
300 G4double sinp3 = 0.;
301 G4double cosp5 = std::cos(phi5);
302 G4double sinp5 = std::sin(phi5);
303 G4double cosp6 = std::cos(phi6);
304 G4double sinp6 = std::sin(phi6);
305
306 G4double muFMomentum = std::sqrt(muFEnergy*muFEnergy - muMass*muMass);
307 G4double pMomentum = std::sqrt(pEnergy*pEnergy - eMass*eMass);
308 G4double eEnergy = muEnergy - muFEnergy - pEnergy;
309 G4double eMomentum = std::sqrt(eEnergy*eEnergy - eMass*eMass);
310
311 G4double A3 = -2.*muMass*muMass + 2.*muEnergy*muFinalEnergy;
312 G4double B3 = -2.*muMomentum*muFMomentum;
313 G4double cost3interval = G4Log((A3 + B3*Cmax)/(A3 + B3*Cmin))/B3;
314 G4double expanCost3r6 = G4Exp(B3*cost3interval*randNumbs[5]);
315 G4double cost3 = A3*(expanCost3r6 - 1.)/B3 + Cmin*expanCost3r6;
316 G4double sint3 = std::sqrt((1. - cost3)*(1. + cost3));
317
318 G4ThreeVector muFinalMomentumVector;
319 muFinalMomentumVector.set(muFMomentum*sint3*cosp3, muFMomentum*sint3*sinp3,
320 muFMomentum*cost3);
321
322 G4LorentzVector muFourMomentum(muMomentum, muMomentumVector);
323 muFinalFourMomentum.set(muFEnergy, muFinalMomentumVector);
324 G4LorentzVector auxVec1 = muFourMomentum - muFinalFourMomentum;
325 G4double A6 = auxVec1.mag2() - 2.*pEnergy*(muEnergy - muFEnergy) +
326 2.*muMomentumVector[2]*pMomentum - 2.*muFMomentum*pMomentum*cost3;
327 G4double B6 = -2.*muFMomentum*pMomentum*(sint3*cosp3*cosp6 + sint3*sinp3*sinp6);
328 G4double absA6 = std::abs(A6);
329 G4double absB6 = std::abs(B6);
330 G4double mint6 = 0.;
331 G4double maxt6 = CLHEP::pi;
332 G4double t6interval = G4Log((absA6 + absB6*maxt6)/(absA6 + absB6*mint6))/absB6;
333 G4double argexp = absB6*t6interval*randNumbs[6] + G4Log(absA6 + absB6*mint6);
334 G4double t6 = -absA6/absB6 + G4Exp(argexp)/absB6;
335 G4double sint6 = std::sin(t6);
336 G4double cost6 = std::cos(t6);
337
338 dirPositron.set(sint6*cosp6, sint6*sinp6, cost6);
339 G4ThreeVector pMomentumVector = pMomentum*dirPositron;
340
341 G4ThreeVector auxVec2 = muMomentumVector - muFinalMomentumVector - pMomentumVector;
342 G4double p1mp3mp62 = auxVec2.dot(auxVec2);
343 G4double Bp = muFMomentum*(sint3*cosp3*cosp5 + sint3*sinp3*sinp5) +
344 pMomentum*(sint6*cosp6*cosp5 + sint6*sinp6*sinp5);
345 G4double Cp = -muMomentum + muFMomentum*cost3 + pMomentum*cost6;
346 G4double A5 = p1mp3mp62 + eMomentum*eMomentum;
347 G4double B5 = 2.*eMomentum*Bp;
348 G4double C5 = 2.*eMomentum*Cp;
349 G4double mint5 = 0.;
350 G4double maxt5 = CLHEP::pi;
351 G4double absA5C5 = std::abs(A5 + C5);
352 G4double absB5 = std::abs(B5);
353 G4double t5interval = (1./(absA5C5 + absB5*mint5) - 1./(absA5C5 + absB5*maxt5))/absB5;
354 G4double t5 = (-absA5C5 + 1./(1./(absA5C5 + absB5*mint5) - absB5*t5interval*randNumbs[8]))/absB5;
355 G4double sint5 = std::sin(t5);
356 G4double cost5 = std::cos(t5);
357
358 dirElectron.set(sint5*cosp5, sint5*sinp5, cost5);
359
360 PhiRotation(dirElectron, phi3);
361 PhiRotation(dirPositron, phi3);
362 }
363 return muFinalFourMomentum;
364}
G4double C(G4double temp)
G4double B(G4double temperature)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:180
G4double G4Log(G4double x)
Definition G4Log.hh:227
CLHEP::Hep3Vector G4ThreeVector
const G4double A[17]
#define C5
Hep3Vector unit() const
double mag() const
void set(double x, double y, double z)
HepLorentzVector & boost(double, double, double)
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
G4ThreeVector GetMomentum() const
void PhiRotation(G4ThreeVector &dir, G4double phi)
G4LorentzVector eDP2(G4double x1, G4double x2, G4double x3, G4double x4, G4double x5)
G4LorentzVector pDP2(G4double x3, const G4LorentzVector &x6)
#define W
Definition crc32.c:85

◆ SampleDirection()

G4ThreeVector & G4RiGeAngularGenerator::SampleDirection ( const G4DynamicParticle * dp,
G4double gEnergy,
G4int Z,
const G4Material * mat = nullptr )
overridevirtual

Implements G4VEmAngularDistribution.

Definition at line 56 of file G4RiGeAngularGenerator.cc.

59{
60 // Sample gamma angle (Z - axis along the parent particle).
61 G4double cost = SampleCosTheta(dp->GetKineticEnergy(), gEnergy,
62 dp->GetDefinition()->GetPDGMass());
63 G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
64 G4double phi = CLHEP::twopi*G4UniformRand();
65
66 fLocalDirection.set(sint*std::cos(phi), sint*std::sin(phi), cost);
68
69 return fLocalDirection;
70}
#define G4UniformRand()
Definition Randomize.hh:52
const G4ThreeVector & GetMomentumDirection() const

The documentation for this class was generated from the following files: