60 _NumOfNeutralFragments(0),
61 _NumOfChargedFragments(0)
66 if (!_theFragments.empty()) {
67 std::for_each(_theFragments.begin(),_theFragments.end(),
74 std::deque<G4StatMFFragment*>::iterator i;
75 for (i = _theFragments.begin();
76 i != _theFragments.end(); ++i)
78 G4int A = (*i)->GetA();
80 if ( (A > 1 && (Z > A || Z <= 0)) || (A==1 && Z > A) || A <= 0 )
return false;
92 _NumOfNeutralFragments++;
95 _NumOfChargedFragments++;
103 G4double Coulomb = std::accumulate(_theFragments.begin(),_theFragments.end(),
117 G4double TranslationalEnergy = (3./2.)*T*
static_cast<G4double>(_theFragments.size());
119 std::deque<G4StatMFFragment*>::const_iterator i;
120 for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
122 Energy += (*i)->GetEnergy(T);
124 return Energy + TranslationalEnergy;
133 CoulombImpulse(anA,anZ,T);
136 FragmentsMomenta(_NumOfNeutralFragments, _NumOfChargedFragments, T);
139 std::deque<G4StatMFFragment*>::iterator i;
140 for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
141 theResult->push_back((*i)->GetFragment(T));
156 FragmentsMomenta(_NumOfChargedFragments, 0, T);
160 SolveEqOfMotion(anA,anZ,T);
165 void G4StatMFChannel::PlaceFragments(
G4int anA)
176 TooMuchIterations =
false;
179 G4double R = (Rsys - R0*g4pow->
Z13(_theFragments[0]->GetA()))*
181 _theFragments[0]->SetPosition(IsotropicVector(R));
185 G4bool ThereAreOverlaps =
false;
186 std::deque<G4StatMFFragment*>::iterator i;
187 for (i = _theFragments.begin()+1; i != _theFragments.end(); ++i)
193 (*i)->SetPosition(IsotropicVector(R));
196 std::deque<G4StatMFFragment*>::iterator j;
197 for (j = _theFragments.begin(); j != i; ++j)
199 G4ThreeVector FragToFragVector = (*i)->GetPosition() - (*j)->GetPosition();
201 g4pow->
Z13((*j)->GetA()));
202 if ( (ThereAreOverlaps = (FragToFragVector.
mag2() < Rmin*Rmin)) )
break;
205 }
while (ThereAreOverlaps && counter < 1000);
209 TooMuchIterations =
true;
213 }
while (TooMuchIterations);
218 void G4StatMFChannel::FragmentsMomenta(
G4int NF,
G4int idx,
234 p = IsotropicVector(std::sqrt(2.0*_theFragments[idx]->GetNuclearMass()*KinE));
235 _theFragments[idx]->SetMomentum(p);
240 G4double M1 = _theFragments[idx]->GetNuclearMass();
241 G4double M2 = _theFragments[idx+1]->GetNuclearMass();
242 p = IsotropicVector(std::sqrt(2.0*KinE*(M1*M2)/(M1+M2)));
243 _theFragments[idx]->SetMomentum(p);
244 _theFragments[idx+1]->SetMomentum(-p);
260 for (
G4int i = idx; i < idx+NF-2; i++)
268 Boltzmann = std::sqrt(E)*std::exp(-E/T);
271 while (RandE > Boltzmann);
272 p = IsotropicVector(std::sqrt(2.0*E*_theFragments[i]->GetNuclearMass()));
273 _theFragments[i]->SetMomentum(p);
282 AvailableE = KinE - SummedE;
285 while (AvailableE <= p.mag2()/(2.0*(_theFragments[i1]->GetNuclearMass()+
286 _theFragments[i2]->GetNuclearMass())));
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());
293 if (CTM12 > 0.9999) {CosTheta1 = 1.;}
301 while (CosTheta1*CosTheta1 < CTM12);
303 while (CTM12 >= 0.0 && CosTheta1 < 0.0);
306 if (CTM12 < 0.0) Sign = 1.0;
311 G4double P1 = (p.mag()*CosTheta1+Sign*std::sqrt(p.mag2()*(CosTheta1*CosTheta1-CTM12)))/H;
312 G4double P2 = std::sqrt(P1*P1+p.mag2() - 2.0*P1*p.mag()*CosTheta1);
314 G4double SinTheta1 = std::sqrt(1.0 - CosTheta1*CosTheta1);
319 G4double CosTheta2 = (p.mag2() + P2*P2 - P1*P1)/(2.0*p.mag()*P2);
321 if (CosTheta2 > -1.0 && CosTheta2 < 1.0) SinTheta2 = std::sqrt(1.0 - CosTheta2*CosTheta2);
323 G4ThreeVector p1(P1*SinTheta1*CosPhi1,P1*SinTheta1*SinPhi1,P1*CosTheta1);
324 G4ThreeVector p2(P2*SinTheta2*CosPhi2,P2*SinTheta2*SinPhi2,P2*CosTheta2);
327 p1 = RotateMomentum(p,
b,p1);
328 p2 = RotateMomentum(p,
b,p2);
331 SummedE += p1.
mag2()/(2.0*_theFragments[i1]->GetNuclearMass()) +
332 p2.mag2()/(2.0*_theFragments[i2]->GetNuclearMass());
334 _theFragments[i1]->SetMomentum(p1);
335 _theFragments[i2]->SetMomentum(p2);
351 if (CoulombEnergy <= 0.0)
return;
353 G4int Iterations = 0;
363 for (i = 0; i < _NumOfChargedFragments; i++)
365 Vel[i] = (1.0/(_theFragments[i]->GetNuclearMass()))*
366 _theFragments[i]->GetMomentum();
367 Pos[i] = _theFragments[i]->GetPosition();
376 for (i = 0; i < _NumOfChargedFragments; i++)
379 for (
G4int j = 0; j < _NumOfChargedFragments; j++)
383 distance = Pos[i] - Pos[j];
385 *_theFragments[j]->GetZ()/
386 (distance.
mag2()*distance.
mag()))*distance;
389 Accel[i] = (1./(_theFragments[i]->GetNuclearMass()))*force;
392 TimeN = TimeS + DeltaTime;
395 for ( i = 0; i < _NumOfChargedFragments; i++)
398 Vel[i] += Accel[i]*(TimeN-TimeS);
399 Pos[i] += (SavedVel+Vel[i])*(TimeN-TimeS)*0.5;
408 while (Iterations++ < 100);
412 for (i = 0; i < _NumOfChargedFragments; i++)
414 TotalKineticEnergy += _theFragments[i]->GetNuclearMass()*
418 G4double KineticEnergy = (3./2.)*_theFragments.size()*T;
419 G4double Eta = ( CoulombEnergy + KineticEnergy ) / TotalKineticEnergy;
420 for (i = 0; i < _NumOfChargedFragments; i++)
426 for (i = 0; i < _NumOfChargedFragments; i++)
428 _theFragments[i]->SetMomentum(_theFragments[i]->GetNuclearMass()*Vel[i]);
449 G4double Alpha2 = std::sqrt(V.mag2() - Alpha1*Alpha1);
454 ( (V.x() - Alpha1*U.
x())/Alpha2 ) * P.
x() + N.x() * P.
y() + U.
x() * P.
z(),
455 ( (V.y() - Alpha1*U.
y())/Alpha2 ) * P.
x() + N.y() * P.
y() + U.
y() * P.
z(),
456 ( (V.z() - Alpha1*U.
z())/Alpha2 ) * P.
x() + N.z() * P.
y() + U.
z() * P.
z()
458 return RotatedMomentum;
470 G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
473 Magnitude*std::cos(Phi)*CosTheta,
474 Magnitude*std::sin(Phi));
static G4Pow * GetInstance()
static G4double GetKappaCoulomb()
G4double GetFragmentsEnergy(G4double T) const
G4double GetCoulombEnergy(void) const
G4double operator()(G4double &, G4StatMFFragment *&frag)
void CreateFragment(G4int A, G4int Z)
G4double Z13(G4int Z) const
std::vector< G4Fragment * > G4FragmentVector
G4double GetFragmentsCoulombEnergy(void)
Hep3Vector cross(const Hep3Vector &) const
short Sign(short a, short b)
G4bool CheckFragments(void)
G4FragmentVector * GetFragments(G4int anA, G4int anZ, G4double T)