Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4StatMF Class Reference

#include <G4StatMF.hh>

+ Inheritance diagram for G4StatMF:

Public Member Functions

 G4StatMF ()
 
 ~G4StatMF ()
 
 G4StatMF (const G4StatMF &right)=delete
 
G4StatMFoperator= (const G4StatMF &right)=delete
 
G4bool operator== (const G4StatMF &right)=delete
 
G4bool operator!= (const G4StatMF &right)=delete
 
G4FragmentVectorBreakItUp (const G4Fragment &theNucleus) override
 
- Public Member Functions inherited from G4VMultiFragmentation
 G4VMultiFragmentation ()
 
virtual ~G4VMultiFragmentation ()
 

Detailed Description

Definition at line 45 of file G4StatMF.hh.

Constructor & Destructor Documentation

◆ G4StatMF() [1/2]

G4StatMF::G4StatMF ( )

Definition at line 37 of file G4StatMF.cc.

38{
39 _secID = G4PhysicsModelCatalog::GetModelID("model_G4StatMF");
40}
static G4int GetModelID(const G4int modelIndex)

◆ ~G4StatMF()

G4StatMF::~G4StatMF ( )

Definition at line 42 of file G4StatMF.cc.

42{}

◆ G4StatMF() [2/2]

G4StatMF::G4StatMF ( const G4StatMF & right)
delete

Member Function Documentation

◆ BreakItUp()

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

Implements G4VMultiFragmentation.

Definition at line 44 of file G4StatMF.cc.

