#include <G4QMDGroundStateNucleus.hh>

Inheritance diagram for G4QMDGroundStateNucleus:

Public Member Functions
void	CalEnergyAndAngularMomentumInCM ()

void	Clear ()

void	DeleteParticipant (G4int i)

G4QMDParticipant *	EraseParticipant (G4int i)

	G4QMDGroundStateNucleus ()

	G4QMDGroundStateNucleus (G4int z, G4int a)

G4LorentzVector	Get4Momentum ()

G4int	GetAngularMomentum ()

G4int	GetAtomicNumber ()

G4double	GetExcitationEnergy ()

G4int	GetMassNumber ()

G4int	GetNOCollision ()

G4double	GetNuclearMass ()

G4QMDParticipant *	GetParticipant (G4int i)

G4int	GetTotalNumberOfParticipant ()

void	IncrementCollisionCounter ()

void	InsertParticipant (G4QMDParticipant *particle, G4int j)

void	SetParticipant (G4QMDParticipant *particle)

void	SetSystem (G4QMDSystem *, G4ThreeVector, G4ThreeVector)

void	SetTotalPotential (G4double x)

void	ShowParticipants ()

void	SubtractSystem (G4QMDSystem *)

	~G4QMDGroundStateNucleus ()

Protected Attributes
std::vector< G4QMDParticipant * >	participants

Private Member Functions
void	killCMMotionAndAngularM ()

void	packNucleons ()

G4bool	samplingMomentum (G4int)

G4bool	samplingPosition (G4int)

Private Attributes
G4double	c0p

G4double	c3p

G4double	cdp

G4double	clp

G4double	cpc

G4double	cph

G4double	cpw

G4double	csp

std::vector< G4double >	d_pot

G4double	ddif

G4double	ddif2

G4double	dsam

G4double	dsam2

G4double	ebini

G4double	edepth

G4double	epse

G4double	epsx

std::vector< G4double >	es

G4double	excitationEnergy

G4double	gamm

G4double	hbc

G4int	jj

G4int	maxTrial

G4QMDMeanField *	meanfield

G4int	numberOfCollision

std::vector< G4ThreeVector >	pcm

std::vector< G4double >	phase_g

G4double	potentialEnergy

G4double	r00

G4double	r01

G4double	rada

G4double	radb

G4double	radm

std::vector< G4ThreeVector >	rcm

std::vector< G4double >	rho_l

G4double	rmax

G4double	rt00

G4double	saa

Detailed Description

Definition at line 45 of file G4QMDGroundStateNucleus.hh.

Constructor & Destructor Documentation

◆ G4QMDGroundStateNucleus() [1/2]

G4QMDGroundStateNucleus::G4QMDGroundStateNucleus ( )

◆ G4QMDGroundStateNucleus() [2/2]

G4QMDGroundStateNucleus::G4QMDGroundStateNucleus	(	G4int	z,
		G4int	a
	)

Definition at line 38 of file G4QMDGroundStateNucleus.cc.

: maxTrial ( 1000 )
, r00 ( 1.124 )  // radius parameter for Woods-Saxon [fm] 
, r01 ( 0.5 )    // radius parameter for Woods-Saxon 
, saa ( 0.2 )    // diffuse parameter for initial Woods-Saxon shape
, rada ( 0.9 )   // cutoff parameter
, radb ( 0.3 )   // cutoff parameter
, dsam ( 1.5 )   // minimum distance for same particle [fm]
, ddif ( 1.0 )   // minimum distance for different particle
, edepth ( 0.0 )
, epse ( 0.000001 )  // torelance for energy in [GeV]
, meanfield ( NULL ) 
{
 
   //std::cout << " G4QMDGroundStateNucleus( G4int z , G4int a ) Begin " << z << " " << a << std::endl;
 
   dsam2 = dsam*dsam;
   ddif2 = ddif*ddif;
 
   G4QMDParameters* parameters = G4QMDParameters::GetInstance();
 
   hbc = parameters->Get_hbc();
   gamm = parameters->Get_gamm();
   cpw = parameters->Get_cpw();
   cph = parameters->Get_cph();
   epsx = parameters->Get_epsx();
   cpc = parameters->Get_cpc();
 
   cdp = parameters->Get_cdp();
   c0p = parameters->Get_c0p();
   c3p = parameters->Get_c3p();
   csp = parameters->Get_csp();
   clp = parameters->Get_clp();
 
   // Following 10 lines should be here, right before the line 90.
   // Otherwise, mass number cannot be conserved if the projectile or
   // the target are nucleons.
   //Nucleon primary or target case;
   if ( z == 1 && a == 1 ) {  // Hydrogen  Case or proton primary 
      SetParticipant( new G4QMDParticipant( G4Proton::Proton() , G4ThreeVector( 0.0 ) , G4ThreeVector( 0.0 ) ) );
      ebini = 0.0; 
      return;
   }
   else if ( z == 0 && a == 1 ) { // Neutron primary 
      SetParticipant( new G4QMDParticipant( G4Neutron::Neutron() , G4ThreeVector( 0.0 ) , G4ThreeVector( 0.0 ) ) );
      ebini = 0.0; 
      return;
   }
 
 
   //edepth = 0.0; 
 
   for ( int i = 0 ; i < a ; i++ )
   {
 
      G4ParticleDefinition* pd; 
 
      if ( i < z )
      { 
         pd = G4Proton::Proton();
      }
      else
      {
         pd = G4Neutron::Neutron();
      }
         
      G4ThreeVector p( 0.0 );
      G4ThreeVector r( 0.0 );
      G4QMDParticipant* aParticipant = new G4QMDParticipant( pd , p , r );
      SetParticipant( aParticipant );
 
   }
 
   G4double radious = r00 * G4Pow::GetInstance()->A13( double ( GetMassNumber() ) ); 
 
   rt00 = radious - r01; 
   radm = radious - rada * ( gamm - 1.0 ) + radb;
   rmax = 1.0 / ( 1.0 + G4Exp ( -rt00/saa ) );
 
   //maxTrial = 1000;
   
   
   meanfield = new G4QMDMeanField();
   meanfield->SetSystem( this );
 
   //std::cout << "G4QMDGroundStateNucleus( G4int z , G4int a ) packNucleons Begin ( z , a ) ( " << z << ", " << a << " )" << std::endl;
   packNucleons();
   //std::cout << "G4QMDGroundStateNucleus( G4int z , G4int a ) packNucleons End" << std::endl;
 
   delete meanfield;
 
}

