1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // |
---|
27 | // $Id: G4StatMFChannel.cc,v 1.10 2008/11/19 14:33:31 vnivanch Exp $ |
---|
28 | // GEANT4 tag $Name: geant4-09-02 $ |
---|
29 | // |
---|
30 | // Hadronic Process: Nuclear De-excitations |
---|
31 | // by V. Lara |
---|
32 | // |
---|
33 | // Modified: |
---|
34 | // 25.07.08 I.Pshenichnov (in collaboration with Alexander Botvina and Igor |
---|
35 | // Mishustin (FIAS, Frankfurt, INR, Moscow and Kurchatov Institute, |
---|
36 | // Moscow, pshenich@fias.uni-frankfurt.de) fixed semi-infinite loop |
---|
37 | |
---|
38 | |
---|
39 | #include "G4StatMFChannel.hh" |
---|
40 | #include "G4HadronicException.hh" |
---|
41 | #include <numeric> |
---|
42 | |
---|
43 | class SumCoulombEnergy : public std::binary_function<G4double,G4double,G4double> |
---|
44 | { |
---|
45 | public: |
---|
46 | SumCoulombEnergy() : total(0.0) {} |
---|
47 | G4double operator() (G4double& , G4StatMFFragment*& frag) |
---|
48 | { |
---|
49 | total += frag->GetCoulombEnergy(); |
---|
50 | return total; |
---|
51 | } |
---|
52 | |
---|
53 | G4double GetTotal() { return total; } |
---|
54 | public: |
---|
55 | G4double total; |
---|
56 | }; |
---|
57 | |
---|
58 | |
---|
59 | |
---|
60 | |
---|
61 | |
---|
62 | // Copy constructor |
---|
63 | G4StatMFChannel::G4StatMFChannel(const G4StatMFChannel & ) |
---|
64 | { |
---|
65 | throw G4HadronicException(__FILE__, __LINE__, "G4StatMFChannel::copy_constructor meant to not be accessable"); |
---|
66 | } |
---|
67 | |
---|
68 | // Operators |
---|
69 | |
---|
70 | G4StatMFChannel & G4StatMFChannel:: |
---|
71 | operator=(const G4StatMFChannel & ) |
---|
72 | { |
---|
73 | throw G4HadronicException(__FILE__, __LINE__, "G4StatMFChannel::operator= meant to not be accessable"); |
---|
74 | return *this; |
---|
75 | } |
---|
76 | |
---|
77 | |
---|
78 | G4bool G4StatMFChannel::operator==(const G4StatMFChannel & ) const |
---|
79 | { |
---|
80 | // throw G4HadronicException(__FILE__, __LINE__, "G4StatMFChannel::operator== meant to not be accessable"); |
---|
81 | return false; |
---|
82 | } |
---|
83 | |
---|
84 | |
---|
85 | G4bool G4StatMFChannel::operator!=(const G4StatMFChannel & ) const |
---|
86 | { |
---|
87 | // throw G4HadronicException(__FILE__, __LINE__, "G4StatMFChannel::operator!= meant to not be accessable"); |
---|
88 | return true; |
---|
89 | } |
---|
90 | |
---|
91 | |
---|
92 | G4bool G4StatMFChannel::CheckFragments(void) |
---|
93 | { |
---|
94 | std::deque<G4StatMFFragment*>::iterator i; |
---|
95 | for (i = _theFragments.begin(); |
---|
96 | i != _theFragments.end(); ++i) |
---|
97 | { |
---|
98 | G4int A = static_cast<G4int>((*i)->GetA()); |
---|
99 | G4int Z = static_cast<G4int>((*i)->GetZ()); |
---|
100 | if ( (A > 1 && (Z > A || Z <= 0)) || (A==1 && Z > A) || A <= 0 ) return false; |
---|
101 | } |
---|
102 | |
---|
103 | return true; |
---|
104 | } |
---|
105 | |
---|
106 | |
---|
107 | |
---|
108 | |
---|
109 | void G4StatMFChannel::CreateFragment(const G4double A, const G4double Z) |
---|
110 | // Create a new fragment. |
---|
111 | // Fragments are automatically sorted: first charged fragments, |
---|
112 | // then neutral ones. |
---|
113 | { |
---|
114 | if (Z <= 0.5) { |
---|
115 | _theFragments.push_back(new G4StatMFFragment(static_cast<G4int>(A),static_cast<G4int>(Z))); |
---|
116 | _NumOfNeutralFragments++; |
---|
117 | } else { |
---|
118 | _theFragments.push_front(new G4StatMFFragment(static_cast<G4int>(A),static_cast<G4int>(Z))); |
---|
119 | _NumOfChargedFragments++; |
---|
120 | } |
---|
121 | |
---|
122 | return; |
---|
123 | } |
---|
124 | |
---|
125 | |
---|
126 | G4double G4StatMFChannel::GetFragmentsCoulombEnergy(void) |
---|
127 | { |
---|
128 | G4double Coulomb = std::accumulate(_theFragments.begin(),_theFragments.end(), |
---|
129 | 0.0,SumCoulombEnergy()); |
---|
130 | // G4double Coulomb = 0.0; |
---|
131 | // for (unsigned int i = 0;i < _theFragments.size(); i++) |
---|
132 | // Coulomb += _theFragments[i]->GetCoulombEnergy(); |
---|
133 | return Coulomb; |
---|
134 | } |
---|
135 | |
---|
136 | |
---|
137 | |
---|
138 | G4double G4StatMFChannel::GetFragmentsEnergy(const G4double T) const |
---|
139 | { |
---|
140 | G4double Energy = 0.0; |
---|
141 | |
---|
142 | G4double TranslationalEnergy = (3./2.)*T*static_cast<G4double>(_theFragments.size()); |
---|
143 | |
---|
144 | std::deque<G4StatMFFragment*>::const_iterator i; |
---|
145 | for (i = _theFragments.begin(); i != _theFragments.end(); ++i) |
---|
146 | { |
---|
147 | Energy += (*i)->GetEnergy(T); |
---|
148 | } |
---|
149 | return Energy + TranslationalEnergy; |
---|
150 | } |
---|
151 | |
---|
152 | G4FragmentVector * G4StatMFChannel::GetFragments(const G4double anA, |
---|
153 | const G4double anZ, |
---|
154 | const G4double T) |
---|
155 | // |
---|
156 | { |
---|
157 | // calculate momenta of charged fragments |
---|
158 | CoulombImpulse(anA,anZ,T); |
---|
159 | |
---|
160 | // calculate momenta of neutral fragments |
---|
161 | FragmentsMomenta(_NumOfNeutralFragments, _NumOfChargedFragments, T); |
---|
162 | |
---|
163 | |
---|
164 | G4FragmentVector * theResult = new G4FragmentVector; |
---|
165 | std::deque<G4StatMFFragment*>::iterator i; |
---|
166 | for (i = _theFragments.begin(); i != _theFragments.end(); ++i) |
---|
167 | theResult->push_back((*i)->GetFragment(T)); |
---|
168 | |
---|
169 | return theResult; |
---|
170 | |
---|
171 | } |
---|
172 | |
---|
173 | |
---|
174 | |
---|
175 | void G4StatMFChannel::CoulombImpulse(const G4double anA, const G4double anZ, const G4double T) |
---|
176 | // Aafter breakup, fragments fly away under Coulomb field. |
---|
177 | // This method calculates asymptotic fragments momenta. |
---|
178 | { |
---|
179 | // First, we have to place the fragments inside of the original nucleus volume |
---|
180 | PlaceFragments(anA); |
---|
181 | |
---|
182 | // Second, we sample initial charged fragments momenta. There are |
---|
183 | // _NumOfChargedFragments charged fragments and they start at the begining |
---|
184 | // of the vector _theFragments (i.e. 0) |
---|
185 | FragmentsMomenta(_NumOfChargedFragments, 0, T); |
---|
186 | |
---|
187 | // Third, we have to figure out the asymptotic momenta of charged fragments |
---|
188 | // For taht we have to solve equations of motion for fragments |
---|
189 | SolveEqOfMotion(anA,anZ,T); |
---|
190 | |
---|
191 | return; |
---|
192 | } |
---|
193 | |
---|
194 | |
---|
195 | |
---|
196 | void G4StatMFChannel::PlaceFragments(const G4double anA) |
---|
197 | // This gives the position of fragments at the breakup instant. |
---|
198 | // Fragments positions are sampled inside prolongated ellipsoid. |
---|
199 | { |
---|
200 | const G4double R0 = G4StatMFParameters::Getr0(); |
---|
201 | const G4double Rsys = 2.0*R0*std::pow(anA,1./3.); |
---|
202 | |
---|
203 | G4bool TooMuchIterations; |
---|
204 | do |
---|
205 | { |
---|
206 | TooMuchIterations = false; |
---|
207 | |
---|
208 | // Sample the position of the first fragment |
---|
209 | G4double R = (Rsys - R0*std::pow(_theFragments[0]->GetA(),1./3.))* |
---|
210 | std::pow(G4UniformRand(),1./3.); |
---|
211 | _theFragments[0]->SetPosition(IsotropicVector(R)); |
---|
212 | |
---|
213 | |
---|
214 | // Sample the position of the remaining fragments |
---|
215 | G4bool ThereAreOverlaps = false; |
---|
216 | std::deque<G4StatMFFragment*>::iterator i; |
---|
217 | for (i = _theFragments.begin()+1; i != _theFragments.end(); ++i) |
---|
218 | { |
---|
219 | G4int counter = 0; |
---|
220 | do |
---|
221 | { |
---|
222 | R = (Rsys - R0*std::pow((*i)->GetA(),1./3.))*std::pow(G4UniformRand(),1./3.); |
---|
223 | (*i)->SetPosition(IsotropicVector(R)); |
---|
224 | |
---|
225 | // Check that there are not overlapping fragments |
---|
226 | std::deque<G4StatMFFragment*>::iterator j; |
---|
227 | for (j = _theFragments.begin(); j != i; ++j) |
---|
228 | { |
---|
229 | G4ThreeVector FragToFragVector = (*i)->GetPosition() - (*j)->GetPosition(); |
---|
230 | G4double Rmin = R0*(std::pow((*i)->GetA(),1./3.) + |
---|
231 | std::pow((*j)->GetA(),1./3)); |
---|
232 | if ( (ThereAreOverlaps = (FragToFragVector.mag2() < Rmin*Rmin)) ) break; |
---|
233 | } |
---|
234 | counter++; |
---|
235 | } while (ThereAreOverlaps && counter < 1000); |
---|
236 | |
---|
237 | if (counter >= 1000) |
---|
238 | { |
---|
239 | TooMuchIterations = true; |
---|
240 | break; |
---|
241 | } |
---|
242 | } |
---|
243 | } while (TooMuchIterations); |
---|
244 | |
---|
245 | return; |
---|
246 | } |
---|
247 | |
---|
248 | |
---|
249 | void G4StatMFChannel::FragmentsMomenta(const G4int NF, const G4int idx, |
---|
250 | const G4double T) |
---|
251 | // Calculate fragments momenta at the breakup instant |
---|
252 | // Fragment kinetic energies are calculated according to the |
---|
253 | // Boltzmann distribution at given temperature. |
---|
254 | // NF is number of fragments |
---|
255 | // idx is index of first fragment |
---|
256 | { |
---|
257 | G4double KinE = (3./2.)*T*static_cast<G4double>(NF); |
---|
258 | |
---|
259 | G4ThreeVector p; |
---|
260 | |
---|
261 | if (NF <= 0) return; |
---|
262 | else if (NF == 1) |
---|
263 | { |
---|
264 | // We have only one fragment to deal with |
---|
265 | p = IsotropicVector(std::sqrt(2.0*_theFragments[idx]->GetNuclearMass()*KinE)); |
---|
266 | _theFragments[idx]->SetMomentum(p); |
---|
267 | } |
---|
268 | else if (NF == 2) |
---|
269 | { |
---|
270 | // We have only two fragment to deal with |
---|
271 | G4double M1 = _theFragments[idx]->GetNuclearMass(); |
---|
272 | G4double M2 = _theFragments[idx+1]->GetNuclearMass(); |
---|
273 | p = IsotropicVector(std::sqrt(2.0*KinE*(M1*M2)/(M1+M2))); |
---|
274 | _theFragments[idx]->SetMomentum(p); |
---|
275 | _theFragments[idx+1]->SetMomentum(-p); |
---|
276 | } |
---|
277 | else |
---|
278 | { |
---|
279 | // We have more than two fragments |
---|
280 | G4double AvailableE; |
---|
281 | G4int i1,i2; |
---|
282 | G4double SummedE; |
---|
283 | G4ThreeVector SummedP; |
---|
284 | do |
---|
285 | { |
---|
286 | // Fisrt sample momenta of NF-2 fragments |
---|
287 | // according to Boltzmann distribution |
---|
288 | AvailableE = 0.0; |
---|
289 | SummedE = 0.0; |
---|
290 | SummedP.setX(0.0);SummedP.setY(0.0);SummedP.setZ(0.0); |
---|
291 | for (G4int i = idx; i < idx+NF-2; i++) |
---|
292 | { |
---|
293 | G4double E; |
---|
294 | G4double RandE; |
---|
295 | G4double Boltzmann; |
---|
296 | do |
---|
297 | { |
---|
298 | E = 9.0*T*G4UniformRand(); |
---|
299 | Boltzmann = std::sqrt(E)*std::exp(-E/T); |
---|
300 | RandE = std::sqrt(T/2.)*std::exp(-0.5)*G4UniformRand(); |
---|
301 | } |
---|
302 | while (RandE > Boltzmann); |
---|
303 | p = IsotropicVector(std::sqrt(2.0*E*_theFragments[i]->GetNuclearMass())); |
---|
304 | _theFragments[i]->SetMomentum(p); |
---|
305 | SummedE += E; |
---|
306 | SummedP += p; |
---|
307 | } |
---|
308 | // Calculate momenta of last two fragments in such a way |
---|
309 | // that constraints are satisfied |
---|
310 | i1 = idx+NF-2; // before last fragment index |
---|
311 | i2 = idx+NF-1; // last fragment index |
---|
312 | p = -SummedP; |
---|
313 | AvailableE = KinE - SummedE; |
---|
314 | // Available Kinetic Energy should be shared between two last fragments |
---|
315 | } |
---|
316 | while (AvailableE <= p.mag2()/(2.0*(_theFragments[i1]->GetNuclearMass()+ |
---|
317 | _theFragments[i2]->GetNuclearMass()))); |
---|
318 | |
---|
319 | G4double H = 1.0 + _theFragments[i2]->GetNuclearMass()/_theFragments[i1]->GetNuclearMass(); |
---|
320 | G4double CTM12 = H*(1.0 - 2.0*_theFragments[i2]->GetNuclearMass()*AvailableE/p.mag2()); |
---|
321 | G4double CosTheta1; |
---|
322 | G4double Sign; |
---|
323 | |
---|
324 | if (CTM12 > 0.9999) {CosTheta1 = 1.;} |
---|
325 | else { |
---|
326 | do |
---|
327 | { |
---|
328 | do |
---|
329 | { |
---|
330 | CosTheta1 = 1.0 - 2.0*G4UniformRand(); |
---|
331 | } |
---|
332 | while (CosTheta1*CosTheta1 < CTM12); |
---|
333 | } |
---|
334 | while (CTM12 >= 0.0 && CosTheta1 < 0.0); |
---|
335 | } |
---|
336 | |
---|
337 | if (CTM12 < 0.0) Sign = 1.0; |
---|
338 | else if (G4UniformRand() <= 0.5) Sign = -1.0; |
---|
339 | else Sign = 1.0; |
---|
340 | |
---|
341 | |
---|
342 | G4double P1 = (p.mag()*CosTheta1+Sign*std::sqrt(p.mag2()*(CosTheta1*CosTheta1-CTM12)))/H; |
---|
343 | G4double P2 = std::sqrt(P1*P1+p.mag2() - 2.0*P1*p.mag()*CosTheta1); |
---|
344 | G4double Phi = twopi*G4UniformRand(); |
---|
345 | G4double SinTheta1 = std::sqrt(1.0 - CosTheta1*CosTheta1); |
---|
346 | G4double CosPhi1 = std::cos(Phi); |
---|
347 | G4double SinPhi1 = std::sin(Phi); |
---|
348 | G4double CosPhi2 = -CosPhi1; |
---|
349 | G4double SinPhi2 = -SinPhi1; |
---|
350 | G4double CosTheta2 = (p.mag2() + P2*P2 - P1*P1)/(2.0*p.mag()*P2); |
---|
351 | G4double SinTheta2 = 0.0; |
---|
352 | if (CosTheta2 > -1.0 && CosTheta2 < 1.0) SinTheta2 = std::sqrt(1.0 - CosTheta2*CosTheta2); |
---|
353 | |
---|
354 | G4ThreeVector p1(P1*SinTheta1*CosPhi1,P1*SinTheta1*SinPhi1,P1*CosTheta1); |
---|
355 | G4ThreeVector p2(P2*SinTheta2*CosPhi2,P2*SinTheta2*SinPhi2,P2*CosTheta2); |
---|
356 | G4ThreeVector b(1.0,0.0,0.0); |
---|
357 | |
---|
358 | p1 = RotateMomentum(p,b,p1); |
---|
359 | p2 = RotateMomentum(p,b,p2); |
---|
360 | |
---|
361 | SummedP += p1 + p2; |
---|
362 | SummedE += p1.mag2()/(2.0*_theFragments[i1]->GetNuclearMass()) + |
---|
363 | p2.mag2()/(2.0*_theFragments[i2]->GetNuclearMass()); |
---|
364 | |
---|
365 | _theFragments[i1]->SetMomentum(p1); |
---|
366 | _theFragments[i2]->SetMomentum(p2); |
---|
367 | |
---|
368 | } |
---|
369 | |
---|
370 | return; |
---|
371 | } |
---|
372 | |
---|
373 | |
---|
374 | void G4StatMFChannel::SolveEqOfMotion(const G4double anA, const G4double anZ, const G4double T) |
---|
375 | // This method will find a solution of Newton's equation of motion |
---|
376 | // for fragments in the self-consistent time-dependent Coulomb field |
---|
377 | { |
---|
378 | G4double CoulombEnergy = (3./5.)*(elm_coupling*anZ*anZ)* |
---|
379 | std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1./3.)/ |
---|
380 | (G4StatMFParameters::Getr0()*std::pow(anA,1./3.)) |
---|
381 | - GetFragmentsCoulombEnergy(); |
---|
382 | if (CoulombEnergy <= 0.0) return; |
---|
383 | |
---|
384 | G4int Iterations = 0; |
---|
385 | G4double TimeN = 0.0; |
---|
386 | G4double TimeS = 0.0; |
---|
387 | G4double DeltaTime = 10.0; |
---|
388 | |
---|
389 | G4ThreeVector * Pos = new G4ThreeVector[_NumOfChargedFragments]; |
---|
390 | G4ThreeVector * Vel = new G4ThreeVector[_NumOfChargedFragments]; |
---|
391 | G4ThreeVector * Accel = new G4ThreeVector[_NumOfChargedFragments]; |
---|
392 | |
---|
393 | G4int i; |
---|
394 | for (i = 0; i < _NumOfChargedFragments; i++) |
---|
395 | { |
---|
396 | Vel[i] = (1.0/(_theFragments[i]->GetNuclearMass()))* |
---|
397 | _theFragments[i]->GetMomentum(); |
---|
398 | Pos[i] = _theFragments[i]->GetPosition(); |
---|
399 | } |
---|
400 | |
---|
401 | do |
---|
402 | { |
---|
403 | |
---|
404 | G4ThreeVector distance; |
---|
405 | G4ThreeVector force; |
---|
406 | |
---|
407 | for (i = 0; i < _NumOfChargedFragments; i++) |
---|
408 | { |
---|
409 | force.setX(0.0); force.setY(0.0); force.setZ(0.0); |
---|
410 | for (G4int j = 0; j < _NumOfChargedFragments; j++) |
---|
411 | { |
---|
412 | if (i != j) |
---|
413 | { |
---|
414 | distance = Pos[i] - Pos[j]; |
---|
415 | force += (elm_coupling*(_theFragments[i]->GetZ()*_theFragments[j]->GetZ())/ |
---|
416 | (distance.mag2()*distance.mag()))*distance; |
---|
417 | } |
---|
418 | } |
---|
419 | Accel[i] = (1./(_theFragments[i]->GetNuclearMass()))*force; |
---|
420 | } |
---|
421 | |
---|
422 | TimeN = TimeS + DeltaTime; |
---|
423 | |
---|
424 | G4ThreeVector SavedVel; |
---|
425 | for ( i = 0; i < _NumOfChargedFragments; i++) |
---|
426 | { |
---|
427 | SavedVel = Vel[i]; |
---|
428 | Vel[i] += Accel[i]*(TimeN-TimeS); |
---|
429 | Pos[i] += (SavedVel+Vel[i])*(TimeN-TimeS)*0.5; |
---|
430 | } |
---|
431 | |
---|
432 | // if (Iterations >= 50 && Iterations < 75) DeltaTime = 4.; |
---|
433 | // else if (Iterations >= 75) DeltaTime = 10.; |
---|
434 | |
---|
435 | TimeS = TimeN; |
---|
436 | |
---|
437 | } |
---|
438 | while (Iterations++ < 100); |
---|
439 | |
---|
440 | // Summed fragment kinetic energy |
---|
441 | G4double TotalKineticEnergy = 0.0; |
---|
442 | for (i = 0; i < _NumOfChargedFragments; i++) |
---|
443 | { |
---|
444 | TotalKineticEnergy += _theFragments[i]->GetNuclearMass()* |
---|
445 | 0.5*Vel[i].mag2(); |
---|
446 | } |
---|
447 | // Scaling of fragment velocities |
---|
448 | G4double KineticEnergy = (3./2.)*static_cast<G4double>(_theFragments.size())*T; |
---|
449 | G4double Eta = ( CoulombEnergy + KineticEnergy ) / TotalKineticEnergy; |
---|
450 | for (i = 0; i < _NumOfChargedFragments; i++) |
---|
451 | { |
---|
452 | Vel[i] *= Eta; |
---|
453 | } |
---|
454 | |
---|
455 | // Finally calculate fragments momenta |
---|
456 | for (i = 0; i < _NumOfChargedFragments; i++) |
---|
457 | { |
---|
458 | _theFragments[i]->SetMomentum(_theFragments[i]->GetNuclearMass()*Vel[i]); |
---|
459 | } |
---|
460 | |
---|
461 | // garbage collection |
---|
462 | delete [] Pos; |
---|
463 | delete [] Vel; |
---|
464 | delete [] Accel; |
---|
465 | |
---|
466 | return; |
---|
467 | } |
---|
468 | |
---|
469 | |
---|
470 | |
---|
471 | G4ThreeVector G4StatMFChannel::RotateMomentum(G4ThreeVector Pa, |
---|
472 | G4ThreeVector V, G4ThreeVector P) |
---|
473 | // Rotates a 3-vector P to close momentum triangle Pa + V + P = 0 |
---|
474 | { |
---|
475 | G4ThreeVector U = Pa.unit(); |
---|
476 | |
---|
477 | G4double Alpha1 = U * V; |
---|
478 | |
---|
479 | G4double Alpha2 = std::sqrt(V.mag2() - Alpha1*Alpha1); |
---|
480 | |
---|
481 | G4ThreeVector N = (1./Alpha2)*U.cross(V); |
---|
482 | |
---|
483 | G4ThreeVector RotatedMomentum( |
---|
484 | ( (V.x() - Alpha1*U.x())/Alpha2 ) * P.x() + N.x() * P.y() + U.x() * P.z(), |
---|
485 | ( (V.y() - Alpha1*U.y())/Alpha2 ) * P.x() + N.y() * P.y() + U.y() * P.z(), |
---|
486 | ( (V.z() - Alpha1*U.z())/Alpha2 ) * P.x() + N.z() * P.y() + U.z() * P.z() |
---|
487 | ); |
---|
488 | return RotatedMomentum; |
---|
489 | } |
---|
490 | |
---|
491 | |
---|
492 | |
---|
493 | |
---|
494 | |
---|
495 | G4ThreeVector G4StatMFChannel::IsotropicVector(const G4double Magnitude) |
---|
496 | // Samples a isotropic random vector with a magnitud given by Magnitude. |
---|
497 | // By default Magnitude = 1 |
---|
498 | { |
---|
499 | G4double CosTheta = 1.0 - 2.0*G4UniformRand(); |
---|
500 | G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta); |
---|
501 | G4double Phi = twopi*G4UniformRand(); |
---|
502 | G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta, |
---|
503 | Magnitude*std::cos(Phi)*CosTheta, |
---|
504 | Magnitude*std::sin(Phi)); |
---|
505 | return Vector; |
---|
506 | } |
---|