G4StatMFChannel Class Reference

#include <G4StatMFChannel.hh>

Collaboration diagram for G4StatMFChannel:

Classes
struct	DeleteFragment

Public Member Functions
	G4StatMFChannel ()
	~G4StatMFChannel ()
void	CreateFragment (G4int A, G4int Z)
size_t	GetMultiplicity (void)
G4bool	CheckFragments (void)
G4double	GetFragmentsCoulombEnergy (void)
G4double	GetFragmentsEnergy (G4double T) const
G4FragmentVector *	GetFragments (G4int anA, G4int anZ, G4double T)

Private Member Functions
	G4StatMFChannel (const G4StatMFChannel &right)
G4StatMFChannel &	operator= (const G4StatMFChannel &right)
G4bool	operator== (const G4StatMFChannel &right) const
G4bool	operator!= (const G4StatMFChannel &right) const
void	CoulombImpulse (G4int anA, G4int anZ, G4double T)
void	PlaceFragments (G4int anA)
void	SolveEqOfMotion (G4int anA, G4int anZ, G4double T)
void	FragmentsMomenta (G4int NF, G4int idx, G4double T)
G4ThreeVector	IsotropicVector (G4double Magnitude=1.0)
G4ThreeVector	RotateMomentum (G4ThreeVector Pa, G4ThreeVector V, G4ThreeVector P)

Private Attributes
std::deque< G4StatMFFragment * >	_theFragments
G4int	_NumOfNeutralFragments
G4int	_NumOfChargedFragments

Detailed Description

Definition at line 41 of file G4StatMFChannel.hh.

Constructor & Destructor Documentation

G4StatMFChannel() [1/2]

G4StatMFChannel::G4StatMFChannel

(

)

Definition at line 61 of file G4StatMFChannel.cc.

                                 : 
  _NumOfNeutralFragments(0), 
  _NumOfChargedFragments(0)
{}

~G4StatMFChannel()

G4StatMFChannel::~G4StatMFChannel

(

)

Definition at line 66 of file G4StatMFChannel.cc.

{ 
  if (!_theFragments.empty()) {
    std::for_each(_theFragments.begin(),_theFragments.end(),
          DeleteFragment());
  }
}

G4StatMFChannel() [2/2]

G4StatMFChannel::G4StatMFChannel ( const G4StatMFChannel &  right )

private

Member Function Documentation

CheckFragments()

G4bool G4StatMFChannel::CheckFragments

(

void

)

Definition at line 74 of file G4StatMFChannel.cc.

{
  std::deque<G4StatMFFragment*>::iterator i;
  for (i = _theFragments.begin(); 
       i != _theFragments.end(); ++i) 
    {
      G4int A = (*i)->GetA();
      G4int Z = (*i)->GetZ();
      if ( (A > 1 && (Z > A || Z <= 0)) || (A==1 && Z > A) || A <= 0 ) return false;
    }
    return true;
}

Here is the call graph for this function:

Here is the caller graph for this function:

CoulombImpulse()

void G4StatMFChannel::CoulombImpulse ( G4int  anA, G4int  anZ, G4double  T  )

private

Definition at line 145 of file G4StatMFChannel.cc.

{
  // First, we have to place the fragments inside of the original nucleus volume
  PlaceFragments(anA);
    
  // Second, we sample initial charged fragments momenta. There are
  // _NumOfChargedFragments charged fragments and they start at the begining
  // of the vector _theFragments (i.e. 0) 
  FragmentsMomenta(_NumOfChargedFragments, 0, T);

  // Third, we have to figure out the asymptotic momenta of charged fragments 
  // For taht we have to solve equations of motion for fragments
  SolveEqOfMotion(anA,anZ,T);

  return;
}

Here is the call graph for this function:

Here is the caller graph for this function:

CreateFragment()

void G4StatMFChannel::CreateFragment	(	G4int	A,
		G4int	Z
	)

Definition at line 87 of file G4StatMFChannel.cc.

