G4StatMF Class Reference

#include <G4StatMF.hh>

Inheritance diagram for G4StatMF:

Public Member Functions
G4FragmentVector *	BreakItUp (const G4Fragment &theNucleus)
	G4StatMF ()
	~G4StatMF ()

Private Member Functions
G4double	CalcEnergy (G4int A, G4int Z, const G4StatMFChannel *aChannel, G4double T)
G4bool	FindTemperatureOfBreakingChannel (const G4Fragment &theFragment, const G4StatMFChannel *aChannel, G4double &Temperature)
	G4StatMF (const G4StatMF &right)
G4bool	operator!= (const G4StatMF &right)
G4StatMF &	operator= (const G4StatMF &right)
G4bool	operator== (const G4StatMF &right)

Private Attributes
G4int	_secID
G4VStatMFEnsemble *	_theEnsemble

Detailed Description

Definition at line 46 of file G4StatMF.hh.

Constructor & Destructor Documentation

G4StatMF() [1/2]

G4StatMF::G4StatMF

(

)

Definition at line 38 of file G4StatMF.cc.

                   : _theEnsemble(0), _secID(-1) {
  _secID = G4PhysicsModelCatalog::GetModelID("model_G4StatMF");
}

References _secID, and G4PhysicsModelCatalog::GetModelID().

~G4StatMF()

G4StatMF::~G4StatMF

(

)

Definition at line 44 of file G4StatMF.cc.

{} //{if (_theEnsemble != 0) delete _theEnsemble;}

G4StatMF() [2/2]

G4StatMF::G4StatMF ( const G4StatMF &  right )

private

Member Function Documentation

BreakItUp()

G4FragmentVector * G4StatMF::BreakItUp ( const G4Fragment &  theNucleus )

virtual

Implements G4VMultiFragmentation.

Definition at line 47 of file G4StatMF.cc.

