G4StatMF Class Reference

#include <G4StatMF.hh>

Inheritance diagram for G4StatMF:

Public Member Functions
	G4StatMF ()
	~G4StatMF ()
G4FragmentVector *	BreakItUp (const G4Fragment &theNucleus)
Public Member Functions inherited from G4VMultiFragmentation
	G4VMultiFragmentation ()
virtual	~G4VMultiFragmentation ()

Detailed Description

Definition at line 47 of file G4StatMF.hh.

Constructor & Destructor Documentation

G4StatMF::G4StatMF

(

)

Definition at line 39 of file G4StatMF.cc.

: _theEnsemble(0) {}

G4StatMF::~G4StatMF

(

)

Definition at line 43 of file G4StatMF.cc.

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

Member Function Documentation

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

virtual

Implements G4VMultiFragmentation.

Definition at line 46 of file G4StatMF.cc.

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

{
  //    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; 
      
    } 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;    

  } 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);       
    }
  } 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);
#ifdef PRECOMPOUND_TEST
    (*i)->SetCreatorModel(G4String("G4StatMF"));
#endif
  }
  
  // garbage collection
  delete theMicrocanonicalEnsemble;
  if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
  delete theChannel;
  
  return theResult;
}

---

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

• G4StatMF.hh
• G4StatMF.cc