◆ ~G4QMDGroundStateNucleus()

G4QMDGroundStateNucleus::~G4QMDGroundStateNucleus ( )

inline

Definition at line 50 of file G4QMDGroundStateNucleus.hh.

      {
         rho_l.clear();  
         d_pot.clear();
      };

References d_pot, and rho_l.

Member Function Documentation

◆ CalEnergyAndAngularMomentumInCM()

void G4QMDNucleus::CalEnergyAndAngularMomentumInCM ( )

inherited

Definition at line 138 of file G4QMDNucleus.cc.

{
 
   //G4cout << "CalEnergyAndAngularMomentumInCM " << this->GetAtomicNumber() << " " << GetMassNumber() << G4endl;
 
   G4double gamma = Get4Momentum().gamma();
   G4ThreeVector beta = Get4Momentum().v()/ Get4Momentum().e();
 
   G4ThreeVector pcm0( 0.0 ) ;
 
   G4int n = GetTotalNumberOfParticipant();
   pcm.resize( n );
 
   for ( G4int i= 0; i < n ; i++ ) 
   {
      G4ThreeVector p_i = GetParticipant( i )->GetMomentum();
 
      G4double trans = gamma / ( gamma + 1.0 ) * p_i * beta; 
      pcm[i] = p_i - trans*beta;
 
      pcm0 += pcm[i];
   }
 
   pcm0 = pcm0 / double ( n );
 
   //G4cout << "pcm0 " << pcm0 << G4endl;
 
   for ( G4int i= 0; i < n ; i++ ) 
   {
      pcm[i] += -pcm0;
      //G4cout << "pcm " << i << " " << pcm[i] << G4endl;
   }
 
 
   G4double tmass = 0;
   G4ThreeVector rcm0( 0.0 ) ;
   rcm.resize( n );
   es.resize( n );
 
   for ( G4int i= 0; i < n ; i++ ) 
   {
      G4ThreeVector ri = GetParticipant( i )->GetPosition();
      G4double trans = gamma / ( gamma + 1.0 ) * ri * beta; 
 
      es[i] = std::sqrt ( G4Pow::GetInstance()->powN ( GetParticipant( i )->GetMass() , 2 ) + pcm[i]*pcm[i] );
 
      rcm[i] = ri + trans*beta;
 
      rcm0 += rcm[i]*es[i];
 
      tmass += es[i];
   }
 
   rcm0 = rcm0/tmass;
 
   for ( G4int i= 0; i < n ; i++ ) 
   {
      rcm[i] += -rcm0;
      //G4cout << "rcm " << i << " " << rcm[i] << G4endl;
   }
 
// Angular momentum
 
   G4ThreeVector rl ( 0.0 ); 
   for ( G4int i= 0; i < n ; i++ ) 
   {
      rl += rcm[i].cross ( pcm[i] );
   }
 
// DHW: move hbc outside of sqrt to get correct units
//  jj = int ( std::sqrt ( rl*rl / hbc ) + 0.5 );
 
   jj = int (std::sqrt(rl*rl)/hbc + 0.5);
 
// kinetic energy per nucleon in CM
 
   G4double totalMass = 0.0;
   for ( G4int i= 0; i < n ; i++ ) 
   {
      // following two lines are equivalent
      //totalMass += GetParticipant( i )->GetDefinition()->GetPDGMass()/GeV;
      totalMass += GetParticipant( i )->GetMass();
   }
 
   //G4double kineticEnergyPerNucleon = ( std::accumulate ( es.begin() , es.end() , 0.0 ) - totalMass )/n;
 
// Total (not per nucleion ) Binding Energy 
   G4double bindingEnergy =  ( std::accumulate ( es.begin() , es.end() , 0.0 ) -totalMass ) + potentialEnergy;
 
   //G4cout << "KineticEnergyPerNucleon in GeV " << kineticEnergyPerNucleon << G4endl;
   //G4cout << "KineticEnergySum in GeV " << std::accumulate ( es.begin() , es.end() , 0.0 ) - totalMass << G4endl;
   //G4cout << "PotentialEnergy in GeV " << potentialEnergy << G4endl;
   //G4cout << "BindingEnergy in GeV " << bindingEnergy << G4endl;
   //G4cout << "G4BindingEnergy in GeV " << G4NucleiProperties::GetBindingEnergy( GetAtomicNumber() , GetMassNumber() )/GeV << G4endl;
 
   excitationEnergy = bindingEnergy + G4NucleiProperties::GetBindingEnergy( GetMassNumber() , GetAtomicNumber() )/GeV;
   //G4cout << "excitationEnergy in GeV " << excitationEnergy << G4endl;
   if ( excitationEnergy < 0 ) excitationEnergy = 0.0; 
 
}

References anonymous_namespace{G4PionRadiativeDecayChannel.cc}::beta, G4InuclSpecialFunctions::bindingEnergy(), CLHEP::HepLorentzVector::e(), G4QMDNucleus::es, G4QMDNucleus::excitationEnergy, CLHEP::HepLorentzVector::gamma(), G4QMDNucleus::Get4Momentum(), G4QMDNucleus::GetAtomicNumber(), G4NucleiProperties::GetBindingEnergy(), G4Pow::GetInstance(), G4QMDParticipant::GetMass(), G4QMDNucleus::GetMassNumber(), G4QMDParticipant::GetMomentum(), G4QMDSystem::GetParticipant(), G4QMDParticipant::GetPosition(), G4QMDSystem::GetTotalNumberOfParticipant(), GeV, G4QMDNucleus::hbc, G4QMDNucleus::jj, CLHEP::detail::n, G4QMDNucleus::pcm, G4QMDNucleus::potentialEnergy, G4QMDNucleus::rcm, and CLHEP::HepLorentzVector::v().