{
  if (Z <= 0.5) {
    _theFragments.push_back(new G4StatMFFragment(A,Z));
    _NumOfNeutralFragments++;
  } else {
    _theFragments.push_front(new G4StatMFFragment(A,Z));
    _NumOfChargedFragments++;
  }
    
  return;
}

Here is the caller graph for this function:

FragmentsMomenta()

void G4StatMFChannel::FragmentsMomenta ( G4int  NF, G4int  idx, G4double  T  )

private

Definition at line 220 of file G4StatMFChannel.cc.

{
  G4double KinE = 1.5*T*NF; 
  G4ThreeVector p(0.,0.,0.);
    
  if (NF <= 0) return;
  else if (NF == 1) 
    {
      // We have only one fragment to deal with
      p = IsotropicVector(std::sqrt(2.0*_theFragments[idx]->GetNuclearMass()*KinE));
      _theFragments[idx]->SetMomentum(p);
    } 
  else if (NF == 2) 
    {
      // We have only two fragment to deal with
      G4double M1 = _theFragments[idx]->GetNuclearMass();
      G4double M2 = _theFragments[idx+1]->GetNuclearMass();
      p = IsotropicVector(std::sqrt(2.0*KinE*(M1*M2)/(M1+M2)));     
      _theFragments[idx]->SetMomentum(p);
      _theFragments[idx+1]->SetMomentum(-p);
    } 
  else 
    {
      // We have more than two fragments
      G4double AvailableE;
      G4int i1,i2;
      G4double SummedE;
      G4ThreeVector SummedP(0.,0.,0.);
      do 
    {
      // Fisrt sample momenta of NF-2 fragments 
      // according to Boltzmann distribution
      AvailableE = 0.0;
      SummedE = 0.0;
      SummedP.setX(0.0);SummedP.setY(0.0);SummedP.setZ(0.0);
      for (G4int i = idx; i < idx+NF-2; ++i) 
        {
          G4double E;
          G4double RandE;
          do 
        {
          E = 9.0*G4UniformRand();
          RandE = std::sqrt(0.5/E)*G4Exp(E-0.5)*G4UniformRand();
        } 
          // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
          while (RandE > 1.0);
          E *= T;
          p = IsotropicVector(std::sqrt(2.0*E*_theFragments[i]->GetNuclearMass()));
          _theFragments[i]->SetMomentum(p);
          SummedE += E;
          SummedP += p;
        }   
      // Calculate momenta of last two fragments in such a way
      // that constraints are satisfied
      i1 = idx+NF-2;  // before last fragment index
      i2 = idx+NF-1;  // last fragment index
      p = -SummedP;
      AvailableE = KinE - SummedE;
      // Available Kinetic Energy should be shared between two last fragments
    } 
      // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
      while (AvailableE <= p.mag2()/(2.0*(_theFragments[i1]->GetNuclearMass()+
                      _theFragments[i2]->GetNuclearMass())));
      G4double H = 1.0 + _theFragments[i2]->GetNuclearMass()
    /_theFragments[i1]->GetNuclearMass();
      G4double CTM12 = H*(1.0 - 2.0*_theFragments[i2]->GetNuclearMass()
              *AvailableE/p.mag2());
      G4double CosTheta1;
      G4double Sign;

      if (CTM12 > 1.) {CosTheta1 = 1.;} 
      else {
    do 
      {
        do 
          {
        CosTheta1 = 1.0 - 2.0*G4UniformRand();
          } 
        // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
        while (CosTheta1*CosTheta1 < CTM12);
      }
    // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
    while (CTM12 >= 0.0 && CosTheta1 < 0.0);
      }

      if (CTM12 < 0.0) Sign = 1.0;
      else if (G4UniformRand() <= 0.5) Sign = -1.0;
      else Sign = 1.0;      
        
      G4double P1 = (p.mag()*CosTheta1+Sign*std::sqrt(p.mag2()
                              *(CosTheta1*CosTheta1-CTM12)))/H;
      G4double P2 = std::sqrt(P1*P1+p.mag2() - 2.0*P1*p.mag()*CosTheta1);
      G4double Phi = twopi*G4UniformRand();
      G4double SinTheta1 = std::sqrt(1.0 - CosTheta1*CosTheta1);
      G4double CosPhi1 = std::cos(Phi);
      G4double SinPhi1 = std::sin(Phi);
      G4double CosPhi2 = -CosPhi1;
      G4double SinPhi2 = -SinPhi1;
      G4double CosTheta2 = (p.mag2() + P2*P2 - P1*P1)/(2.0*p.mag()*P2);
      G4double SinTheta2 = 0.0;
      if (CosTheta2 > -1.0 && CosTheta2 < 1.0) {
    SinTheta2 = std::sqrt(1.0 - CosTheta2*CosTheta2);
      }
      G4ThreeVector p1(P1*SinTheta1*CosPhi1,P1*SinTheta1*SinPhi1,P1*CosTheta1);
      G4ThreeVector p2(P2*SinTheta2*CosPhi2,P2*SinTheta2*SinPhi2,P2*CosTheta2);
      G4ThreeVector b(1.0,0.0,0.0);
    
      p1 = RotateMomentum(p,b,p1);
      p2 = RotateMomentum(p,b,p2);
    
      SummedP += p1 + p2;
      SummedE += p1.mag2()/(2.0*_theFragments[i1]->GetNuclearMass()) + 
    p2.mag2()/(2.0*_theFragments[i2]->GetNuclearMass());        
        
      _theFragments[i1]->SetMomentum(p1);
      _theFragments[i2]->SetMomentum(p2);
        
    }
  return;
}

