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