Referenced by G4QMDReaction::ApplyYourself(), and G4QMDMeanField::SetNucleus().

◆ Clear()

void G4QMDSystem::Clear ( )

inherited

Definition at line 68 of file G4QMDSystem.cc.

{
   for ( G4int i = 0 ; i < this->GetTotalNumberOfParticipant() ; i++ )
   {
      delete participants[i];
   }
   participants.clear();
}

References G4QMDSystem::GetTotalNumberOfParticipant(), and G4QMDSystem::participants.

Referenced by G4QMDReaction::ApplyYourself(), and G4QMDSystem::~G4QMDSystem().

◆ DeleteParticipant()

void G4QMDSystem::DeleteParticipant ( G4int i )

inlineinherited

Definition at line 57 of file G4QMDSystem.hh.

57{ delete participants[ i ] ; participants.erase( std::find ( participants.begin() , participants.end() , participants[ i ] ) ); };

References G4QMDSystem::participants.

Referenced by G4QMDCollision::CalKinematicsOfBinaryCollisions().

◆ EraseParticipant()

G4QMDParticipant * G4QMDSystem::EraseParticipant ( G4int i )

inlineinherited

Definition at line 56 of file G4QMDSystem.hh.

56{ G4QMDParticipant* particle = participants[ i ]; participants.erase( std::find ( participants.begin() , participants.end() , participants[ i ] ) ) ; return particle; };

References G4QMDSystem::participants.

Referenced by G4QMDCollision::CalFinalStateOfTheBinaryCollision().

◆ Get4Momentum()

G4LorentzVector G4QMDNucleus::Get4Momentum ( )

inherited

Definition at line 56 of file G4QMDNucleus.cc.

{
   G4LorentzVector p( 0 );
   std::vector< G4QMDParticipant* >::iterator it; 
   for ( it = participants.begin() ; it != participants.end() ; it++ ) 
      p += (*it)->Get4Momentum();   
 
   return p;
}

References G4QMDSystem::participants.

Referenced by G4QMDNucleus::CalEnergyAndAngularMomentumInCM().

◆ GetAngularMomentum()

G4int G4QMDNucleus::GetAngularMomentum ( )

inlineinherited

Definition at line 65 of file G4QMDNucleus.hh.

65{ return jj; };

References G4QMDNucleus::jj.

◆ GetAtomicNumber()

G4int G4QMDNucleus::GetAtomicNumber ( )

inherited

Definition at line 89 of file G4QMDNucleus.cc.

{
   G4int Z = 0; 
   std::vector< G4QMDParticipant* >::iterator it; 
   for ( it = participants.begin() ; it != participants.end() ; it++ ) 
   {
      if ( (*it)->GetDefinition() == G4Proton::Proton() ) 
         Z++; 
   }
   return Z;
}

References G4QMDSystem::participants, G4Proton::Proton(), and Z.

Referenced by G4QMDNucleus::CalEnergyAndAngularMomentumInCM(), G4QMDNucleus::GetNuclearMass(), and packNucleons().

◆ GetExcitationEnergy()

G4double G4QMDNucleus::GetExcitationEnergy ( void )

inlineinherited

Definition at line 63 of file G4QMDNucleus.hh.

63{ return excitationEnergy; };

References G4QMDNucleus::excitationEnergy.

◆ GetMassNumber()

G4int G4QMDNucleus::GetMassNumber ( )

inherited

Definition at line 68 of file G4QMDNucleus.cc.

{
 
   G4int A = 0; 
   std::vector< G4QMDParticipant* >::iterator it; 
   for ( it = participants.begin() ; it != participants.end() ; it++ ) 
   {
      if ( (*it)->GetDefinition() == G4Proton::Proton() 
        || (*it)->GetDefinition() == G4Neutron::Neutron() ) 
         A++; 
   }
 
   if ( A == 0 ) {
      throw G4HadronicException(__FILE__, __LINE__, "G4QMDNucleus has the mass number of 0!");
   }
 
   return A;
}

References A, G4Neutron::Neutron(), G4QMDSystem::participants, and G4Proton::Proton().

Referenced by G4QMDNucleus::CalEnergyAndAngularMomentumInCM(), G4QMDGroundStateNucleus(), G4QMDNucleus::GetNuclearMass(), killCMMotionAndAngularM(), packNucleons(), and samplingMomentum().

◆ GetNOCollision()

G4int G4QMDSystem::GetNOCollision ( )

inlineinherited

Definition at line 65 of file G4QMDSystem.hh.

65{ return numberOfCollision; };

G4QMDSystem::numberOfCollision

G4int numberOfCollision

Definition: G4QMDSystem.hh:75

References G4QMDSystem::numberOfCollision.

Referenced by G4QMDReaction::ApplyYourself().

◆ GetNuclearMass()

G4double G4QMDNucleus::GetNuclearMass ( void )

inherited

Definition at line 103 of file G4QMDNucleus.cc.

{
 
   G4double mass = G4NucleiProperties::GetNuclearMass( GetMassNumber() , GetAtomicNumber() );
   
   if ( mass == 0.0 )
   {
 
      G4int Z = GetAtomicNumber();
      G4int A = GetMassNumber();
      G4int N = A - Z;
 
// Weizsacker-Bethe 
 
      G4double Av = 16*MeV; 
      G4double As = 17*MeV; 
      G4double Ac = 0.7*MeV; 
      G4double Asym = 23*MeV; 
 
      G4double BE = Av * A 
                  - As * G4Pow::GetInstance()->A23 ( G4double ( A ) ) 
                  - Ac * Z*Z/G4Pow::GetInstance()->A13 ( G4double ( A ) )
                  - Asym * ( N - Z )* ( N - Z ) / A; 
 
      mass = Z * G4Proton::Proton()->GetPDGMass() 
           + N * G4Neutron::Neutron()->GetPDGMass()
           - BE;
 
   }
 
   return mass;
}