Here is the call graph for this function:

Here is the caller graph for this function:

GetFragments()

G4FragmentVector * G4StatMFChannel::GetFragments	(	G4int	anA,
		G4int	anZ,
		G4double	T
	)

Definition at line 127 of file G4StatMFChannel.cc.

{
  // calculate momenta of charged fragments  
  CoulombImpulse(anA,anZ,T);
    
  // calculate momenta of neutral fragments
  FragmentsMomenta(_NumOfNeutralFragments, _NumOfChargedFragments, T);

  G4FragmentVector * theResult = new G4FragmentVector;
  std::deque<G4StatMFFragment*>::iterator i;
  for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
    theResult->push_back((*i)->GetFragment(T));

  return theResult;
}

Here is the call graph for this function:

Here is the caller graph for this function:

GetFragmentsCoulombEnergy()

G4double G4StatMFChannel::GetFragmentsCoulombEnergy

(

void

)

Definition at line 103 of file G4StatMFChannel.cc.

{
  G4double Coulomb = std::accumulate(_theFragments.begin(),_theFragments.end(),
                     0.0,SumCoulombEnergy());
  //      G4double Coulomb = 0.0;
  //      for (unsigned int i = 0;i < _theFragments.size(); i++)
  //    Coulomb += _theFragments[i]->GetCoulombEnergy();
  return Coulomb;
}

Here is the caller graph for this function:

GetFragmentsEnergy()

G4double G4StatMFChannel::GetFragmentsEnergy

(

G4double

T

)

const

Definition at line 113 of file G4StatMFChannel.cc.

{
  G4double Energy = 0.0;
    
  G4double TranslationalEnergy = 1.5*T*_theFragments.size();

  std::deque<G4StatMFFragment*>::const_iterator i;
  for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
    {
      Energy += (*i)->GetEnergy(T);
    }
  return Energy + TranslationalEnergy;  
}

Here is the caller graph for this function:

GetMultiplicity()

size_t G4StatMFChannel::GetMultiplicity ( void  )

inline

Definition at line 65 of file G4StatMFChannel.hh.

{ return _theFragments.size();}

Here is the call graph for this function:

Here is the caller graph for this function:

IsotropicVector()

G4ThreeVector G4StatMFChannel::IsotropicVector ( G4double  Magnitude = 1.0 )

private