45{
46 if (theFragment.GetExcitationEnergy() <= 0.0) {
47 return nullptr;
48 }
49
50 // Maximun average multiplicity: M_0 = 2.6 for A ~ 200
51 // and M_0 = 3.3 for A <= 110
52 G4double MaxAverageMultiplicity =
53 G4StatMFParameters::GetMaxAverageMultiplicity(theFragment.GetA_asInt());
54
55
56 // We'll use two kinds of ensembles
57 G4StatMFMicroCanonical * theMicrocanonicalEnsemble = 0;
58 G4StatMFMacroCanonical * theMacrocanonicalEnsemble = 0;
59
60 //-------------------------------------------------------
61 // Direct simulation part (Microcanonical ensemble)
62 //-------------------------------------------------------
63
64 // Microcanonical ensemble initialization
65 theMicrocanonicalEnsemble = new G4StatMFMicroCanonical(theFragment);
66
67 G4int Iterations = 0;
68 G4int IterationsLimit = 100000;
69 G4double Temperature = 0.0;
70
71 G4bool FirstTime = true;
72 G4StatMFChannel * theChannel = 0;
73
74 G4bool ChannelOk;
75 do { // Try to de-excite as much as IterationLimit permits
76 do {
77
78 G4double theMeanMult = theMicrocanonicalEnsemble->GetMeanMultiplicity();
79 if (theMeanMult <= MaxAverageMultiplicity) {
80 // G4cout << "MICROCANONICAL" << G4endl;
81 // Choose fragments atomic numbers and charges from direct simulation
82 theChannel = theMicrocanonicalEnsemble->ChooseAandZ(theFragment);
83 _theEnsemble = theMicrocanonicalEnsemble;
84 } else {
85 //-----------------------------------------------------
86 // Non direct simulation part (Macrocanonical Ensemble)
87 //-----------------------------------------------------
88 if (FirstTime) {
89 // Macrocanonical ensemble initialization
90 theMacrocanonicalEnsemble = new G4StatMFMacroCanonical(theFragment);
91 _theEnsemble = theMacrocanonicalEnsemble;
92 FirstTime = false;
93 }
94 // G4cout << "MACROCANONICAL" << G4endl;
95 // Select calculated fragment total multiplicity,
96 // fragment atomic numbers and fragment charges.
97 theChannel = theMacrocanonicalEnsemble->ChooseAandZ(theFragment);
98 }
99
100 ChannelOk = theChannel->CheckFragments();
101 if (!ChannelOk) delete theChannel;
102
103 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
104 } while (!ChannelOk);
105
106
107 if (theChannel->GetMultiplicity() <= 1) {
108 G4FragmentVector * theResult = new G4FragmentVector;
109 theResult->push_back(new G4Fragment(theFragment));
110 delete theMicrocanonicalEnsemble;
111 if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
112 delete theChannel;
113 return theResult;
114 }
115
116 //--------------------------------------
117 // Second part of simulation procedure.
118 //--------------------------------------
119
120 // Find temperature of breaking channel.
121 Temperature = _theEnsemble->GetMeanTemperature(); // Initial guess for Temperature
122
123 if (FindTemperatureOfBreakingChannel(theFragment,theChannel,Temperature)) break;
124
125 // Do not forget to delete this unusable channel, for which we failed to find the temperature,
126 // otherwise for very proton-reach nuclei it would lead to memory leak due to large
127 // number of iterations. N.B. "theChannel" is created in G4StatMFMacroCanonical::ChooseZ()
128
129 // G4cout << " Iteration # " << Iterations << " Mean Temperature = " << Temperature << G4endl;
130
131 delete theChannel;
132
133 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
134 } while (Iterations++ < IterationsLimit );
135
136 // If Iterations >= IterationsLimit means that we couldn't solve for temperature
137 if (Iterations >= IterationsLimit)
138 throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::BreakItUp: Was not possible to solve for temperature of breaking channel");
139
140 G4FragmentVector * theResult = theChannel->
141 GetFragments(theFragment.GetA_asInt(),theFragment.GetZ_asInt(),Temperature);
142
143 // ~~~~~~ Energy conservation Patch !!!!!!!!!!!!!!!!!!!!!!
144 // Original nucleus 4-momentum in CM system
145 G4LorentzVector InitialMomentum(theFragment.GetMomentum());
146 InitialMomentum.boost(-InitialMomentum.boostVector());
147 G4double ScaleFactor = 0.0;
148 G4double SavedScaleFactor = 0.0;
149 do {
150 G4double FragmentsEnergy = 0.0;
151 G4FragmentVector::iterator j;
152 for (j = theResult->begin(); j != theResult->end(); j++)
153 FragmentsEnergy += (*j)->GetMomentum().e();
154 SavedScaleFactor = ScaleFactor;
155 ScaleFactor = InitialMomentum.e()/FragmentsEnergy;
156 G4ThreeVector ScaledMomentum(0.0,0.0,0.0);
157 for (j = theResult->begin(); j != theResult->end(); j++) {
158 ScaledMomentum = ScaleFactor * (*j)->GetMomentum().vect();
159 G4double Mass = (*j)->GetMomentum().m();
160 G4LorentzVector NewMomentum;
161 NewMomentum.setVect(ScaledMomentum);
162 NewMomentum.setE(std::sqrt(ScaledMomentum.mag2()+Mass*Mass));
163 (*j)->SetMomentum(NewMomentum);
164 }
165 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
166 } while (ScaleFactor > 1.0+1.e-5 && std::abs(ScaleFactor-SavedScaleFactor)/ScaleFactor > 1.e-10);
167 // ~~~~~~ End of patch !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
168
169 // Perform Lorentz boost
170 G4FragmentVector::iterator i;
171 for (i = theResult->begin(); i != theResult->end(); i++) {
172 G4LorentzVector FourMom = (*i)->GetMomentum();
173 FourMom.boost(theFragment.GetMomentum().boostVector());
174 (*i)->SetMomentum(FourMom);
175 (*i)->SetCreatorModelID(_secID);
176 }
177
178 // garbage collection
179 delete theMicrocanonicalEnsemble;
180 if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
181 delete theChannel;
182
183 return theResult;
184}
std::vector< G4Fragment * > G4FragmentVector
Definition G4Fragment.hh:65
double G4double
Definition G4Types.hh:83
bool G4bool
Definition G4Types.hh:86
int G4int
Definition G4Types.hh:85
HepLorentzVector & boost(double, double, double)
void setVect(const Hep3Vector &)
G4bool CheckFragments(void)
size_t GetMultiplicity(void)
G4StatMFChannel * ChooseAandZ(const G4Fragment &theFragment)
G4StatMFChannel * ChooseAandZ(const G4Fragment &theFragment)
static G4double GetMaxAverageMultiplicity(G4int A)
G4double GetMeanTemperature(void) const
G4double GetMeanMultiplicity(void) const

◆ operator!=()

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

◆ operator=()

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

◆ operator==()

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

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