References A, G4Pow::A13(), G4Pow::A23(), G4QMDNucleus::GetAtomicNumber(), G4Pow::GetInstance(), G4QMDNucleus::GetMassNumber(), G4NucleiProperties::GetNuclearMass(), G4ParticleDefinition::GetPDGMass(), MeV, G4Neutron::Neutron(), G4Proton::Proton(), and Z.

◆ GetParticipant()

G4QMDParticipant * G4QMDSystem::GetParticipant ( G4int i )

inlineinherited

Definition at line 62 of file G4QMDSystem.hh.

62{ return participants[i]; };

References G4QMDSystem::participants.

◆ GetTotalNumberOfParticipant()

G4int G4QMDSystem::GetTotalNumberOfParticipant ( )

inlineinherited

Definition at line 60 of file G4QMDSystem.hh.

60{ return participants.size(); };

References G4QMDSystem::participants.

Referenced by G4QMDReaction::ApplyYourself(), G4QMDMeanField::Cal2BodyQuantities(), G4QMDNucleus::CalEnergyAndAngularMomentumInCM(), G4QMDMeanField::CalGraduate(), G4QMDCollision::CalKinematicsOfBinaryCollisions(), G4QMDMeanField::calPauliBlockingFactor(), G4QMDSystem::Clear(), G4QMDMeanField::DoClusterJudgment(), G4QMDMeanField::DoPropagation(), G4QMDMeanField::GetPotential(), G4QMDMeanField::GetTotalPotential(), G4QMDMeanField::SetSystem(), and G4QMDSystem::SubtractSystem().

◆ IncrementCollisionCounter()

void G4QMDSystem::IncrementCollisionCounter ( )

inlineinherited

Definition at line 64 of file G4QMDSystem.hh.

64{ numberOfCollision++; };

References G4QMDSystem::numberOfCollision.

Referenced by G4QMDCollision::CalKinematicsOfBinaryCollisions().

◆ InsertParticipant()

void G4QMDSystem::InsertParticipant	(	G4QMDParticipant *	particle,
		G4int	j
	)

inherited

Definition at line 110 of file G4QMDSystem.cc.

{
 
   if ( (size_t) n > participants.size()+1 )
      G4cout << "G4QMDSystem::InsertParticipant size error" << G4endl;
 
   std::vector< G4QMDParticipant* >::iterator it; 
   it = participants.begin();
 
   for ( G4int i = 0; i < n ; i++ )
      it++;
 
   participants.insert( it, particle );
}

References G4cout, G4endl, CLHEP::detail::n, and G4QMDSystem::participants.

Referenced by G4QMDCollision::CalFinalStateOfTheBinaryCollision().

◆ killCMMotionAndAngularM()

void G4QMDGroundStateNucleus::killCMMotionAndAngularM ( )

private

Definition at line 730 of file G4QMDGroundStateNucleus.cc.

{
 
//   CalEnergyAndAngularMomentumInCM();
 
   //std::vector< G4ThreeVector > p ( GetMassNumber() , 0.0 );
   //std::vector< G4ThreeVector > r ( GetMassNumber() , 0.0 );
 
// Move to cm system
 
   G4ThreeVector pcm_tmp ( 0.0 );
   G4ThreeVector rcm_tmp ( 0.0 );
   G4double sumMass = 0.0;
 
   for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
   {
      pcm_tmp += participants[i]->GetMomentum();
      rcm_tmp += participants[i]->GetPosition() * participants[i]->GetMass();
      sumMass += participants[i]->GetMass();
   }
 
   pcm_tmp = pcm_tmp/GetMassNumber();
   rcm_tmp = rcm_tmp/sumMass;
 
   for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
   {
      participants[i]->SetMomentum( participants[i]->GetMomentum() - pcm_tmp );
      participants[i]->SetPosition( participants[i]->GetPosition() - rcm_tmp );
   }
 
// kill the angular momentum
 
   G4ThreeVector ll ( 0.0 );
   for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
   {
      ll += participants[i]->GetPosition().cross ( participants[i]->GetMomentum() );
   }
 
   G4double rr[3][3];
   G4double ss[3][3];
   for ( G4int i = 0 ; i < 3 ; i++ )
   {
      for ( G4int j = 0 ; j < 3 ; j++ )
      {
         rr [i][j] = 0.0;
 
         if ( i == j ) 
         {
            ss [i][j] = 1.0;
         }
         else
         {
            ss [i][j] = 0.0;
         }
      } 
   }
 
   for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
   {
      G4ThreeVector r = participants[i]->GetPosition();
      rr[0][0] += r.y() * r.y() + r.z() * r.z(); 
      rr[1][0] += - r.y() * r.x();
      rr[2][0] += - r.z() * r.x();
      rr[0][1] += - r.x() * r.y();
      rr[1][1] += r.z() * r.z() + r.x() * r.x(); 
      rr[2][1] += - r.z() * r.y();
      rr[2][0] += - r.x() * r.z();
      rr[2][1] += - r.y() * r.z();
      rr[2][2] += r.x() * r.x() + r.y() * r.y(); 
   }
 
   for ( G4int i = 0 ; i < 3 ; i++ )
   {
      G4double x = rr [i][i];
      for ( G4int j = 0 ; j < 3 ; j++ )
      {
         rr[i][j] = rr[i][j] / x;
         ss[i][j] = ss[i][j] / x;
      }
      for ( G4int j = 0 ; j < 3 ; j++ )
      {
         if ( j != i ) 
         {
            G4double y = rr [j][i];
            for ( G4int k = 0 ; k < 3 ; k++ )
            {
               rr[j][k] += -y * rr[i][k];
               ss[j][k] += -y * ss[i][k];
            }
         }
      }
   }
 
   G4double opl[3];
   G4double rll[3];
 
   rll[0] = ll.x();
   rll[1] = ll.y();
   rll[2] = ll.z();
   
   for ( G4int i = 0 ; i < 3 ; i++ )
   {
      opl[i] = 0.0;
 
      for ( G4int j = 0 ; j < 3 ; j++ )
      {
     opl[i] += ss[i][j]*rll[j];
      }
   }
 
   for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
   {
      G4ThreeVector p_i = participants[i]->GetMomentum() ;
      G4ThreeVector ri = participants[i]->GetPosition() ;
      G4ThreeVector opl_v ( opl[0] , opl[1] , opl[2] );  
 
      p_i += ri.cross(opl_v);
 
      participants[i]->SetMomentum( p_i );
   }
 
}