Definition at line 456 of file G4StatMFChannel.cc.

{
    G4double CosTheta = 1.0 - 2.0*G4UniformRand();
    G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
    G4double Phi = twopi*G4UniformRand();
    G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
             Magnitude*std::cos(Phi)*CosTheta,
             Magnitude*std::sin(Phi));
    return Vector;
}

Here is the caller graph for this function:

operator!=()

G4bool G4StatMFChannel::operator!= ( const G4StatMFChannel &  right ) const

private

operator=()

G4StatMFChannel& G4StatMFChannel::operator= ( const G4StatMFChannel &  right )

private

operator==()

G4bool G4StatMFChannel::operator== ( const G4StatMFChannel &  right ) const

private

PlaceFragments()

void G4StatMFChannel::PlaceFragments ( G4int  anA )

private

Definition at line 164 of file G4StatMFChannel.cc.

{
  G4Pow* g4pow = G4Pow::GetInstance();
  const G4double R0 = G4StatMFParameters::Getr0();
  G4double Rsys = 2.0*R0*g4pow->Z13(anA);

  G4bool TooMuchIterations;
  do 
    {
      TooMuchIterations = false;
    
      // Sample the position of the first fragment
      G4double R = (Rsys - R0*g4pow->Z13(_theFragments[0]->GetA()))*
    g4pow->A13(G4UniformRand());
      _theFragments[0]->SetPosition(IsotropicVector(R));
    
    
      // Sample the position of the remaining fragments
      G4bool ThereAreOverlaps = false;
      std::deque<G4StatMFFragment*>::iterator i;
      for (i = _theFragments.begin()+1; i != _theFragments.end(); ++i) 
    {
      G4int counter = 0;
      do 
        {
          R = (Rsys - R0*g4pow->Z13((*i)->GetA()))*g4pow->A13(G4UniformRand());
          (*i)->SetPosition(IsotropicVector(R));
        
          // Check that there are not overlapping fragments
          std::deque<G4StatMFFragment*>::iterator j;
          for (j = _theFragments.begin(); j != i; ++j) 
        {
          G4ThreeVector FragToFragVector = 
            (*i)->GetPosition() - (*j)->GetPosition();
          G4double Rmin = R0*(g4pow->Z13((*i)->GetA()) +
                      g4pow->Z13((*j)->GetA()));
          if ( (ThereAreOverlaps = (FragToFragVector.mag2() < Rmin*Rmin))) 
            { break; }
        }
          counter++;
          // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
        } while (ThereAreOverlaps && counter < 1000);
        
      if (counter >= 1000) 
        {
          TooMuchIterations = true;
          break;
        } 
    }
      // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
    } while (TooMuchIterations);
    return;
}

Here is the call graph for this function:

Here is the caller graph for this function:

RotateMomentum()

G4ThreeVector G4StatMFChannel::RotateMomentum ( G4ThreeVector  Pa, G4ThreeVector  V, G4ThreeVector  P  )

private

Definition at line 436 of file G4StatMFChannel.cc.

{
  G4ThreeVector U = Pa.unit();
  
  G4double Alpha1 = U * V;
  
  G4double Alpha2 = std::sqrt(V.mag2() - Alpha1*Alpha1);

  G4ThreeVector N = (1./Alpha2)*U.cross(V);
  
  G4ThreeVector RotatedMomentum(
                ( (V.x() - Alpha1*U.x())/Alpha2 ) * P.x() + N.x() * P.y() + U.x() * P.z(),
                ( (V.y() - Alpha1*U.y())/Alpha2 ) * P.x() + N.y() * P.y() + U.y() * P.z(),
                ( (V.z() - Alpha1*U.z())/Alpha2 ) * P.x() + N.z() * P.y() + U.z() * P.z()
                );
  return RotatedMomentum;
}

Here is the call graph for this function:

Here is the caller graph for this function:

SolveEqOfMotion()

void G4StatMFChannel::SolveEqOfMotion ( G4int  anA, G4int  anZ, G4double  T  )

