#include <G4StatMF.hh>
Inheritance diagram for G4StatMF:
Public Member Functions | |
G4StatMF () | |
~G4StatMF () | |
G4FragmentVector * | BreakItUp (const G4Fragment &theNucleus) |
Definition at line 47 of file G4StatMF.hh.
G4StatMF::G4StatMF | ( | ) |
G4StatMF::~G4StatMF | ( | ) |
G4FragmentVector * G4StatMF::BreakItUp | ( | const G4Fragment & | theNucleus | ) | [virtual] |
Implements G4VMultiFragmentation.
Definition at line 46 of file G4StatMF.cc.
References G4StatMFChannel::CheckFragments(), G4StatMFMacroCanonical::ChooseAandZ(), G4StatMFMicroCanonical::ChooseAandZ(), G4Fragment::GetA(), G4Fragment::GetA_asInt(), G4Fragment::GetExcitationEnergy(), G4StatMFParameters::GetMaxAverageMultiplicity(), G4VStatMFEnsemble::GetMeanMultiplicity(), G4VStatMFEnsemble::GetMeanTemperature(), G4Fragment::GetMomentum(), G4StatMFChannel::GetMultiplicity(), G4Fragment::GetZ_asInt(), and G4Chips::Temperature.
00047 { 00048 // G4FragmentVector * theResult = new G4FragmentVector; 00049 00050 if (theFragment.GetExcitationEnergy() <= 0.0) { 00051 //G4FragmentVector * theResult = new G4FragmentVector; 00052 //theResult->push_back(new G4Fragment(theFragment)); 00053 return 0; 00054 } 00055 00056 00057 // Maximun average multiplicity: M_0 = 2.6 for A ~ 200 00058 // and M_0 = 3.3 for A <= 110 00059 G4double MaxAverageMultiplicity = 00060 G4StatMFParameters::GetMaxAverageMultiplicity(static_cast<G4int>(theFragment.GetA())); 00061 00062 00063 // We'll use two kinds of ensembles 00064 G4StatMFMicroCanonical * theMicrocanonicalEnsemble = 0; 00065 G4StatMFMacroCanonical * theMacrocanonicalEnsemble = 0; 00066 00067 00068 //------------------------------------------------------- 00069 // Direct simulation part (Microcanonical ensemble) 00070 //------------------------------------------------------- 00071 00072 // Microcanonical ensemble initialization 00073 theMicrocanonicalEnsemble = new G4StatMFMicroCanonical(theFragment); 00074 00075 G4int Iterations = 0; 00076 G4int IterationsLimit = 100000; 00077 G4double Temperature = 0.0; 00078 00079 G4bool FirstTime = true; 00080 G4StatMFChannel * theChannel = 0; 00081 00082 G4bool ChannelOk; 00083 do { // Try to de-excite as much as IterationLimit permits 00084 do { 00085 00086 G4double theMeanMult = theMicrocanonicalEnsemble->GetMeanMultiplicity(); 00087 if (theMeanMult <= MaxAverageMultiplicity) { 00088 // G4cout << "MICROCANONICAL" << G4endl; 00089 // Choose fragments atomic numbers and charges from direct simulation 00090 theChannel = theMicrocanonicalEnsemble->ChooseAandZ(theFragment); 00091 _theEnsemble = theMicrocanonicalEnsemble; 00092 } else { 00093 //----------------------------------------------------- 00094 // Non direct simulation part (Macrocanonical Ensemble) 00095 //----------------------------------------------------- 00096 if (FirstTime) { 00097 // Macrocanonical ensemble initialization 00098 theMacrocanonicalEnsemble = new G4StatMFMacroCanonical(theFragment); 00099 _theEnsemble = theMacrocanonicalEnsemble; 00100 FirstTime = false; 00101 } 00102 // G4cout << "MACROCANONICAL" << G4endl; 00103 // Select calculated fragment total multiplicity, 00104 // fragment atomic numbers and fragment charges. 00105 theChannel = theMacrocanonicalEnsemble->ChooseAandZ(theFragment); 00106 } 00107 00108 ChannelOk = theChannel->CheckFragments(); 00109 if (!ChannelOk) delete theChannel; 00110 00111 } while (!ChannelOk); 00112 00113 00114 if (theChannel->GetMultiplicity() <= 1) { 00115 G4FragmentVector * theResult = new G4FragmentVector; 00116 theResult->push_back(new G4Fragment(theFragment)); 00117 delete theMicrocanonicalEnsemble; 00118 if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble; 00119 delete theChannel; 00120 return theResult; 00121 } 00122 00123 //-------------------------------------- 00124 // Second part of simulation procedure. 00125 //-------------------------------------- 00126 00127 // Find temperature of breaking channel. 00128 Temperature = _theEnsemble->GetMeanTemperature(); // Initial guess for Temperature 00129 00130 if (FindTemperatureOfBreakingChannel(theFragment,theChannel,Temperature)) break; 00131 00132 // Do not forget to delete this unusable channel, for which we failed to find the temperature, 00133 // otherwise for very proton-reach nuclei it would lead to memory leak due to large 00134 // number of iterations. N.B. "theChannel" is created in G4StatMFMacroCanonical::ChooseZ() 00135 00136 // G4cout << " Iteration # " << Iterations << " Mean Temperature = " << Temperature << G4endl; 00137 00138 delete theChannel; 00139 00140 } while (Iterations++ < IterationsLimit ); 00141 00142 00143 00144 // If Iterations >= IterationsLimit means that we couldn't solve for temperature 00145 if (Iterations >= IterationsLimit) 00146 throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::BreakItUp: Was not possible to solve for temperature of breaking channel"); 00147 00148 00149 G4FragmentVector * theResult = theChannel-> 00150 GetFragments(theFragment.GetA_asInt(),theFragment.GetZ_asInt(),Temperature); 00151 00152 00153 00154 // ~~~~~~ Energy conservation Patch !!!!!!!!!!!!!!!!!!!!!! 00155 // Original nucleus 4-momentum in CM system 00156 G4LorentzVector InitialMomentum(theFragment.GetMomentum()); 00157 InitialMomentum.boost(-InitialMomentum.boostVector()); 00158 G4double ScaleFactor = 0.0; 00159 G4double SavedScaleFactor = 0.0; 00160 do { 00161 G4double FragmentsEnergy = 0.0; 00162 G4FragmentVector::iterator j; 00163 for (j = theResult->begin(); j != theResult->end(); j++) 00164 FragmentsEnergy += (*j)->GetMomentum().e(); 00165 SavedScaleFactor = ScaleFactor; 00166 ScaleFactor = InitialMomentum.e()/FragmentsEnergy; 00167 G4ThreeVector ScaledMomentum(0.0,0.0,0.0); 00168 for (j = theResult->begin(); j != theResult->end(); j++) { 00169 ScaledMomentum = ScaleFactor * (*j)->GetMomentum().vect(); 00170 G4double Mass = (*j)->GetMomentum().m(); 00171 G4LorentzVector NewMomentum; 00172 NewMomentum.setVect(ScaledMomentum); 00173 NewMomentum.setE(std::sqrt(ScaledMomentum.mag2()+Mass*Mass)); 00174 (*j)->SetMomentum(NewMomentum); 00175 } 00176 } while (ScaleFactor > 1.0+1.e-5 && std::abs(ScaleFactor-SavedScaleFactor)/ScaleFactor > 1.e-10); 00177 // ~~~~~~ End of patch !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00178 00179 // Perform Lorentz boost 00180 G4FragmentVector::iterator i; 00181 for (i = theResult->begin(); i != theResult->end(); i++) { 00182 G4LorentzVector FourMom = (*i)->GetMomentum(); 00183 FourMom.boost(theFragment.GetMomentum().boostVector()); 00184 (*i)->SetMomentum(FourMom); 00185 #ifdef PRECOMPOUND_TEST 00186 (*i)->SetCreatorModel(G4String("G4StatMF")); 00187 #endif 00188 } 00189 00190 // garbage collection 00191 delete theMicrocanonicalEnsemble; 00192 if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble; 00193 delete theChannel; 00194 00195 return theResult; 00196 }