| Geant4 Cross Reference |
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 //
28 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara
30 //
31 // Modified:
32 // 25.07.08 I.Pshenichnov (in collaboration with Alexander Botvina and Igor
33 // Mishustin (FIAS, Frankfurt, INR, Moscow and Kurchatov Institute,
34 // Moscow, pshenich@fias.uni-frankfurt.de) fixed semi-infinite loop
35
36 #include <numeric>
37
38 #include "G4StatMFChannel.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4HadronicException.hh"
41 #include "Randomize.hh"
42 #include "G4Pow.hh"
43 #include "G4Exp.hh"
44 #include "G4RandomDirection.hh"
45
46 G4StatMFChannel::G4StatMFChannel() :
47 _NumOfNeutralFragments(0),
48 _NumOfChargedFragments(0)
49 {
50 Pos.resize(8);
51 Vel.resize(8);
52 Accel.resize(8);
53 }
54
55 G4StatMFChannel::~G4StatMFChannel()
56 {
57 if (!_theFragments.empty()) {
58 std::for_each(_theFragments.begin(),_theFragments.end(),
59 DeleteFragment());
60 }
61 }
62
63 G4bool G4StatMFChannel::CheckFragments(void)
64 {
65 std::deque<G4StatMFFragment*>::iterator i;
66 for (i = _theFragments.begin();
67 i != _theFragments.end(); ++i)
68 {
69 G4int A = (*i)->GetA();
70 G4int Z = (*i)->GetZ();
71 if ( (A > 1 && (Z > A || Z <= 0)) || (A==1 && Z > A) || A <= 0 ) return false;
72 }
73 return true;
74 }
75
76 void G4StatMFChannel::CreateFragment(G4int A, G4int Z)
77 // Create a new fragment.
78 // Fragments are automatically sorted: first charged fragments,
79 // then neutral ones.
80 {
81 if (Z <= 0.5) {
82 _theFragments.push_back(new G4StatMFFragment(A,Z));
83 _NumOfNeutralFragments++;
84 } else {
85 _theFragments.push_front(new G4StatMFFragment(A,Z));
86 _NumOfChargedFragments++;
87 }
88
89 return;
90 }
91
92 G4double G4StatMFChannel::GetFragmentsCoulombEnergy(void)
93 {
94 G4double Coulomb =
95 std::accumulate(_theFragments.begin(),_theFragments.end(),
96 0.0,
97 [](const G4double& running_total,
98 G4StatMFFragment*& fragment)
99 {
100 return running_total + fragment->GetCoulombEnergy();
101 } );
102 // G4double Coulomb = 0.0;
103 // for (unsigned int i = 0;i < _theFragments.size(); i++)
104 // Coulomb += _theFragments[i]->GetCoulombEnergy();
105 return Coulomb;
106 }
107
108 G4double G4StatMFChannel::GetFragmentsEnergy(G4double T) const
109 {
110 G4double Energy = 0.0;
111
112 G4double TranslationalEnergy = 1.5*T*_theFragments.size();
113
114 std::deque<G4StatMFFragment*>::const_iterator i;
115 for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
116 {
117 Energy += (*i)->GetEnergy(T);
118 }
119 return Energy + TranslationalEnergy;
120 }
121
122 G4FragmentVector * G4StatMFChannel::GetFragments(G4int anA,
123 G4int anZ,
124 G4double T)
125 {
126 // calculate momenta of charged fragments
127 CoulombImpulse(anA,anZ,T);
128
129 // calculate momenta of neutral fragments
130 FragmentsMomenta(_NumOfNeutralFragments, _NumOfChargedFragments, T);
131
132 G4FragmentVector * theResult = new G4FragmentVector;
133 std::deque<G4StatMFFragment*>::iterator i;
134 for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
135 theResult->push_back((*i)->GetFragment(T));
136
137 return theResult;
138 }
139
140 void G4StatMFChannel::CoulombImpulse(G4int anA, G4int anZ, G4double T)
141 // Aafter breakup, fragments fly away under Coulomb field.
142 // This method calculates asymptotic fragments momenta.
143 {
144 // First, we have to place the fragments inside of the original nucleus volume
145 PlaceFragments(anA);
146
147 // Second, we sample initial charged fragments momenta. There are
148 // _NumOfChargedFragments charged fragments and they start at the begining
149 // of the vector _theFragments (i.e. 0)
150 FragmentsMomenta(_NumOfChargedFragments, 0, T);
151
152 // Third, we have to figure out the asymptotic momenta of charged fragments
153 // For taht we have to solve equations of motion for fragments
154 SolveEqOfMotion(anA,anZ,T);
155
156 return;
157 }
158
159 void G4StatMFChannel::PlaceFragments(G4int anA)
160 // This gives the position of fragments at the breakup instant.
161 // Fragments positions are sampled inside prolongated ellipsoid.
162 {
163 G4Pow* g4calc = G4Pow::GetInstance();
164 const G4double R0 = G4StatMFParameters::Getr0();
165 G4double Rsys = 2.0*R0*g4calc->Z13(anA);
166
167 G4bool TooMuchIterations;
168 do
169 {
170 TooMuchIterations = false;
171
172 // Sample the position of the first fragment
173 G4double R = (Rsys - R0*g4calc->Z13(_theFragments[0]->GetA()))*
174 g4calc->A13(G4UniformRand());
175 _theFragments[0]->SetPosition(R*G4RandomDirection());
176
177
178 // Sample the position of the remaining fragments
179 G4bool ThereAreOverlaps = false;
180 std::deque<G4StatMFFragment*>::iterator i;
181 for (i = _theFragments.begin()+1; i != _theFragments.end(); ++i)
182 {
183 G4int counter = 0;
184 do
185 {
186 R = (Rsys - R0*g4calc->Z13((*i)->GetA()))*g4calc->A13(G4UniformRand());
187 (*i)->SetPosition(R*G4RandomDirection());
188
189 // Check that there are not overlapping fragments
190 std::deque<G4StatMFFragment*>::iterator j;
191 for (j = _theFragments.begin(); j != i; ++j)
192 {
193 G4ThreeVector FragToFragVector =
194 (*i)->GetPosition() - (*j)->GetPosition();
195 G4double Rmin = R0*(g4calc->Z13((*i)->GetA()) +
196 g4calc->Z13((*j)->GetA()));
197 if ( (ThereAreOverlaps = (FragToFragVector.mag2() < Rmin*Rmin)))
198 { break; }
199 }
200 counter++;
201 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
202 } while (ThereAreOverlaps && counter < 1000);
203
204 if (counter >= 1000)
205 {
206 TooMuchIterations = true;
207 break;
208 }
209 }
210 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
211 } while (TooMuchIterations);
212 return;
213 }
214
215 void G4StatMFChannel::FragmentsMomenta(G4int NF, G4int idx,
216 G4double T)
217 // Calculate fragments momenta at the breakup instant
218 // Fragment kinetic energies are calculated according to the
219 // Boltzmann distribution at given temperature.
220 // NF is number of fragments
221 // idx is index of first fragment
222 {
223 G4double KinE = 1.5*T*NF;
224 G4ThreeVector p(0.,0.,0.);
225
226 if (NF <= 0) return;
227 else if (NF == 1)
228 {
229 // We have only one fragment to deal with
230 p = std::sqrt(2.0*_theFragments[idx]->GetNuclearMass()*KinE)
231 *G4RandomDirection();
232 _theFragments[idx]->SetMomentum(p);
233 }
234 else if (NF == 2)
235 {
236 // We have only two fragment to deal with
237 G4double M1 = _theFragments[idx]->GetNuclearMass();
238 G4double M2 = _theFragments[idx+1]->GetNuclearMass();
239 p = std::sqrt(2.0*KinE*(M1*M2)/(M1+M2))*G4RandomDirection();
240 _theFragments[idx]->SetMomentum(p);
241 _theFragments[idx+1]->SetMomentum(-p);
242 }
243 else
244 {
245 // We have more than two fragments
246 G4double AvailableE;
247 G4int i1,i2;
248 G4double SummedE;
249 G4ThreeVector SummedP(0.,0.,0.);
250 do
251 {
252 // Fisrt sample momenta of NF-2 fragments
253 // according to Boltzmann distribution
254 AvailableE = 0.0;
255 SummedE = 0.0;
256 SummedP.setX(0.0);SummedP.setY(0.0);SummedP.setZ(0.0);
257 for (G4int i = idx; i < idx+NF-2; ++i)
258 {
259 G4double E;
260 G4double RandE;
261 do
262 {
263 E = 9.0*G4UniformRand();
264 RandE = std::sqrt(0.5/E)*G4Exp(E-0.5)*G4UniformRand();
265 }
266 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
267 while (RandE > 1.0);
268 E *= T;
269 p = std::sqrt(2.0*E*_theFragments[i]->GetNuclearMass())
270 *G4RandomDirection();
271 _theFragments[i]->SetMomentum(p);
272 SummedE += E;
273 SummedP += p;
274 }
275 // Calculate momenta of last two fragments in such a way
276 // that constraints are satisfied
277 i1 = idx+NF-2; // before last fragment index
278 i2 = idx+NF-1; // last fragment index
279 p = -SummedP;
280 AvailableE = KinE - SummedE;
281 // Available Kinetic Energy should be shared between two last fragments
282 }
283 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
284 while (AvailableE <= p.mag2()/(2.0*(_theFragments[i1]->GetNuclearMass()+
285 _theFragments[i2]->GetNuclearMass())));
286 G4double H = 1.0 + _theFragments[i2]->GetNuclearMass()
287 /_theFragments[i1]->GetNuclearMass();
288 G4double CTM12 = H*(1.0 - 2.0*_theFragments[i2]->GetNuclearMass()
289 *AvailableE/p.mag2());
290 G4double CosTheta1;
291 G4double Sign;
292
293 if (CTM12 > 1.) {CosTheta1 = 1.;}
294 else {
295 do
296 {
297 do
298 {
299 CosTheta1 = 1.0 - 2.0*G4UniformRand();
300 }
301 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
302 while (CosTheta1*CosTheta1 < CTM12);
303 }
304 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
305 while (CTM12 >= 0.0 && CosTheta1 < 0.0);
306 }
307
308 if (CTM12 < 0.0) Sign = 1.0;
309 else if (G4UniformRand() <= 0.5) Sign = -1.0;
310 else Sign = 1.0;
311
312 G4double P1 = (p.mag()*CosTheta1+Sign*std::sqrt(p.mag2()
313 *(CosTheta1*CosTheta1-CTM12)))/H;
314 G4double P2 = std::sqrt(P1*P1+p.mag2() - 2.0*P1*p.mag()*CosTheta1);
315 G4double Phi = twopi*G4UniformRand();
316 G4double SinTheta1 = std::sqrt(1.0 - CosTheta1*CosTheta1);
317 G4double CosPhi1 = std::cos(Phi);
318 G4double SinPhi1 = std::sin(Phi);
319 G4double CosPhi2 = -CosPhi1;
320 G4double SinPhi2 = -SinPhi1;
321 G4double CosTheta2 = (p.mag2() + P2*P2 - P1*P1)/(2.0*p.mag()*P2);
322 G4double SinTheta2 = 0.0;
323 if (CosTheta2 > -1.0 && CosTheta2 < 1.0) {
324 SinTheta2 = std::sqrt(1.0 - CosTheta2*CosTheta2);
325 }
326 G4ThreeVector p1(P1*SinTheta1*CosPhi1,P1*SinTheta1*SinPhi1,P1*CosTheta1);
327 G4ThreeVector p2(P2*SinTheta2*CosPhi2,P2*SinTheta2*SinPhi2,P2*CosTheta2);
328 G4ThreeVector b(1.0,0.0,0.0);
329
330 p1 = RotateMomentum(p,b,p1);
331 p2 = RotateMomentum(p,b,p2);
332
333 SummedP += p1 + p2;
334 SummedE += p1.mag2()/(2.0*_theFragments[i1]->GetNuclearMass()) +
335 p2.mag2()/(2.0*_theFragments[i2]->GetNuclearMass());
336
337 _theFragments[i1]->SetMomentum(p1);
338 _theFragments[i2]->SetMomentum(p2);
339
340 }
341 return;
342 }
343
344 void G4StatMFChannel::SolveEqOfMotion(G4int anA, G4int anZ, G4double T)
345 // This method will find a solution of Newton's equation of motion
346 // for fragments in the self-consistent time-dependent Coulomb field
347 {
348 G4Pow* g4calc = G4Pow::GetInstance();
349 G4double CoulombEnergy = 0.6*CLHEP::elm_coupling*anZ*anZ*
350 g4calc->A13(1.0+G4StatMFParameters::GetKappaCoulomb())/
351 (G4StatMFParameters::Getr0()*g4calc->Z13(anA)) - GetFragmentsCoulombEnergy();
352 if (CoulombEnergy <= 0.0) return;
353
354 G4int Iterations = 0;
355 G4double TimeN = 0.0;
356 G4double TimeS = 0.0;
357 G4double DeltaTime = 10.0;
358 if (_NumOfChargedFragments > (G4int)Pos.size()) {
359 Pos.resize(_NumOfChargedFragments);
360 Vel.resize(_NumOfChargedFragments);
361 Accel.resize(_NumOfChargedFragments);
362 }
363
364 G4int i;
365 for (i = 0; i < _NumOfChargedFragments; ++i)
366 {
367 Vel[i] = (1.0/(_theFragments[i]->GetNuclearMass()))*
368 _theFragments[i]->GetMomentum();
369 Pos[i] = _theFragments[i]->GetPosition();
370 }
371
372 G4ThreeVector distance(0.,0.,0.);
373 G4ThreeVector force(0.,0.,0.);
374 G4ThreeVector SavedVel(0.,0.,0.);
375 do {
376 for (i = 0; i < _NumOfChargedFragments; ++i)
377 {
378 force.set(0.,0.,0.);
379 for (G4int j = 0; j < _NumOfChargedFragments; ++j)
380 {
381 if (i != j)
382 {
383 distance = Pos[i] - Pos[j];
384 force += (_theFragments[i]->GetZ()*_theFragments[j]->GetZ()/
385 (distance.mag2()*distance.mag()))*distance;
386 }
387 }
388 Accel[i] = CLHEP::elm_coupling*CLHEP::fermi*force/_theFragments[i]->GetNuclearMass();
389 }
390
391 TimeN = TimeS + DeltaTime;
392
393 for ( i = 0; i < _NumOfChargedFragments; ++i)
394 {
395 SavedVel = Vel[i];
396 Vel[i] += Accel[i]*(TimeN-TimeS);
397 Pos[i] += (SavedVel+Vel[i])*(TimeN-TimeS)*0.5;
398 }
399 TimeS = TimeN;
400
401 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
402 } while (Iterations++ < 100);
403
404 // Summed fragment kinetic energy
405 G4double TotalKineticEnergy = 0.0;
406 for (i = 0; i < _NumOfChargedFragments; ++i)
407 {
408 TotalKineticEnergy += _theFragments[i]->GetNuclearMass()*
409 0.5*Vel[i].mag2();
410 }
411 // Scaling of fragment velocities
412 G4double KineticEnergy = 1.5*_NumOfChargedFragments*T;
413 G4double Eta = std::pow(( CoulombEnergy + KineticEnergy ) / TotalKineticEnergy, 2);
414
415 // Finally calculate fragments momenta
416 for (i = 0; i < _NumOfChargedFragments; ++i)
417 {
418 _theFragments[i]->SetMomentum((_theFragments[i]->GetNuclearMass()*Eta)*Vel[i]);
419 }
420 return;
421 }
422
423 G4ThreeVector G4StatMFChannel::RotateMomentum(G4ThreeVector Pa,
424 G4ThreeVector V, G4ThreeVector P)
425 // Rotates a 3-vector P to close momentum triangle Pa + V + P = 0
426 {
427 G4ThreeVector U = Pa.unit();
428
429 G4double Alpha1 = U * V;
430
431 G4double Alpha2 = std::sqrt(V.mag2() - Alpha1*Alpha1);
432
433 G4ThreeVector N = (1./Alpha2)*U.cross(V);
434
435 G4ThreeVector RotatedMomentum(
436 ( (V.x() - Alpha1*U.x())/Alpha2 ) * P.x() + N.x() * P.y() + U.x() * P.z(),
437 ( (V.y() - Alpha1*U.y())/Alpha2 ) * P.x() + N.y() * P.y() + U.y() * P.z(),
438 ( (V.z() - Alpha1*U.z())/Alpha2 ) * P.x() + N.z() * P.y() + U.z() * P.z()
439 );
440 return RotatedMomentum;
441 }
442
443