private

Definition at line 347 of file G4StatMFChannel.cc.

{
  G4Pow* g4pow = G4Pow::GetInstance();
  G4double CoulombEnergy = 0.6*elm_coupling*anZ*anZ*
    g4pow->A13(1.0+G4StatMFParameters::GetKappaCoulomb())/
    (G4StatMFParameters::Getr0()*g4pow->Z13(anA)) - GetFragmentsCoulombEnergy();
  if (CoulombEnergy <= 0.0) return;
  
  G4int Iterations = 0;
  G4double TimeN = 0.0;
  G4double TimeS = 0.0;
  G4double DeltaTime = 10.0;
  
  G4ThreeVector * Pos = new G4ThreeVector[_NumOfChargedFragments];
  G4ThreeVector * Vel = new G4ThreeVector[_NumOfChargedFragments];
  G4ThreeVector * Accel = new G4ThreeVector[_NumOfChargedFragments];
  
  G4int i;
  for (i = 0; i < _NumOfChargedFragments; i++) 
    {
      Vel[i] = (1.0/(_theFragments[i]->GetNuclearMass()))*
    _theFragments[i]->GetMomentum();
      Pos[i] = _theFragments[i]->GetPosition();
    }

  G4ThreeVector distance(0.,0.,0.);
  G4ThreeVector force(0.,0.,0.);
  G4ThreeVector SavedVel(0.,0.,0.);
  do {
    for (i = 0; i < _NumOfChargedFragments; i++) 
      {
    force.set(0.,0.,0.); 
    for (G4int j = 0; j < _NumOfChargedFragments; j++) 
      {
        if (i != j) 
          {
        distance = Pos[i] - Pos[j];
        force += (elm_coupling*_theFragments[i]->GetZ()
              *_theFragments[j]->GetZ()/
              (distance.mag2()*distance.mag()))*distance;
          }
      }
    Accel[i] = (1./(_theFragments[i]->GetNuclearMass()))*force;
      }

    TimeN = TimeS + DeltaTime;
    
    for ( i = 0; i < _NumOfChargedFragments; i++) 
      {
    SavedVel = Vel[i];
    Vel[i] += Accel[i]*(TimeN-TimeS);
    Pos[i] += (SavedVel+Vel[i])*(TimeN-TimeS)*0.5;
      }
    TimeS = TimeN;
    
    // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
  } while (Iterations++ < 100);
    
  // Summed fragment kinetic energy
  G4double TotalKineticEnergy = 0.0;
  for (i = 0; i < _NumOfChargedFragments; i++) 
    {
      TotalKineticEnergy += _theFragments[i]->GetNuclearMass()*
    0.5*Vel[i].mag2();
    }
  // Scaling of fragment velocities
  G4double KineticEnergy = 1.5*_theFragments.size()*T;
  G4double Eta = ( CoulombEnergy + KineticEnergy ) / TotalKineticEnergy;
  for (i = 0; i < _NumOfChargedFragments; i++) 
    {
      Vel[i] *= Eta;
    }
  
  // Finally calculate fragments momenta
  for (i = 0; i < _NumOfChargedFragments; i++) 
    {
      _theFragments[i]->SetMomentum(_theFragments[i]->GetNuclearMass()*Vel[i]);
    }
  
  // garbage collection
  delete [] Pos;
  delete [] Vel;
  delete [] Accel;
  
  return;
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

_NumOfChargedFragments

G4int G4StatMFChannel::_NumOfChargedFragments

private

Definition at line 106 of file G4StatMFChannel.hh.

_NumOfNeutralFragments

G4int G4StatMFChannel::_NumOfNeutralFragments

private

Definition at line 104 of file G4StatMFChannel.hh.

_theFragments

std::deque<G4StatMFFragment*> G4StatMFChannel::_theFragments

private

Definition at line 102 of file G4StatMFChannel.hh.

---

The documentation for this class was generated from the following files:

• G4StatMFChannel.hh
• G4StatMFChannel.cc