00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include <numeric>
00038
00039 #include "G4StatMFChannel.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4HadronicException.hh"
00042 #include "G4Pow.hh"
00043
00044 class SumCoulombEnergy : public std::binary_function<G4double,G4double,G4double>
00045 {
00046 public:
00047 SumCoulombEnergy() : total(0.0) {}
00048 G4double operator() (G4double& , G4StatMFFragment*& frag)
00049 {
00050 total += frag->GetCoulombEnergy();
00051 return total;
00052 }
00053
00054 G4double GetTotal() { return total; }
00055 public:
00056 G4double total;
00057 };
00058
00059 G4StatMFChannel::G4StatMFChannel() :
00060 _NumOfNeutralFragments(0),
00061 _NumOfChargedFragments(0)
00062 {}
00063
00064 G4StatMFChannel::~G4StatMFChannel()
00065 {
00066 if (!_theFragments.empty()) {
00067 std::for_each(_theFragments.begin(),_theFragments.end(),
00068 DeleteFragment());
00069 }
00070 }
00071
00072 G4bool G4StatMFChannel::CheckFragments(void)
00073 {
00074 std::deque<G4StatMFFragment*>::iterator i;
00075 for (i = _theFragments.begin();
00076 i != _theFragments.end(); ++i)
00077 {
00078 G4int A = (*i)->GetA();
00079 G4int Z = (*i)->GetZ();
00080 if ( (A > 1 && (Z > A || Z <= 0)) || (A==1 && Z > A) || A <= 0 ) return false;
00081 }
00082 return true;
00083 }
00084
00085 void G4StatMFChannel::CreateFragment(G4int A, G4int Z)
00086
00087
00088
00089 {
00090 if (Z <= 0.5) {
00091 _theFragments.push_back(new G4StatMFFragment(A,Z));
00092 _NumOfNeutralFragments++;
00093 } else {
00094 _theFragments.push_front(new G4StatMFFragment(A,Z));
00095 _NumOfChargedFragments++;
00096 }
00097
00098 return;
00099 }
00100
00101 G4double G4StatMFChannel::GetFragmentsCoulombEnergy(void)
00102 {
00103 G4double Coulomb = std::accumulate(_theFragments.begin(),_theFragments.end(),
00104 0.0,SumCoulombEnergy());
00105
00106
00107
00108 return Coulomb;
00109 }
00110
00111
00112
00113 G4double G4StatMFChannel::GetFragmentsEnergy(G4double T) const
00114 {
00115 G4double Energy = 0.0;
00116
00117 G4double TranslationalEnergy = (3./2.)*T*static_cast<G4double>(_theFragments.size());
00118
00119 std::deque<G4StatMFFragment*>::const_iterator i;
00120 for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
00121 {
00122 Energy += (*i)->GetEnergy(T);
00123 }
00124 return Energy + TranslationalEnergy;
00125 }
00126
00127 G4FragmentVector * G4StatMFChannel::GetFragments(G4int anA,
00128 G4int anZ,
00129 G4double T)
00130
00131 {
00132
00133 CoulombImpulse(anA,anZ,T);
00134
00135
00136 FragmentsMomenta(_NumOfNeutralFragments, _NumOfChargedFragments, T);
00137
00138 G4FragmentVector * theResult = new G4FragmentVector;
00139 std::deque<G4StatMFFragment*>::iterator i;
00140 for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
00141 theResult->push_back((*i)->GetFragment(T));
00142
00143 return theResult;
00144 }
00145
00146 void G4StatMFChannel::CoulombImpulse(G4int anA, G4int anZ, G4double T)
00147
00148
00149 {
00150
00151 PlaceFragments(anA);
00152
00153
00154
00155
00156 FragmentsMomenta(_NumOfChargedFragments, 0, T);
00157
00158
00159
00160 SolveEqOfMotion(anA,anZ,T);
00161
00162 return;
00163 }
00164
00165 void G4StatMFChannel::PlaceFragments(G4int anA)
00166
00167
00168 {
00169 G4Pow* g4pow = G4Pow::GetInstance();
00170 const G4double R0 = G4StatMFParameters::Getr0();
00171 const G4double Rsys = 2.0*R0*g4pow->Z13(anA);
00172
00173 G4bool TooMuchIterations;
00174 do
00175 {
00176 TooMuchIterations = false;
00177
00178
00179 G4double R = (Rsys - R0*g4pow->Z13(_theFragments[0]->GetA()))*
00180 std::pow(G4UniformRand(),1./3.);
00181 _theFragments[0]->SetPosition(IsotropicVector(R));
00182
00183
00184
00185 G4bool ThereAreOverlaps = false;
00186 std::deque<G4StatMFFragment*>::iterator i;
00187 for (i = _theFragments.begin()+1; i != _theFragments.end(); ++i)
00188 {
00189 G4int counter = 0;
00190 do
00191 {
00192 R = (Rsys - R0*g4pow->Z13((*i)->GetA()))*std::pow(G4UniformRand(),1./3.);
00193 (*i)->SetPosition(IsotropicVector(R));
00194
00195
00196 std::deque<G4StatMFFragment*>::iterator j;
00197 for (j = _theFragments.begin(); j != i; ++j)
00198 {
00199 G4ThreeVector FragToFragVector = (*i)->GetPosition() - (*j)->GetPosition();
00200 G4double Rmin = R0*(g4pow->Z13((*i)->GetA()) +
00201 g4pow->Z13((*j)->GetA()));
00202 if ( (ThereAreOverlaps = (FragToFragVector.mag2() < Rmin*Rmin)) ) break;
00203 }
00204 counter++;
00205 } while (ThereAreOverlaps && counter < 1000);
00206
00207 if (counter >= 1000)
00208 {
00209 TooMuchIterations = true;
00210 break;
00211 }
00212 }
00213 } while (TooMuchIterations);
00214
00215 return;
00216 }
00217
00218 void G4StatMFChannel::FragmentsMomenta(G4int NF, G4int idx,
00219 G4double T)
00220
00221
00222
00223
00224
00225 {
00226 G4double KinE = (3./2.)*T*static_cast<G4double>(NF);
00227
00228 G4ThreeVector p;
00229
00230 if (NF <= 0) return;
00231 else if (NF == 1)
00232 {
00233
00234 p = IsotropicVector(std::sqrt(2.0*_theFragments[idx]->GetNuclearMass()*KinE));
00235 _theFragments[idx]->SetMomentum(p);
00236 }
00237 else if (NF == 2)
00238 {
00239
00240 G4double M1 = _theFragments[idx]->GetNuclearMass();
00241 G4double M2 = _theFragments[idx+1]->GetNuclearMass();
00242 p = IsotropicVector(std::sqrt(2.0*KinE*(M1*M2)/(M1+M2)));
00243 _theFragments[idx]->SetMomentum(p);
00244 _theFragments[idx+1]->SetMomentum(-p);
00245 }
00246 else
00247 {
00248
00249 G4double AvailableE;
00250 G4int i1,i2;
00251 G4double SummedE;
00252 G4ThreeVector SummedP;
00253 do
00254 {
00255
00256
00257 AvailableE = 0.0;
00258 SummedE = 0.0;
00259 SummedP.setX(0.0);SummedP.setY(0.0);SummedP.setZ(0.0);
00260 for (G4int i = idx; i < idx+NF-2; i++)
00261 {
00262 G4double E;
00263 G4double RandE;
00264 G4double Boltzmann;
00265 do
00266 {
00267 E = 9.0*T*G4UniformRand();
00268 Boltzmann = std::sqrt(E)*std::exp(-E/T);
00269 RandE = std::sqrt(T/2.)*std::exp(-0.5)*G4UniformRand();
00270 }
00271 while (RandE > Boltzmann);
00272 p = IsotropicVector(std::sqrt(2.0*E*_theFragments[i]->GetNuclearMass()));
00273 _theFragments[i]->SetMomentum(p);
00274 SummedE += E;
00275 SummedP += p;
00276 }
00277
00278
00279 i1 = idx+NF-2;
00280 i2 = idx+NF-1;
00281 p = -SummedP;
00282 AvailableE = KinE - SummedE;
00283
00284 }
00285 while (AvailableE <= p.mag2()/(2.0*(_theFragments[i1]->GetNuclearMass()+
00286 _theFragments[i2]->GetNuclearMass())));
00287
00288 G4double H = 1.0 + _theFragments[i2]->GetNuclearMass()/_theFragments[i1]->GetNuclearMass();
00289 G4double CTM12 = H*(1.0 - 2.0*_theFragments[i2]->GetNuclearMass()*AvailableE/p.mag2());
00290 G4double CosTheta1;
00291 G4double Sign;
00292
00293 if (CTM12 > 0.9999) {CosTheta1 = 1.;}
00294 else {
00295 do
00296 {
00297 do
00298 {
00299 CosTheta1 = 1.0 - 2.0*G4UniformRand();
00300 }
00301 while (CosTheta1*CosTheta1 < CTM12);
00302 }
00303 while (CTM12 >= 0.0 && CosTheta1 < 0.0);
00304 }
00305
00306 if (CTM12 < 0.0) Sign = 1.0;
00307 else if (G4UniformRand() <= 0.5) Sign = -1.0;
00308 else Sign = 1.0;
00309
00310
00311 G4double P1 = (p.mag()*CosTheta1+Sign*std::sqrt(p.mag2()*(CosTheta1*CosTheta1-CTM12)))/H;
00312 G4double P2 = std::sqrt(P1*P1+p.mag2() - 2.0*P1*p.mag()*CosTheta1);
00313 G4double Phi = twopi*G4UniformRand();
00314 G4double SinTheta1 = std::sqrt(1.0 - CosTheta1*CosTheta1);
00315 G4double CosPhi1 = std::cos(Phi);
00316 G4double SinPhi1 = std::sin(Phi);
00317 G4double CosPhi2 = -CosPhi1;
00318 G4double SinPhi2 = -SinPhi1;
00319 G4double CosTheta2 = (p.mag2() + P2*P2 - P1*P1)/(2.0*p.mag()*P2);
00320 G4double SinTheta2 = 0.0;
00321 if (CosTheta2 > -1.0 && CosTheta2 < 1.0) SinTheta2 = std::sqrt(1.0 - CosTheta2*CosTheta2);
00322
00323 G4ThreeVector p1(P1*SinTheta1*CosPhi1,P1*SinTheta1*SinPhi1,P1*CosTheta1);
00324 G4ThreeVector p2(P2*SinTheta2*CosPhi2,P2*SinTheta2*SinPhi2,P2*CosTheta2);
00325 G4ThreeVector b(1.0,0.0,0.0);
00326
00327 p1 = RotateMomentum(p,b,p1);
00328 p2 = RotateMomentum(p,b,p2);
00329
00330 SummedP += p1 + p2;
00331 SummedE += p1.mag2()/(2.0*_theFragments[i1]->GetNuclearMass()) +
00332 p2.mag2()/(2.0*_theFragments[i2]->GetNuclearMass());
00333
00334 _theFragments[i1]->SetMomentum(p1);
00335 _theFragments[i2]->SetMomentum(p2);
00336
00337 }
00338
00339 return;
00340 }
00341
00342
00343 void G4StatMFChannel::SolveEqOfMotion(G4int anA, G4int anZ, G4double T)
00344
00345
00346 {
00347 G4double CoulombEnergy = (3./5.)*(elm_coupling*anZ*anZ)*
00348 std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1./3.)/
00349 (G4StatMFParameters::Getr0()*G4Pow::GetInstance()->Z13(anA))
00350 - GetFragmentsCoulombEnergy();
00351 if (CoulombEnergy <= 0.0) return;
00352
00353 G4int Iterations = 0;
00354 G4double TimeN = 0.0;
00355 G4double TimeS = 0.0;
00356 G4double DeltaTime = 10.0;
00357
00358 G4ThreeVector * Pos = new G4ThreeVector[_NumOfChargedFragments];
00359 G4ThreeVector * Vel = new G4ThreeVector[_NumOfChargedFragments];
00360 G4ThreeVector * Accel = new G4ThreeVector[_NumOfChargedFragments];
00361
00362 G4int i;
00363 for (i = 0; i < _NumOfChargedFragments; i++)
00364 {
00365 Vel[i] = (1.0/(_theFragments[i]->GetNuclearMass()))*
00366 _theFragments[i]->GetMomentum();
00367 Pos[i] = _theFragments[i]->GetPosition();
00368 }
00369
00370 do
00371 {
00372
00373 G4ThreeVector distance;
00374 G4ThreeVector force;
00375
00376 for (i = 0; i < _NumOfChargedFragments; i++)
00377 {
00378 force.setX(0.0); force.setY(0.0); force.setZ(0.0);
00379 for (G4int j = 0; j < _NumOfChargedFragments; j++)
00380 {
00381 if (i != j)
00382 {
00383 distance = Pos[i] - Pos[j];
00384 force += (elm_coupling*_theFragments[i]->GetZ()
00385 *_theFragments[j]->GetZ()/
00386 (distance.mag2()*distance.mag()))*distance;
00387 }
00388 }
00389 Accel[i] = (1./(_theFragments[i]->GetNuclearMass()))*force;
00390 }
00391
00392 TimeN = TimeS + DeltaTime;
00393
00394 G4ThreeVector SavedVel;
00395 for ( i = 0; i < _NumOfChargedFragments; i++)
00396 {
00397 SavedVel = Vel[i];
00398 Vel[i] += Accel[i]*(TimeN-TimeS);
00399 Pos[i] += (SavedVel+Vel[i])*(TimeN-TimeS)*0.5;
00400 }
00401
00402
00403
00404
00405 TimeS = TimeN;
00406
00407 }
00408 while (Iterations++ < 100);
00409
00410
00411 G4double TotalKineticEnergy = 0.0;
00412 for (i = 0; i < _NumOfChargedFragments; i++)
00413 {
00414 TotalKineticEnergy += _theFragments[i]->GetNuclearMass()*
00415 0.5*Vel[i].mag2();
00416 }
00417
00418 G4double KineticEnergy = (3./2.)*_theFragments.size()*T;
00419 G4double Eta = ( CoulombEnergy + KineticEnergy ) / TotalKineticEnergy;
00420 for (i = 0; i < _NumOfChargedFragments; i++)
00421 {
00422 Vel[i] *= Eta;
00423 }
00424
00425
00426 for (i = 0; i < _NumOfChargedFragments; i++)
00427 {
00428 _theFragments[i]->SetMomentum(_theFragments[i]->GetNuclearMass()*Vel[i]);
00429 }
00430
00431
00432 delete [] Pos;
00433 delete [] Vel;
00434 delete [] Accel;
00435
00436 return;
00437 }
00438
00439
00440
00441 G4ThreeVector G4StatMFChannel::RotateMomentum(G4ThreeVector Pa,
00442 G4ThreeVector V, G4ThreeVector P)
00443
00444 {
00445 G4ThreeVector U = Pa.unit();
00446
00447 G4double Alpha1 = U * V;
00448
00449 G4double Alpha2 = std::sqrt(V.mag2() - Alpha1*Alpha1);
00450
00451 G4ThreeVector N = (1./Alpha2)*U.cross(V);
00452
00453 G4ThreeVector RotatedMomentum(
00454 ( (V.x() - Alpha1*U.x())/Alpha2 ) * P.x() + N.x() * P.y() + U.x() * P.z(),
00455 ( (V.y() - Alpha1*U.y())/Alpha2 ) * P.x() + N.y() * P.y() + U.y() * P.z(),
00456 ( (V.z() - Alpha1*U.z())/Alpha2 ) * P.x() + N.z() * P.y() + U.z() * P.z()
00457 );
00458 return RotatedMomentum;
00459 }
00460
00461
00462
00463
00464
00465 G4ThreeVector G4StatMFChannel::IsotropicVector(const G4double Magnitude)
00466
00467
00468 {
00469 G4double CosTheta = 1.0 - 2.0*G4UniformRand();
00470 G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
00471 G4double Phi = twopi*G4UniformRand();
00472 G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
00473 Magnitude*std::cos(Phi)*CosTheta,
00474 Magnitude*std::sin(Phi));
00475 return Vector;
00476 }