98{
103
104
105 G4double eMass = CLHEP::electron_mass_c2;
107
108
110
115
116 if (randNumbs[7] <
W[0]) {
117 G4double A1 = -(q2 - 2.*muEnergy*gEnergy);
118 G4double B1 = -(2.*gMomentum*muMomentum);
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];
126
128 dirGamma.
set(sintg*cospg, sintg*sinpg, costg);
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;
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;
140
142 dirMuon.
set(sint1, 0., cost1);
143
144 G4double cost5 = -1. + 2.*randNumbs[6];
145 G4double phi5 = CLHEP::twopi*randNumbs[8];
146
149
152
153 dirElectron = eFourMomentum.
vect().
unit();
154 dirPositron = pFourMomentum.
vect().
unit();
155
156 G4double phi = CLHEP::twopi*randNumbs[3];
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
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);
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
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;
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;
187
189 dirGamma.
set(sint*cosp, sint*sinp, cost);
191
193 dirMuon.
set(sint3, 0., cost3);
194
195 G4double cost5 = -1. + 2.*randNumbs[6];
196 G4double phi5 = CLHEP::twopi*randNumbs[8];
197
200
203
204 dirElectron = eFourMomentum.
vect().
unit();
205 dirPositron = pFourMomentum.
vect().
unit();
206
207 G4double phi = CLHEP::twopi*randNumbs[3];
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];
216 G4double muEnergyInterval = muEnergy - 2.*eMass - minmuFinalEnergy;
217 G4double muFEnergy = minmuFinalEnergy + muEnergyInterval*randNumbs[3];
218
220 G4double maxeEnergy = muEnergy - muFEnergy - eMass;
221 G4double eEnergyInterval = maxeEnergy - mineEnergy;
222 G4double eEnergy = mineEnergy + eEnergyInterval*randNumbs[4];
223
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
246 muFinalFourMomentum.set(muFEnergy, muFinalMomentumVector);
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);
255 G4double t5interval =
G4Log((absA5 + absB5*maxt5)/(absA5 + absB5*mint5))/absB5;
256 G4double argexp = absB5*t5interval*randNumbs[6] +
G4Log(absA5 + absB5*mint5);
260
261 dirElectron.
set(sint5*cosp5, sint5*sinp5, cost5);
263
264 G4ThreeVector auxVec2 = muMomentumVector - muFinalMomentumVector - eMomentumVector;
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;
274 G4double absA6C6 = std::abs(A6 + C6);
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;
280
281 dirPositron.
set(sint6*cosp6, sint6*sinp6, cost6);
282
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];
291 G4double muFinalEnergyinterval = muEnergy - 2.*eMass - minmuFinalEnergy;
292 G4double muFEnergy = minmuFinalEnergy + muFinalEnergyinterval*randNumbs[3];
293
295 G4double maxpEnergy = muEnergy - muFEnergy - eMass;
296 G4double pEnergyinterval = maxpEnergy - minpEnergy;
297 G4double pEnergy = minpEnergy + pEnergyinterval*randNumbs[4];
298
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
319 muFinalMomentumVector.
set(muFMomentum*sint3*cosp3, muFMomentum*sint3*sinp3,
320 muFMomentum*cost3);
321
323 muFinalFourMomentum.set(muFEnergy, muFinalMomentumVector);
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);
332 G4double t6interval =
G4Log((absA6 + absB6*maxt6)/(absA6 + absB6*mint6))/absB6;
333 G4double argexp = absB6*t6interval*randNumbs[6] +
G4Log(absA6 + absB6*mint6);
337
338 dirPositron.
set(sint6*cosp6, sint6*sinp6, cost6);
340
341 G4ThreeVector auxVec2 = muMomentumVector - muFinalMomentumVector - pMomentumVector;
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;
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;
357
358 dirElectron.set(sint5*cosp5, sint5*sinp5, cost5);
359
362 }
363 return muFinalFourMomentum;
364}
G4double C(G4double temp)
G4double B(G4double temperature)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double G4Log(G4double x)
CLHEP::Hep3Vector G4ThreeVector
void set(double x, double y, double z)
HepLorentzVector & boost(double, double, double)
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
G4ThreeVector GetMomentum() const
G4double GetPDGMass() 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)