{
  //    G4FragmentVector * theResult = new G4FragmentVector;
 
  if (theFragment.GetExcitationEnergy() <= 0.0) {
    //G4FragmentVector * theResult = new G4FragmentVector;
    //theResult->push_back(new G4Fragment(theFragment));
    return 0;
  }
 
 
  // Maximun average multiplicity: M_0 = 2.6 for A ~ 200 
  // and M_0 = 3.3 for A <= 110
  G4double MaxAverageMultiplicity = 
    G4StatMFParameters::GetMaxAverageMultiplicity(theFragment.GetA_asInt());
 
    
    // We'll use two kinds of ensembles
  G4StatMFMicroCanonical * theMicrocanonicalEnsemble = 0;
  G4StatMFMacroCanonical * theMacrocanonicalEnsemble = 0;
    
  //-------------------------------------------------------
  // Direct simulation part (Microcanonical ensemble)
  //-------------------------------------------------------
  
  // Microcanonical ensemble initialization 
  theMicrocanonicalEnsemble = new G4StatMFMicroCanonical(theFragment);
 
  G4int Iterations = 0;
  G4int IterationsLimit = 100000;
  G4double Temperature = 0.0;
  
  G4bool FirstTime = true;
  G4StatMFChannel * theChannel = 0;
 
  G4bool ChannelOk;
  do {  // Try to de-excite as much as IterationLimit permits
    do {
      
      G4double theMeanMult = theMicrocanonicalEnsemble->GetMeanMultiplicity();
      if (theMeanMult <= MaxAverageMultiplicity) {
    // G4cout << "MICROCANONICAL" << G4endl;
    // Choose fragments atomic numbers and charges from direct simulation
    theChannel = theMicrocanonicalEnsemble->ChooseAandZ(theFragment);
    _theEnsemble = theMicrocanonicalEnsemble;
      } else {
    //-----------------------------------------------------
    // Non direct simulation part (Macrocanonical Ensemble)
    //-----------------------------------------------------
    if (FirstTime) {
      // Macrocanonical ensemble initialization 
      theMacrocanonicalEnsemble = new G4StatMFMacroCanonical(theFragment);
      _theEnsemble = theMacrocanonicalEnsemble;
      FirstTime = false;
    }
    // G4cout << "MACROCANONICAL" << G4endl;
    // Select calculated fragment total multiplicity, 
    // fragment atomic numbers and fragment charges.
    theChannel = theMacrocanonicalEnsemble->ChooseAandZ(theFragment);
      }
      
      ChannelOk = theChannel->CheckFragments();
      if (!ChannelOk) delete theChannel; 
      
      // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
    } while (!ChannelOk);
    
    
    if (theChannel->GetMultiplicity() <= 1) {
      G4FragmentVector * theResult = new G4FragmentVector;
      theResult->push_back(new G4Fragment(theFragment));
      delete theMicrocanonicalEnsemble;
      if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
      delete theChannel;
      return theResult;
    }
    
    //--------------------------------------
    // Second part of simulation procedure.
    //--------------------------------------
    
    // Find temperature of breaking channel.
    Temperature = _theEnsemble->GetMeanTemperature(); // Initial guess for Temperature 
 
    if (FindTemperatureOfBreakingChannel(theFragment,theChannel,Temperature)) break;
 
    // Do not forget to delete this unusable channel, for which we failed to find the temperature,
    // otherwise for very proton-reach nuclei it would lead to memory leak due to large 
    // number of iterations. N.B. "theChannel" is created in G4StatMFMacroCanonical::ChooseZ()
 
    // G4cout << " Iteration # " << Iterations << " Mean Temperature = " << Temperature << G4endl;    
 
    delete theChannel;    
 
    // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
  } while (Iterations++ < IterationsLimit );
  
 
 
  // If Iterations >= IterationsLimit means that we couldn't solve for temperature
  if (Iterations >= IterationsLimit) 
    throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::BreakItUp: Was not possible to solve for temperature of breaking channel");
  
  
  G4FragmentVector * theResult = theChannel->
    GetFragments(theFragment.GetA_asInt(),theFragment.GetZ_asInt(),Temperature);
  
  
    
  // ~~~~~~ Energy conservation Patch !!!!!!!!!!!!!!!!!!!!!!
  // Original nucleus 4-momentum in CM system
  G4LorentzVector InitialMomentum(theFragment.GetMomentum());
  InitialMomentum.boost(-InitialMomentum.boostVector());
  G4double ScaleFactor = 0.0;
  G4double SavedScaleFactor = 0.0;
  do {
    G4double FragmentsEnergy = 0.0;
    G4FragmentVector::iterator j;
    for (j = theResult->begin(); j != theResult->end(); j++) 
      FragmentsEnergy += (*j)->GetMomentum().e();
    SavedScaleFactor = ScaleFactor;
    ScaleFactor = InitialMomentum.e()/FragmentsEnergy;
    G4ThreeVector ScaledMomentum(0.0,0.0,0.0);
    for (j = theResult->begin(); j != theResult->end(); j++) {
      ScaledMomentum = ScaleFactor * (*j)->GetMomentum().vect();
      G4double Mass = (*j)->GetMomentum().m();
      G4LorentzVector NewMomentum;
      NewMomentum.setVect(ScaledMomentum);
      NewMomentum.setE(std::sqrt(ScaledMomentum.mag2()+Mass*Mass));
      (*j)->SetMomentum(NewMomentum);       
    }
    // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
  } while (ScaleFactor > 1.0+1.e-5 && std::abs(ScaleFactor-SavedScaleFactor)/ScaleFactor > 1.e-10);
  // ~~~~~~ End of patch !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  // Perform Lorentz boost
  G4FragmentVector::iterator i;
  for (i = theResult->begin(); i != theResult->end(); i++) {
    G4LorentzVector FourMom = (*i)->GetMomentum();
    FourMom.boost(theFragment.GetMomentum().boostVector());
    (*i)->SetMomentum(FourMom);
    (*i)->SetCreatorModelID(_secID);
  }
  
  // garbage collection
  delete theMicrocanonicalEnsemble;
  if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
  delete theChannel;
  
  return theResult;
}

