| 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,v 1.10 2008-11-19 14:33:31 vnivanch Exp $
>> 28 // GEANT4 tag $Name: geant4-09-04-patch-02 $
27 // 29 //
28 // Hadronic Process: Nuclear De-excitations 30 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara 31 // by V. Lara
30 // 32 //
31 // Modified: 33 // Modified:
32 // 25.07.08 I.Pshenichnov (in collaboration wi 34 // 25.07.08 I.Pshenichnov (in collaboration with Alexander Botvina and Igor
33 // Mishustin (FIAS, Frankfurt, INR, M 35 // Mishustin (FIAS, Frankfurt, INR, Moscow and Kurchatov Institute,
34 // Moscow, pshenich@fias.uni-frankfur 36 // Moscow, pshenich@fias.uni-frankfurt.de) fixed semi-infinite loop
35 37
36 #include <numeric> <<
37 38
38 #include "G4StatMFChannel.hh" 39 #include "G4StatMFChannel.hh"
39 #include "G4PhysicalConstants.hh" <<
40 #include "G4HadronicException.hh" 40 #include "G4HadronicException.hh"
41 #include "Randomize.hh" << 41 #include <numeric>
42 #include "G4Pow.hh" << 42
43 #include "G4Exp.hh" << 43 class SumCoulombEnergy : public std::binary_function<G4double,G4double,G4double>
44 #include "G4RandomDirection.hh" << 44 {
45 << 45 public:
46 G4StatMFChannel::G4StatMFChannel() : << 46 SumCoulombEnergy() : total(0.0) {}
47 _NumOfNeutralFragments(0), << 47 G4double operator() (G4double& , G4StatMFFragment*& frag)
48 _NumOfChargedFragments(0) << 48 {
49 { << 49 total += frag->GetCoulombEnergy();
50 Pos.resize(8); << 50 return total;
51 Vel.resize(8); << 51 }
52 Accel.resize(8); << 52
53 } << 53 G4double GetTotal() { return total; }
54 << 54 public:
55 G4StatMFChannel::~G4StatMFChannel() << 55 G4double total;
56 { << 56 };
57 if (!_theFragments.empty()) { << 57
58 std::for_each(_theFragments.begin(),_theFr << 58
59 DeleteFragment()); << 59
60 } << 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;
61 } 75 }
62 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
63 G4bool G4StatMFChannel::CheckFragments(void) 92 G4bool G4StatMFChannel::CheckFragments(void)
64 { 93 {
65 std::deque<G4StatMFFragment*>::iterator i; << 94 std::deque<G4StatMFFragment*>::iterator i;
66 for (i = _theFragments.begin(); << 95 for (i = _theFragments.begin();
67 i != _theFragments.end(); ++i) << 96 i != _theFragments.end(); ++i)
68 { << 97 {
69 G4int A = (*i)->GetA(); << 98 G4int A = static_cast<G4int>((*i)->GetA());
70 G4int Z = (*i)->GetZ(); << 99 G4int Z = static_cast<G4int>((*i)->GetZ());
71 if ( (A > 1 && (Z > A || Z <= 0)) || (A= << 100 if ( (A > 1 && (Z > A || Z <= 0)) || (A==1 && Z > A) || A <= 0 ) return false;
72 } 101 }
>> 102
73 return true; 103 return true;
74 } 104 }
75 105
76 void G4StatMFChannel::CreateFragment(G4int A, << 106
77 // Create a new fragment. << 107
78 // Fragments are automatically sorted: first c << 108
79 // then neutral ones. << 109 void G4StatMFChannel::CreateFragment(const G4double A, const G4double Z)
80 { << 110 // Create a new fragment.
81 if (Z <= 0.5) { << 111 // Fragments are automatically sorted: first charged fragments,
82 _theFragments.push_back(new G4StatMFFragme << 112 // then neutral ones.
83 _NumOfNeutralFragments++; << 113 {
84 } else { << 114 if (Z <= 0.5) {
85 _theFragments.push_front(new G4StatMFFragm << 115 _theFragments.push_back(new G4StatMFFragment(static_cast<G4int>(A),static_cast<G4int>(Z)));
86 _NumOfChargedFragments++; << 116 _NumOfNeutralFragments++;
87 } << 117 } else {
>> 118 _theFragments.push_front(new G4StatMFFragment(static_cast<G4int>(A),static_cast<G4int>(Z)));
>> 119 _NumOfChargedFragments++;
>> 120 }
88 121
89 return; << 122 return;
90 } 123 }
91 124
>> 125
92 G4double G4StatMFChannel::GetFragmentsCoulombE 126 G4double G4StatMFChannel::GetFragmentsCoulombEnergy(void)
93 { 127 {
94 G4double Coulomb = << 128 G4double Coulomb = std::accumulate(_theFragments.begin(),_theFragments.end(),
95 std::accumulate(_theFragments.begin(),_the << 129 0.0,SumCoulombEnergy());
96 0.0, << 130 // G4double Coulomb = 0.0;
97 [](const G4double& running << 131 // for (unsigned int i = 0;i < _theFragments.size(); i++)
98 G4StatMFFragment*& frag << 132 // Coulomb += _theFragments[i]->GetCoulombEnergy();
99 { << 133 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 } 134 }
107 135
108 G4double G4StatMFChannel::GetFragmentsEnergy(G << 136
>> 137
>> 138 G4double G4StatMFChannel::GetFragmentsEnergy(const G4double T) const
109 { 139 {
110 G4double Energy = 0.0; << 140 G4double Energy = 0.0;
111 141
112 G4double TranslationalEnergy = 1.5*T*_theFra << 142 G4double TranslationalEnergy = (3./2.)*T*static_cast<G4double>(_theFragments.size());
113 143
114 std::deque<G4StatMFFragment*>::const_iterato << 144 std::deque<G4StatMFFragment*>::const_iterator i;
115 for (i = _theFragments.begin(); i != _theFra << 145 for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
116 { << 146 {
117 Energy += (*i)->GetEnergy(T); << 147 Energy += (*i)->GetEnergy(T);
118 } << 148 }
119 return Energy + TranslationalEnergy; << 149 return Energy + TranslationalEnergy;
120 } 150 }
121 151
122 G4FragmentVector * G4StatMFChannel::GetFragmen << 152 G4FragmentVector * G4StatMFChannel::GetFragments(const G4double anA,
123 G4int anZ, << 153 const G4double anZ,
124 G4double T) << 154 const G4double T)
>> 155 //
125 { 156 {
126 // calculate momenta of charged fragments << 157 // calculate momenta of charged fragments
127 CoulombImpulse(anA,anZ,T); << 158 CoulombImpulse(anA,anZ,T);
128 159
129 // calculate momenta of neutral fragments << 160 // calculate momenta of neutral fragments
130 FragmentsMomenta(_NumOfNeutralFragments, _Nu << 161 FragmentsMomenta(_NumOfNeutralFragments, _NumOfChargedFragments, T);
131 162
132 G4FragmentVector * theResult = new G4Fragmen <<
133 std::deque<G4StatMFFragment*>::iterator i; <<
134 for (i = _theFragments.begin(); i != _theFra <<
135 theResult->push_back((*i)->GetFragment(T)) <<
136 163
137 return theResult; << 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
138 } 171 }
139 172
140 void G4StatMFChannel::CoulombImpulse(G4int anA << 173
141 // Aafter breakup, fragments fly away under Co << 174
142 // This method calculates asymptotic fragments << 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.
143 { 178 {
144 // First, we have to place the fragments ins << 179 // First, we have to place the fragments inside of the original nucleus volume
145 PlaceFragments(anA); << 180 PlaceFragments(anA);
146 181
147 // Second, we sample initial charged fragmen << 182 // Second, we sample initial charged fragments momenta. There are
148 // _NumOfChargedFragments charged fragments << 183 // _NumOfChargedFragments charged fragments and they start at the begining
149 // of the vector _theFragments (i.e. 0) << 184 // of the vector _theFragments (i.e. 0)
150 FragmentsMomenta(_NumOfChargedFragments, 0, << 185 FragmentsMomenta(_NumOfChargedFragments, 0, T);
151 186
152 // Third, we have to figure out the asymptot << 187 // Third, we have to figure out the asymptotic momenta of charged fragments
153 // For taht we have to solve equations of mo << 188 // For taht we have to solve equations of motion for fragments
154 SolveEqOfMotion(anA,anZ,T); << 189 SolveEqOfMotion(anA,anZ,T);
155 190
156 return; << 191 return;
157 } 192 }
158 193
159 void G4StatMFChannel::PlaceFragments(G4int anA <<
160 // This gives the position of fragments at the <<
161 // Fragments positions are sampled inside prol <<
162 { <<
163 G4Pow* g4calc = G4Pow::GetInstance(); <<
164 const G4double R0 = G4StatMFParameters::Getr <<
165 G4double Rsys = 2.0*R0*g4calc->Z13(anA); <<
166 194
167 G4bool TooMuchIterations; << 195
168 do << 196 void G4StatMFChannel::PlaceFragments(const G4double anA)
169 { << 197 // This gives the position of fragments at the breakup instant.
170 TooMuchIterations = false; << 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;
171 207
172 // Sample the position of the first frag << 208 // Sample the position of the first fragment
173 G4double R = (Rsys - R0*g4calc->Z13(_the << 209 G4double R = (Rsys - R0*std::pow(_theFragments[0]->GetA(),1./3.))*
174 g4calc->A13(G4UniformRand()); << 210 std::pow(G4UniformRand(),1./3.);
175 _theFragments[0]->SetPosition(R*G4Random << 211 _theFragments[0]->SetPosition(IsotropicVector(R));
176 212
177 213
178 // Sample the position of the remaining << 214 // Sample the position of the remaining fragments
179 G4bool ThereAreOverlaps = false; << 215 G4bool ThereAreOverlaps = false;
180 std::deque<G4StatMFFragment*>::iterator << 216 std::deque<G4StatMFFragment*>::iterator i;
181 for (i = _theFragments.begin()+1; i != _ << 217 for (i = _theFragments.begin()+1; i != _theFragments.end(); ++i)
182 { << 218 {
183 G4int counter = 0; << 219 G4int counter = 0;
184 do << 220 do
185 { << 221 {
186 R = (Rsys - R0*g4calc->Z13((*i)->GetA( << 222 R = (Rsys - R0*std::pow((*i)->GetA(),1./3.))*std::pow(G4UniformRand(),1./3.);
187 (*i)->SetPosition(R*G4RandomDirection( << 223 (*i)->SetPosition(IsotropicVector(R));
188 224
189 // Check that there are not overlappin << 225 // Check that there are not overlapping fragments
190 std::deque<G4StatMFFragment*>::iterato << 226 std::deque<G4StatMFFragment*>::iterator j;
191 for (j = _theFragments.begin(); j != i << 227 for (j = _theFragments.begin(); j != i; ++j)
192 { << 228 {
193 G4ThreeVector FragToFragVector = << 229 G4ThreeVector FragToFragVector = (*i)->GetPosition() - (*j)->GetPosition();
194 (*i)->GetPosition() - (*j)->GetPositio << 230 G4double Rmin = R0*(std::pow((*i)->GetA(),1./3.) +
195 G4double Rmin = R0*(g4calc->Z13((*i)->Ge << 231 std::pow((*j)->GetA(),1./3));
196 g4calc->Z13((*j)->GetA())); << 232 if ( (ThereAreOverlaps = (FragToFragVector.mag2() < Rmin*Rmin)) ) break;
197 if ( (ThereAreOverlaps = (FragToFragVect << 233 }
198 { break; } << 234 counter++;
199 } << 235 } while (ThereAreOverlaps && counter < 1000);
200 counter++; <<
201 // Loop checking, 05-Aug-2015, Vladimi <<
202 } while (ThereAreOverlaps && counter < 1 <<
203 236
204 if (counter >= 1000) << 237 if (counter >= 1000)
205 { << 238 {
206 TooMuchIterations = true; << 239 TooMuchIterations = true;
207 break; << 240 break;
208 } << 241 }
209 } << 242 }
210 // Loop checking, 05-Aug-2015, Vladimir <<
211 } while (TooMuchIterations); 243 } while (TooMuchIterations);
>> 244
212 return; 245 return;
213 } 246 }
214 247
215 void G4StatMFChannel::FragmentsMomenta(G4int N << 248
216 G4double T) << 249 void G4StatMFChannel::FragmentsMomenta(const G4int NF, const G4int idx,
217 // Calculate fragments momenta at the breakup << 250 const G4double T)
218 // Fragment kinetic energies are calculated ac << 251 // Calculate fragments momenta at the breakup instant
219 // Boltzmann distribution at given temperature << 252 // Fragment kinetic energies are calculated according to the
220 // NF is number of fragments << 253 // Boltzmann distribution at given temperature.
221 // idx is index of first fragment << 254 // NF is number of fragments
>> 255 // idx is index of first fragment
222 { 256 {
223 G4double KinE = 1.5*T*NF; << 257 G4double KinE = (3./2.)*T*static_cast<G4double>(NF);
224 G4ThreeVector p(0.,0.,0.); <<
225 258
226 if (NF <= 0) return; << 259 G4ThreeVector p;
227 else if (NF == 1) << 260
228 { << 261 if (NF <= 0) return;
229 // We have only one fragment to deal wit << 262 else if (NF == 1)
230 p = std::sqrt(2.0*_theFragments[idx]->Ge << 263 {
231 *G4RandomDirection(); << 264 // We have only one fragment to deal with
232 _theFragments[idx]->SetMomentum(p); << 265 p = IsotropicVector(std::sqrt(2.0*_theFragments[idx]->GetNuclearMass()*KinE));
233 } << 266 _theFragments[idx]->SetMomentum(p);
234 else if (NF == 2) << 267 }
235 { << 268 else if (NF == 2)
236 // We have only two fragment to deal wit << 269 {
237 G4double M1 = _theFragments[idx]->GetNuc << 270 // We have only two fragment to deal with
238 G4double M2 = _theFragments[idx+1]->GetN << 271 G4double M1 = _theFragments[idx]->GetNuclearMass();
239 p = std::sqrt(2.0*KinE*(M1*M2)/(M1+M2))* << 272 G4double M2 = _theFragments[idx+1]->GetNuclearMass();
240 _theFragments[idx]->SetMomentum(p); << 273 p = IsotropicVector(std::sqrt(2.0*KinE*(M1*M2)/(M1+M2)));
241 _theFragments[idx+1]->SetMomentum(-p); << 274 _theFragments[idx]->SetMomentum(p);
242 } << 275 _theFragments[idx+1]->SetMomentum(-p);
243 else << 276 }
244 { << 277 else
245 // We have more than two fragments << 278 {
246 G4double AvailableE; << 279 // We have more than two fragments
247 G4int i1,i2; << 280 G4double AvailableE;
248 G4double SummedE; << 281 G4int i1,i2;
249 G4ThreeVector SummedP(0.,0.,0.); << 282 G4double SummedE;
250 do << 283 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 284 do
296 { 285 {
297 do << 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++)
298 { 292 {
299 CosTheta1 = 1.0 - 2.0*G4UniformRand(); << 293 G4double E;
300 } << 294 G4double RandE;
301 // Loop checking, 05-Aug-2015, Vladimir << 295 G4double Boltzmann;
302 while (CosTheta1*CosTheta1 < CTM12); << 296 do
303 } << 297 {
304 // Loop checking, 05-Aug-2015, Vladimir Ivan << 298 E = 9.0*T*G4UniformRand();
305 while (CTM12 >= 0.0 && CosTheta1 < 0.0); << 299 Boltzmann = std::sqrt(E)*std::exp(-E/T);
306 } << 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 }
307 336
308 if (CTM12 < 0.0) Sign = 1.0; << 337 if (CTM12 < 0.0) Sign = 1.0;
309 else if (G4UniformRand() <= 0.5) Sign = << 338 else if (G4UniformRand() <= 0.5) Sign = -1.0;
310 else Sign = 1.0; << 339 else Sign = 1.0;
311 340
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 341
337 _theFragments[i1]->SetMomentum(p1); << 342 G4double P1 = (p.mag()*CosTheta1+Sign*std::sqrt(p.mag2()*(CosTheta1*CosTheta1-CTM12)))/H;
338 _theFragments[i2]->SetMomentum(p2); << 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());
339 364
340 } << 365 _theFragments[i1]->SetMomentum(p1);
341 return; << 366 _theFragments[i2]->SetMomentum(p2);
>> 367
>> 368 }
>> 369
>> 370 return;
342 } 371 }
343 372
344 void G4StatMFChannel::SolveEqOfMotion(G4int an << 373
345 // This method will find a solution of Newton' << 374 void G4StatMFChannel::SolveEqOfMotion(const G4double anA, const G4double anZ, const G4double T)
346 // for fragments in the self-consistent time-d << 375 // This method will find a solution of Newton's equation of motion
347 { << 376 // for fragments in the self-consistent time-dependent Coulomb field
348 G4Pow* g4calc = G4Pow::GetInstance(); << 377 {
349 G4double CoulombEnergy = 0.6*CLHEP::elm_coup << 378 G4double CoulombEnergy = (3./5.)*(elm_coupling*anZ*anZ)*
350 g4calc->A13(1.0+G4StatMFParameters::GetKap << 379 std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1./3.)/
351 (G4StatMFParameters::Getr0()*g4calc->Z13(a << 380 (G4StatMFParameters::Getr0()*std::pow(anA,1./3.))
>> 381 - GetFragmentsCoulombEnergy();
352 if (CoulombEnergy <= 0.0) return; 382 if (CoulombEnergy <= 0.0) return;
353 383
354 G4int Iterations = 0; 384 G4int Iterations = 0;
355 G4double TimeN = 0.0; 385 G4double TimeN = 0.0;
356 G4double TimeS = 0.0; 386 G4double TimeS = 0.0;
357 G4double DeltaTime = 10.0; 387 G4double DeltaTime = 10.0;
358 if (_NumOfChargedFragments > (G4int)Pos.size << 388
359 Pos.resize(_NumOfChargedFragments); << 389 G4ThreeVector * Pos = new G4ThreeVector[_NumOfChargedFragments];
360 Vel.resize(_NumOfChargedFragments); << 390 G4ThreeVector * Vel = new G4ThreeVector[_NumOfChargedFragments];
361 Accel.resize(_NumOfChargedFragments); << 391 G4ThreeVector * Accel = new G4ThreeVector[_NumOfChargedFragments];
362 } <<
363 392
364 G4int i; 393 G4int i;
365 for (i = 0; i < _NumOfChargedFragments; ++i) << 394 for (i = 0; i < _NumOfChargedFragments; i++)
366 { 395 {
367 Vel[i] = (1.0/(_theFragments[i]->GetNucl 396 Vel[i] = (1.0/(_theFragments[i]->GetNuclearMass()))*
368 _theFragments[i]->GetMomentum(); 397 _theFragments[i]->GetMomentum();
369 Pos[i] = _theFragments[i]->GetPosition() 398 Pos[i] = _theFragments[i]->GetPosition();
370 } 399 }
>> 400
>> 401 do
>> 402 {
371 403
372 G4ThreeVector distance(0.,0.,0.); << 404 G4ThreeVector distance;
373 G4ThreeVector force(0.,0.,0.); << 405 G4ThreeVector force;
374 G4ThreeVector SavedVel(0.,0.,0.); << 406
375 do { << 407 for (i = 0; i < _NumOfChargedFragments; i++)
376 for (i = 0; i < _NumOfChargedFragments; ++ << 408 {
377 { << 409 force.setX(0.0); force.setY(0.0); force.setZ(0.0);
378 force.set(0.,0.,0.); << 410 for (G4int j = 0; j < _NumOfChargedFragments; j++)
379 for (G4int j = 0; j < _NumOfChargedFragments << 411 {
380 { << 412 if (i != j)
381 if (i != j) << 413 {
382 { << 414 distance = Pos[i] - Pos[j];
383 distance = Pos[i] - Pos[j]; << 415 force += (elm_coupling*(_theFragments[i]->GetZ()*_theFragments[j]->GetZ())/
384 force += (_theFragments[i]->GetZ()*_theFra << 416 (distance.mag2()*distance.mag()))*distance;
385 (distance.mag2()*distance.mag()))*dist << 417 }
386 } << 418 }
387 } << 419 Accel[i] = (1./(_theFragments[i]->GetNuclearMass()))*force;
388 Accel[i] = CLHEP::elm_coupling*CLHEP::fermi* << 420 }
389 } <<
390 421
391 TimeN = TimeS + DeltaTime; << 422 TimeN = TimeS + DeltaTime;
392 423
393 for ( i = 0; i < _NumOfChargedFragments; + << 424 G4ThreeVector SavedVel;
394 { << 425 for ( i = 0; i < _NumOfChargedFragments; i++)
395 SavedVel = Vel[i]; << 426 {
396 Vel[i] += Accel[i]*(TimeN-TimeS); << 427 SavedVel = Vel[i];
397 Pos[i] += (SavedVel+Vel[i])*(TimeN-TimeS)*0. << 428 Vel[i] += Accel[i]*(TimeN-TimeS);
398 } << 429 Pos[i] += (SavedVel+Vel[i])*(TimeN-TimeS)*0.5;
399 TimeS = TimeN; << 430 }
400 << 431
401 // Loop checking, 05-Aug-2015, Vladimir Iv << 432 // if (Iterations >= 50 && Iterations < 75) DeltaTime = 4.;
402 } while (Iterations++ < 100); << 433 // else if (Iterations >= 75) DeltaTime = 10.;
>> 434
>> 435 TimeS = TimeN;
>> 436
>> 437 }
>> 438 while (Iterations++ < 100);
403 439
404 // Summed fragment kinetic energy 440 // Summed fragment kinetic energy
405 G4double TotalKineticEnergy = 0.0; 441 G4double TotalKineticEnergy = 0.0;
406 for (i = 0; i < _NumOfChargedFragments; ++i) << 442 for (i = 0; i < _NumOfChargedFragments; i++)
407 { 443 {
408 TotalKineticEnergy += _theFragments[i]-> 444 TotalKineticEnergy += _theFragments[i]->GetNuclearMass()*
409 0.5*Vel[i].mag2(); 445 0.5*Vel[i].mag2();
410 } 446 }
411 // Scaling of fragment velocities 447 // Scaling of fragment velocities
412 G4double KineticEnergy = 1.5*_NumOfChargedFr << 448 G4double KineticEnergy = (3./2.)*static_cast<G4double>(_theFragments.size())*T;
413 G4double Eta = std::pow(( CoulombEnergy + Ki << 449 G4double Eta = ( CoulombEnergy + KineticEnergy ) / TotalKineticEnergy;
414 << 450 for (i = 0; i < _NumOfChargedFragments; i++)
>> 451 {
>> 452 Vel[i] *= Eta;
>> 453 }
>> 454
415 // Finally calculate fragments momenta 455 // Finally calculate fragments momenta
416 for (i = 0; i < _NumOfChargedFragments; ++i) << 456 for (i = 0; i < _NumOfChargedFragments; i++)
417 { 457 {
418 _theFragments[i]->SetMomentum((_theFragm << 458 _theFragments[i]->SetMomentum(_theFragments[i]->GetNuclearMass()*Vel[i]);
419 } 459 }
>> 460
>> 461 // garbage collection
>> 462 delete [] Pos;
>> 463 delete [] Vel;
>> 464 delete [] Accel;
>> 465
420 return; 466 return;
421 } 467 }
422 468
>> 469
>> 470
423 G4ThreeVector G4StatMFChannel::RotateMomentum( 471 G4ThreeVector G4StatMFChannel::RotateMomentum(G4ThreeVector Pa,
424 G4ThreeVector V, G4ThreeVector 472 G4ThreeVector V, G4ThreeVector P)
425 // Rotates a 3-vector P to close momentum 473 // Rotates a 3-vector P to close momentum triangle Pa + V + P = 0
426 { 474 {
427 G4ThreeVector U = Pa.unit(); 475 G4ThreeVector U = Pa.unit();
428 476
429 G4double Alpha1 = U * V; 477 G4double Alpha1 = U * V;
430 478
431 G4double Alpha2 = std::sqrt(V.mag2() - Alpha 479 G4double Alpha2 = std::sqrt(V.mag2() - Alpha1*Alpha1);
432 480
433 G4ThreeVector N = (1./Alpha2)*U.cross(V); 481 G4ThreeVector N = (1./Alpha2)*U.cross(V);
434 482
435 G4ThreeVector RotatedMomentum( 483 G4ThreeVector RotatedMomentum(
436 ( (V.x() - Alpha1*U.x())/Alpha2 ) * P. 484 ( (V.x() - Alpha1*U.x())/Alpha2 ) * P.x() + N.x() * P.y() + U.x() * P.z(),
437 ( (V.y() - Alpha1*U.y())/Alpha2 ) * P. 485 ( (V.y() - Alpha1*U.y())/Alpha2 ) * P.x() + N.y() * P.y() + U.y() * P.z(),
438 ( (V.z() - Alpha1*U.z())/Alpha2 ) * P. 486 ( (V.z() - Alpha1*U.z())/Alpha2 ) * P.x() + N.z() * P.y() + U.z() * P.z()
439 ); 487 );
440 return RotatedMomentum; 488 return RotatedMomentum;
441 } 489 }
442 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 }
443 507