References CLHEP::Hep3Vector::cross(), G4QMDNucleus::GetMassNumber(), G4QMDSystem::participants, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

Referenced by packNucleons().

◆ packNucleons()

void G4QMDGroundStateNucleus::packNucleons ( )

private

Definition at line 133 of file G4QMDGroundStateNucleus.cc.

{
 
   //std::cout << "G4QMDGroundStateNucleus::packNucleons" << std::endl;
 
   ebini = - G4NucleiProperties::GetBindingEnergy( GetMassNumber() , GetAtomicNumber() ) / GetMassNumber();
 
   G4double ebin00 = ebini * 0.001;
 
   G4double ebin0 = 0.0;  
   G4double ebin1 = 0.0;  
 
   if ( GetMassNumber() != 4  )
   {
      ebin0 = ( ebini - 0.5 ) * 0.001;
      ebin1 = ( ebini + 0.5 ) * 0.001;
   }
   else
   {
      ebin0 = ( ebini - 1.5 ) * 0.001;
      ebin1 = ( ebini + 1.5 ) * 0.001;
   }
 
   G4int n0Try = 0; 
   G4bool isThisOK = false;
   while ( n0Try < maxTrial ) // Loop checking, 11.03.2015, T. Koi
   {
      n0Try++;
      //std::cout << "TKDB packNucleons n0Try " << n0Try << std::endl;
 
//    Sampling Position
 
      //std::cout << "TKDB Sampling Position " << std::endl;
 
      G4bool areThesePsOK = false;
      G4int npTry = 0;
      while ( npTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi
      {
      //std::cout << "TKDB Sampling Position npTry " << npTry << std::endl;
         npTry++; 
         G4int i = 0; 
         if ( samplingPosition( i ) ) 
         {
            //std::cout << "packNucleons samplingPosition 0 succeed " << std::endl;
            for ( i = 1 ; i < GetMassNumber() ; i++ )
            {
               //std::cout << "packNucleons samplingPosition " << i  << " trying " << std::endl;
               if ( !( samplingPosition( i ) ) ) 
               {
                  //std::cout << "packNucleons samplingPosition " << i << " failed" << std::endl;
                  break;
               }
            }
            if ( i == GetMassNumber() ) 
            {
               //std::cout << "packNucleons samplingPosition all scucceed " << std::endl;
               areThesePsOK = true;
               break; 
            }
         }
      }
      if ( areThesePsOK == false ) continue;
 
      //std::cout << "TKDB Sampling Position End" << std::endl;
 
//    Calculate Two-body quantities
 
      meanfield->Cal2BodyQuantities(); 
      std::vector< G4double > rho_a ( GetMassNumber() , 0.0 );
      std::vector< G4double > rho_s ( GetMassNumber() , 0.0 );
      std::vector< G4double > rho_c ( GetMassNumber() , 0.0 );
 
      rho_l.resize ( GetMassNumber() , 0.0 );
      d_pot.resize ( GetMassNumber() , 0.0 );
 
      for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
      {
         for ( G4int j = 0 ; j < GetMassNumber() ; j++ )
         {
 
            rho_a[ i ] += meanfield->GetRHA( j , i ); 
            G4int k = 0; 
            if ( participants[j]->GetDefinition() != participants[i]->GetDefinition() )
            {
               k = 1;
            } 
 
            rho_s[ i ] += meanfield->GetRHA( j , i )*( 1.0 - 2.0 * k ); // OK?  
 
            rho_c[ i ] += meanfield->GetRHE( j , i ); 
         }
 
      }
 
      for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
      {
         rho_l[i] = cdp * ( rho_a[i] + 1.0 );     
         d_pot[i] = c0p * rho_a[i]
                  + c3p * G4Pow::GetInstance()->powA ( rho_a[i] , gamm )
                  + csp * rho_s[i]
                  + clp * rho_c[i];
 
         //std::cout << "d_pot[i] " << i << " " << d_pot[i] << std::endl; 
      }
 
 
// Sampling Momentum
 
      //std::cout << "TKDB Sampling Momentum " << std::endl;
 
      phase_g.clear();
      phase_g.resize( GetMassNumber() , 0.0 );
      
      //std::cout << "TKDB Sampling Momentum 1st " << std::endl;
      G4bool isThis1stMOK = false;
      G4int nmTry = 0; 
      while ( nmTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi
      {
         nmTry++;
         G4int i = 0; 
         if ( samplingMomentum( i ) ) 
         {
           isThis1stMOK = true;
           break; 
         }
      }
      if ( isThis1stMOK == false ) continue;
 
      //std::cout << "TKDB Sampling Momentum 2nd so on" << std::endl;
 
      G4bool areTheseMsOK = true;
      nmTry = 0; 
      while ( nmTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi
      {
         nmTry++;
            G4int i = 0; 
            for ( i = 1 ; i < GetMassNumber() ; i++ )
            {
               //std::cout << "TKDB packNucleons samplingMomentum try " << i << std::endl;
               if ( !( samplingMomentum( i ) ) ) 
               {
                  //std::cout << "TKDB packNucleons samplingMomentum " << i << " failed" << std::endl;
                  areTheseMsOK = false;
                  break;
               }
            }
            if ( i == GetMassNumber() ) 
            {
               areTheseMsOK = true;
            }
 
            break;
      }
      if ( areTheseMsOK == false ) continue;
     
// Kill Angluar Momentum
 
      //std::cout << "TKDB Sampling Kill Angluar Momentum " << std::endl;
 
      killCMMotionAndAngularM();    
 
 
// Check Binding Energy
 
      //std::cout << "packNucleons Check Binding Energy Begin " << std::endl;
 
      G4double ekinal = 0.0;
      for ( int i = 0 ; i < GetMassNumber() ; i++ )
      {
         ekinal += participants[i]->GetKineticEnergy();
      }
 
      meanfield->Cal2BodyQuantities();
 
      G4double totalPotentialE = meanfield->GetTotalPotential(); 
      G4double ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() );
 
      //std::cout << "packNucleons totalPotentialE " << totalPotentialE << std::endl;
      //std::cout << "packNucleons ebinal " << ebinal << std::endl;
      //std::cout << "packNucleons ekinal " << ekinal << std::endl;
 
      if ( ebinal < ebin0 || ebinal > ebin1 ) 
      {
         //std::cout << "packNucleons ebin0 " << ebin0 << std::endl;
         //std::cout << "packNucleons ebin1 " << ebin1 << std::endl;
         //std::cout << "packNucleons ebinal " << ebinal << std::endl;
      //std::cout << "packNucleons Check Binding Energy Failed " << std::endl;
         continue;
      }
 
      //std::cout << "packNucleons Check Binding Energy End = OK " << std::endl;
 
 
// Energy Adujstment
 
      G4double dtc = 1.0;
      G4double frg = -0.1;
      G4double rdf0 = 0.5;
      
      G4double edif0 = ebinal - ebin00;
 
      G4double cfrc = 0.0;
      if ( 0 < edif0 ) 
      {
         cfrc = frg;
      }
      else
      {
         cfrc = -frg;
      }
 
      G4int ifrc = 1;
 
      G4int neaTry = 0;
 
      G4bool isThisEAOK = false;
      while ( neaTry < maxTrial )  // Loop checking, 11.03.2015, T. Koi
      {
         neaTry++;
 
         G4double  edif = ebinal - ebin00; 
 
         //std::cout << "TKDB edif " << edif << std::endl; 
         if ( std::abs ( edif ) < epse )
         {
   
            isThisEAOK = true;
            //std::cout << "isThisEAOK " << isThisEAOK << std::endl; 
            break;
         }
 
         G4int jfrc = 0;
         if ( edif < 0.0 ) 
         {
            jfrc = 1;
         }
         else
         {
            jfrc = -1;
         }
 
         if ( jfrc != ifrc ) 
         {
            cfrc = -rdf0 * cfrc;
            dtc = rdf0 * dtc;
         }
 
         if ( jfrc == ifrc && std::abs( edif0 ) < std::abs( edif ) )
         {
            cfrc = -rdf0 * cfrc;
            dtc = rdf0 * dtc;
         }
 
         ifrc = jfrc;
         edif0 = edif;
 
         meanfield->CalGraduate();
 
         for ( int i = 0 ; i < GetMassNumber() ; i++ )
         {
            G4ThreeVector ri = participants[i]->GetPosition(); 
            G4ThreeVector p_i = participants[i]->GetMomentum(); 
 
            ri += dtc * ( meanfield->GetFFr(i) - cfrc * ( meanfield->GetFFp(i) ) );
            p_i += dtc * ( meanfield->GetFFp(i) + cfrc * ( meanfield->GetFFr(i) ) );
 
            participants[i]->SetPosition( ri ); 
            participants[i]->SetMomentum( p_i ); 
         }
 
         ekinal = 0.0;     
 
         for ( int i = 0 ; i < GetMassNumber() ; i++ )
         {
            ekinal += participants[i]->GetKineticEnergy(); 
         }
 
         meanfield->Cal2BodyQuantities(); 
         totalPotentialE = meanfield->GetTotalPotential(); 
 
         ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() );
 
      }
      //std::cout << "isThisEAOK " << isThisEAOK << std::endl; 
      if ( isThisEAOK == false ) continue;
   
      isThisOK = true;
      //std::cout << "isThisOK " << isThisOK << std::endl; 
      break; 
 
   }
 
   if ( isThisOK == false )
   {
      G4cout << "GroundStateNucleus state cannot be created. Try again with another parameters." << G4endl;
   } 
 
   //std::cout << "packNucleons End" << std::endl;
   return;
}

References c0p, c3p, G4QMDMeanField::Cal2BodyQuantities(), G4QMDMeanField::CalGraduate(), cdp, clp, csp, d_pot, ebini, epse, G4cout, G4endl, gamm, G4QMDNucleus::GetAtomicNumber(), G4NucleiProperties::GetBindingEnergy(), G4QMDMeanField::GetFFp(), G4QMDMeanField::GetFFr(), G4Pow::GetInstance(), G4QMDNucleus::GetMassNumber(), G4QMDMeanField::GetRHA(), G4QMDMeanField::GetRHE(), G4QMDMeanField::GetTotalPotential(), killCMMotionAndAngularM(), maxTrial, meanfield, G4QMDSystem::participants, phase_g, G4Pow::powA(), rho_l, samplingMomentum(), and samplingPosition().

Referenced by G4QMDGroundStateNucleus().

◆ samplingMomentum()

G4bool G4QMDGroundStateNucleus::samplingMomentum ( G4int i )

private

Definition at line 536 of file G4QMDGroundStateNucleus.cc.

{
 
   //std::cout << "TKDB samplingMomentum for " << i << std::endl;
   
   G4bool result = false;
 
   G4double pfm = hbc * G4Pow::GetInstance()->A13 ( ( 3.0 / 2.0 * pi*pi * rho_l[i] ) );
 
   if ( 10 < GetMassNumber() &&  -5.5 < ebini ) 
   {
      pfm = pfm * ( 1.0 + 0.2 * std::sqrt( std::abs( 8.0 + ebini ) / 8.0 ) );
   }
 
   //std::cout << "TKDB samplingMomentum pfm " << pfm << std::endl;
 
   std::vector< G4double > phase; 
   phase.resize( i+1 ); // i start from 0
 
   G4int ntry = 0;
// 710 
   while ( ntry < maxTrial )  // Loop checking, 11.03.2015, T. Koi
   {
 
      //std::cout << " TKDB ntry " << ntry << std::endl;
      ntry++;
 
      G4double ke = DBL_MAX;
 
      G4int tkdb_i =0;
// 700
      G4int icounter = 0;
      G4int icounter_max = 1024;
      while ( ke + d_pot [i] > edepth ) // Loop checking, 11.03.2015, T. Koi
      {
         icounter++;
         if ( icounter > icounter_max ) {
        G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
            break;
         }
      
         G4double psqr = 10.0;
         G4double px = 0.0;
         G4double py = 0.0;
         G4double pz = 0.0;
 
         G4int jcounter = 0;
         G4int jcounter_max = 1024;
         while ( psqr > 1.0 ) // Loop checking, 11.03.2015, T. Koi
         {
            jcounter++;
            if ( jcounter > jcounter_max ) {
           G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
               break;
            }
            px = 1.0 - 2.0*G4UniformRand();
            py = 1.0 - 2.0*G4UniformRand();
            pz = 1.0 - 2.0*G4UniformRand();
 
            psqr = px*px + py*py + pz*pz;
         }
 
         G4ThreeVector p ( px , py , pz ); 
         p = pfm * p;
         participants[i]->SetMomentum( p );
         G4LorentzVector p4 = participants[i]->Get4Momentum();
         //ke = p4.e() - p4.restMass();
         ke = participants[i]->GetKineticEnergy();
   
 
         tkdb_i++;  
         if ( tkdb_i > maxTrial ) return result; // return false
 
      }
 
      //std::cout << "TKDB ke d_pot[i] " << ke << " " << d_pot[i] << std::endl;
 
 
      if ( i == 0 )   
      {
         result = true; 
         return result;
      }
 
      G4bool isThisOK = true;
 
      // Check Pauli principle
 
      phase[ i ] = 0.0; 
 
      //std::cout << "TKDB Check Pauli Principle " << i << std::endl;
 
      for ( G4int j = 0 ; j < i ; j++ )
      {
         phase[ j ] = 0.0;
         //std::cout << "TKDB Check Pauli Principle  i , j " << i << " , " << j << std::endl;
         G4double expa = 0.0;
         if ( participants[j]->GetDefinition() ==  participants[i]->GetDefinition() )
         {
 
            expa = - meanfield->GetRR2(i,j) * cpw;
 
            if ( expa > epsx ) 
            {
               G4ThreeVector p_i = participants[i]->GetMomentum();  
               G4ThreeVector pj = participants[j]->GetMomentum();  
               G4double dist2_p = p_i.diff2( pj ); 
 
               dist2_p = dist2_p*cph;
               expa = expa - dist2_p; 
 
               if ( expa > epsx ) 
               {
 
                  phase[j] = G4Exp ( expa );
 
                  if ( phase[j] * cpc > 0.2 ) 
                  { 
/*
         G4cout << "TKDB Check Pauli Principle A i , j " << i << " , " << j << G4endl;
         G4cout << "TKDB Check Pauli Principle phase[j] " << phase[j] << G4endl;
         G4cout << "TKDB Check Pauli Principle phase[j]*cpc > 0.2 " << phase[j]*cpc << G4endl;
*/
                     isThisOK = false;
                     break;
                  }
                  if ( ( phase_g[j] + phase[j] ) * cpc > 0.5 ) 
                  { 
/*
         G4cout << "TKDB Check Pauli Principle B i , j " << i << " , " << j << G4endl;
         G4cout << "TKDB Check Pauli Principle B phase_g[j] " << phase_g[j] << G4endl;
         G4cout << "TKDB Check Pauli Principle B phase[j] " << phase[j] << G4endl;
         G4cout << "TKDB Check Pauli Principle B phase_g[j] + phase[j] ) * cpc > 0.5  " <<  ( phase_g[j] + phase[j] ) * cpc << G4endl;
*/
                     isThisOK = false;
                     break;
                  }
 
                  phase[i] += phase[j];
                  if ( phase[i] * cpc > 0.3 ) 
                  { 
/*
         G4cout << "TKDB Check Pauli Principle C i , j " << i << " , " << j << G4endl;
         G4cout << "TKDB Check Pauli Principle C phase[i] " << phase[i] << G4endl;
         G4cout << "TKDB Check Pauli Principle C phase[i] * cpc > 0.3 " <<  phase[i] * cpc << G4endl;
*/
                     isThisOK = false;
                     break;
                  }
 
                  //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl;
 
               }
               else
               {
                  //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl;
               }
 
            }
            else
            {
               //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl;
            }
 
         }
         else
         {
            //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl;
         }
 
      }
 
      if ( isThisOK == true )
      {
 
         phase_g[i] = phase[i];
 
         for ( int j = 0 ; j < i ; j++ )
         {
            phase_g[j] += phase[j]; 
         }
 
         result = true; 
         return result;
      }
 
   }
 
   return result;
 
}

References G4Pow::A13(), cpc, cph, cpw, d_pot, DBL_MAX, CLHEP::Hep3Vector::diff2(), ebini, edepth, epsx, G4cout, G4endl, G4Exp(), G4UniformRand, G4Pow::GetInstance(), G4QMDNucleus::GetMassNumber(), G4QMDMeanField::GetRR2(), hbc, maxTrial, meanfield, G4QMDSystem::participants, phase_g, pi, and rho_l.

Referenced by packNucleons().

◆ samplingPosition()

G4bool G4QMDGroundStateNucleus::samplingPosition ( G4int i )

private

Definition at line 435 of file G4QMDGroundStateNucleus.cc.

{
 
   G4bool result = false;
 
   G4int nTry = 0; 
   
   while ( nTry < maxTrial )  // Loop checking, 11.03.2015, T. Koi
   {
 
      //std::cout << "samplingPosition i th particle, nTtry " << i << " " << nTry << std::endl;  
 
      G4double rwod = -1.0;        
      G4double rrr = 0.0; 
 
      G4double rx = 0.0;
      G4double ry = 0.0;
      G4double rz = 0.0;
 
      G4int icounter = 0;
      G4int icounter_max = 1024;
      while ( G4UniformRand() * rmax > rwod ) // Loop checking, 11.03.2015, T. Koi
      {
         icounter++;
         if ( icounter > icounter_max ) {
        G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
            break;
         }
 
         G4double rsqr = 10.0; 
         G4int jcounter = 0;
         G4int jcounter_max = 1024;
         while ( rsqr > 1.0 ) // Loop checking, 11.03.2015, T. Koi
         {
            jcounter++;
            if ( jcounter > jcounter_max ) {
           G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
               break;
            }
            rx = 1.0 - 2.0 * G4UniformRand();
            ry = 1.0 - 2.0 * G4UniformRand();
            rz = 1.0 - 2.0 * G4UniformRand();
            rsqr = rx*rx + ry*ry + rz*rz; 
         }
         rrr = radm * std::sqrt ( rsqr );
         rwod = 1.0 / ( 1.0 + G4Exp ( ( rrr - rt00 ) / saa ) );
 
      } 
 
      participants[i]->SetPosition( G4ThreeVector( rx , ry , rz )*radm );
 
      if ( i == 0 )   
      {
         result = true; 
         return result;
      }
 
//    i > 1 ( Second Particle or later )
//    Check Distance to others 
 
      G4bool isThisOK = true;
      for ( G4int j = 0 ; j < i ; j++ )
      {
 
         G4double r2 =  participants[j]->GetPosition().diff2( participants[i]->GetPosition() );   
         G4double dmin2 = 0.0;
 
         if ( participants[j]->GetDefinition() == participants[i]->GetDefinition() )
         {
            dmin2 = dsam2;
         }
         else
         {
            dmin2 = ddif2;
         }
 
         //std::cout << "distance between j and i " << j << " " << i << " " << r2 << " " << dmin2 << std::endl;
         if ( r2 < dmin2 )
         {
            isThisOK = false;
            break; 
         }
 
      }
 
      if ( isThisOK == true )
      {
         result = true;
         return result; 
      }
 
      nTry++; 
 
   }
 
// Here return "false" 
   return result;
}

References ddif2, dsam2, G4cout, G4endl, G4Exp(), G4UniformRand, maxTrial, G4QMDSystem::participants, radm, rmax, rt00, and saa.

Referenced by packNucleons().

◆ SetParticipant()

void G4QMDSystem::SetParticipant ( G4QMDParticipant * particle )

inlineinherited

Definition at line 51 of file G4QMDSystem.hh.

51{ participants.push_back ( particle ); };

References G4QMDSystem::participants.

Referenced by G4QMDReaction::ApplyYourself(), G4QMDCollision::CalKinematicsOfBinaryCollisions(), G4QMDMeanField::DoClusterJudgment(), G4QMDGroundStateNucleus(), and G4QMDSystem::SetSystem().

◆ SetSystem()

void G4QMDSystem::SetSystem	(	G4QMDSystem *	nucleus,
		G4ThreeVector	dp,
		G4ThreeVector	dr
	)

inherited

Definition at line 46 of file G4QMDSystem.cc.

{
   std::vector< G4QMDParticipant* >::iterator it; 
   for ( it = nucleus->participants.begin() ; it != nucleus->participants.end() ; it++ ) 
   {
      G4ThreeVector r = (*it)->GetPosition() + dr;
      (*it)->SetPosition ( r );
      G4ThreeVector p = (*it)->GetMomentum() + dp;
      (*it)->SetMomentum ( p );
      this->SetParticipant( *it );
   }
}

References G4QMDSystem::participants, and G4QMDSystem::SetParticipant().

◆ SetTotalPotential()

void G4QMDNucleus::SetTotalPotential ( G4double x )

inlineinherited

Definition at line 62 of file G4QMDNucleus.hh.

62{ potentialEnergy = x; };

References G4QMDNucleus::potentialEnergy.

Referenced by G4QMDReaction::ApplyYourself(), and G4QMDMeanField::SetNucleus().

◆ ShowParticipants()

void G4QMDSystem::ShowParticipants ( )

inherited

Definition at line 79 of file G4QMDSystem.cc.

{
   //store orginal precision
   std::ios::fmtflags oldform = G4cout.flags();
 
   G4ThreeVector p_sum( 0.0 );
   std::vector< G4QMDParticipant* >::iterator it; 
   G4cout << "Momentum and Position of each participant " << G4endl; 
   G4int i = 0; 
   for ( it = participants.begin() ; it != participants.end() ; it++ ) 
   { 
      G4cout << i
             << " " 
             << (*it)->GetDefinition()->GetParticleName() 
             << " " 
             << std::setprecision( 8 )
             << (*it)->GetMomentum()
             << " " 
             << (*it)->GetPosition() 
             << G4endl;
      p_sum += (*it)->GetMomentum();
      i++;
   }
   G4cout << "Sum upped Momentum and its mag " << p_sum << " " << p_sum.mag() << G4endl;
 
   //restore orginal precision
   G4cout.flags( oldform );
}

References G4cout, G4endl, CLHEP::Hep3Vector::mag(), and G4QMDSystem::participants.

◆ SubtractSystem()

void G4QMDSystem::SubtractSystem ( G4QMDSystem * nucleus )

inherited

Definition at line 59 of file G4QMDSystem.cc.

{
   
   for ( G4int i = 0 ; i < nucleus->GetTotalNumberOfParticipant() ; i++ )
   {
      participants.erase ( std::find ( participants.begin() , participants.end() , nucleus->GetParticipant( i ) ) );
   }
}