References _secID, _theEnsemble, CLHEP::HepLorentzVector::boost(), CLHEP::HepLorentzVector::boostVector(), G4StatMFChannel::CheckFragments(), G4StatMFMacroCanonical::ChooseAandZ(), G4StatMFMicroCanonical::ChooseAandZ(), CLHEP::HepLorentzVector::e(), FindTemperatureOfBreakingChannel(), G4Fragment::GetA_asInt(), G4Fragment::GetExcitationEnergy(), G4StatMFParameters::GetMaxAverageMultiplicity(), G4VStatMFEnsemble::GetMeanMultiplicity(), G4VStatMFEnsemble::GetMeanTemperature(), G4Fragment::GetMomentum(), G4StatMFChannel::GetMultiplicity(), G4Fragment::GetZ_asInt(), CLHEP::Hep3Vector::mag2(), CLHEP::HepLorentzVector::setE(), and CLHEP::HepLorentzVector::setVect().

CalcEnergy()

G4double G4StatMF::CalcEnergy ( G4int  A, G4int  Z, const G4StatMFChannel *  aChannel, G4double  T  )

private

Definition at line 277 of file G4StatMF.cc.

{
  G4double MassExcess0 = G4NucleiProperties::GetMassExcess(A,Z);
  G4double ChannelEnergy = aChannel->GetFragmentsEnergy(T);
  return -MassExcess0 + G4StatMFParameters::GetCoulomb() + ChannelEnergy;
}

References A, G4StatMFParameters::GetCoulomb(), G4StatMFChannel::GetFragmentsEnergy(), G4NucleiProperties::GetMassExcess(), and Z.

Referenced by FindTemperatureOfBreakingChannel().

FindTemperatureOfBreakingChannel()

G4bool G4StatMF::FindTemperatureOfBreakingChannel ( const G4Fragment &  theFragment, const G4StatMFChannel *  aChannel, G4double &  Temperature  )

private

Definition at line 200 of file G4StatMF.cc.

{
  G4int A = theFragment.GetA_asInt();
  G4int Z = theFragment.GetZ_asInt();
  G4double U = theFragment.GetExcitationEnergy();
  
  G4double T = std::max(Temperature,0.0012*MeV);  
  G4double Ta = T;
  G4double TotalEnergy = CalcEnergy(A,Z,aChannel,T);
  
  G4double Da = (U - TotalEnergy)/U;
  G4double Db = 0.0;
  
  // bracketing the solution
  if (Da == 0.0) {
    Temperature = T;
    return true;
  } else if (Da < 0.0) {
    do {
      T *= 0.5;
      if (T < 0.001*MeV) return false;
      
      TotalEnergy = CalcEnergy(A,Z,aChannel,T);
      
      Db = (U - TotalEnergy)/U;
      // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
    } while (Db < 0.0);
    
  } else {
    do {
      T *= 1.5;
      
      TotalEnergy = CalcEnergy(A,Z,aChannel,T);
      
      Db = (U - TotalEnergy)/U;
      // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
    } while (Db > 0.0); 
  }
  
  G4double eps = 1.0e-14 * std::abs(T-Ta);
  //G4double eps = 1.0e-3 ;
  
  // Start the bisection method
  for (G4int j = 0; j < 1000; j++) {
    G4double Tc =  (Ta+T)*0.5;
    if (std::abs(Ta-Tc) <= eps) {
      Temperature = Tc;
      return true;
    }
    
    T = Tc;
    
    TotalEnergy = CalcEnergy(A,Z,aChannel,T);
    
    G4double Dc = (U - TotalEnergy)/U; 
    
    if (Dc == 0.0) {
      Temperature  = Tc;
      return true;
    }
    
    if (Da*Dc < 0.0) {
      T  = Tc;
      Db = Dc;
    } else {
      Ta = Tc;
      Da = Dc;
    }
  }
  
  Temperature  = (Ta+T)*0.5;
  return false;
}

References A, CalcEnergy(), eps, G4Fragment::GetA_asInt(), G4Fragment::GetExcitationEnergy(), G4Fragment::GetZ_asInt(), G4INCL::Math::max(), MeV, and Z.

Referenced by BreakItUp().

operator!=()

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

private

operator=()

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

private

operator==()

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

private

Field Documentation

_secID

G4int G4StatMF::_secID

private

Definition at line 84 of file G4StatMF.hh.

Referenced by BreakItUp(), and G4StatMF().

_theEnsemble

G4VStatMFEnsemble* G4StatMF::_theEnsemble

private

Definition at line 82 of file G4StatMF.hh.

Referenced by BreakItUp().

---

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

• source/processes/hadronic/models/de_excitation/multifragmentation/include/G4StatMF.hh
• source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMF.cc