| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara
30
31 #include "G4StatMF.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4Pow.hh"
35 #include "G4PhysicsModelCatalog.hh"
36
37 G4StatMF::G4StatMF()
38 {
39 _secID = G4PhysicsModelCatalog::GetModelID("
40 }
41
42 G4StatMF::~G4StatMF() {}
43
44 G4FragmentVector* G4StatMF::BreakItUp(const G4
45 {
46 if (theFragment.GetExcitationEnergy() <= 0.0
47 return nullptr;
48 }
49
50 // Maximun average multiplicity: M_0 = 2.6 f
51 // and M_0 = 3.3 for A <= 110
52 G4double MaxAverageMultiplicity =
53 G4StatMFParameters::GetMaxAverageMultiplic
54
55
56 // We'll use two kinds of ensembles
57 G4StatMFMicroCanonical * theMicrocanonicalEn
58 G4StatMFMacroCanonical * theMacrocanonicalEn
59
60 //------------------------------------------
61 // Direct simulation part (Microcanonical en
62 //------------------------------------------
63
64 // Microcanonical ensemble initialization
65 theMicrocanonicalEnsemble = new G4StatMFMicr
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 Iterati
76 do {
77
78 G4double theMeanMult = theMicrocanonical
79 if (theMeanMult <= MaxAverageMultiplicit
80 // G4cout << "MICROCANONICAL" << G4endl;
81 // Choose fragments atomic numbers and charg
82 theChannel = theMicrocanonicalEnsemble->Choo
83 _theEnsemble = theMicrocanonicalEnsemble;
84 } else {
85 //------------------------------------------
86 // Non direct simulation part (Macrocanonica
87 //------------------------------------------
88 if (FirstTime) {
89 // Macrocanonical ensemble initialization
90 theMacrocanonicalEnsemble = new G4StatMFMa
91 _theEnsemble = theMacrocanonicalEnsemble;
92 FirstTime = false;
93 }
94 // G4cout << "MACROCANONICAL" << G4endl;
95 // Select calculated fragment total multipli
96 // fragment atomic numbers and fragment char
97 theChannel = theMacrocanonicalEnsemble->Choo
98 }
99
100 ChannelOk = theChannel->CheckFragments()
101 if (!ChannelOk) delete theChannel;
102
103 // Loop checking, 05-Aug-2015, Vladimir
104 } while (!ChannelOk);
105
106
107 if (theChannel->GetMultiplicity() <= 1) {
108 G4FragmentVector * theResult = new G4Fra
109 theResult->push_back(new G4Fragment(theF
110 delete theMicrocanonicalEnsemble;
111 if (theMacrocanonicalEnsemble != 0) dele
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->GetMeanTempera
122
123 if (FindTemperatureOfBreakingChannel(theFr
124
125 // Do not forget to delete this unusable c
126 // otherwise for very proton-reach nuclei
127 // number of iterations. N.B. "theChannel"
128
129 // G4cout << " Iteration # " << Iterations
130
131 delete theChannel;
132
133 // Loop checking, 05-Aug-2015, Vladimir Iv
134 } while (Iterations++ < IterationsLimit );
135
136 // If Iterations >= IterationsLimit means th
137 if (Iterations >= IterationsLimit)
138 throw G4HadronicException(__FILE__, __LINE
139
140 G4FragmentVector * theResult = theChannel->
141 GetFragments(theFragment.GetA_asInt(),theF
142
143 // ~~~~~~ Energy conservation Patch !!!!!!!!
144 // Original nucleus 4-momentum in CM system
145 G4LorentzVector InitialMomentum(theFragment.
146 InitialMomentum.boost(-InitialMomentum.boost
147 G4double ScaleFactor = 0.0;
148 G4double SavedScaleFactor = 0.0;
149 do {
150 G4double FragmentsEnergy = 0.0;
151 for (auto const & ptr : *theResult) {
152 FragmentsEnergy += ptr->GetMomentum().e(
153 }
154 if (0.0 == FragmentsEnergy) { break; }
155 SavedScaleFactor = ScaleFactor;
156 ScaleFactor = InitialMomentum.e()/Fragment
157 G4ThreeVector ScaledMomentum(0.0,0.0,0.0);
158 for (auto const & ptr : *theResult) {
159 ScaledMomentum = ScaleFactor * ptr->GetM
160 G4double Mass = ptr->GetMomentum().mag()
161 G4LorentzVector NewMomentum;
162 NewMomentum.setVect(ScaledMomentum);
163 NewMomentum.setE(std::sqrt(ScaledMomentu
164 ptr->SetMomentum(NewMomentum);
165 }
166 // Loop checking, 05-Aug-2015, Vladimir Iv
167 } while (ScaleFactor > 1.0+1.e-5 && std::abs
168 // ~~~~~~ End of patch !!!!!!!!!!!!!!!!!!!!!
169
170 // Perform Lorentz boost
171 G4FragmentVector::iterator i;
172 for (i = theResult->begin(); i != theResult-
173 G4LorentzVector FourMom = (*i)->GetMomentu
174 FourMom.boost(theFragment.GetMomentum().bo
175 (*i)->SetMomentum(FourMom);
176 (*i)->SetCreatorModelID(_secID);
177 }
178
179 // garbage collection
180 delete theMicrocanonicalEnsemble;
181 if (theMacrocanonicalEnsemble != 0) delete t
182 delete theChannel;
183
184 return theResult;
185 }
186
187
188 G4bool G4StatMF::FindTemperatureOfBreakingChan
189 const G4StatMFChannel * aChannel
190 G4double & Temperature)
191 // This finds temperature of breaking channe
192 {
193 G4int A = theFragment.GetA_asInt();
194 G4int Z = theFragment.GetZ_asInt();
195 G4double U = theFragment.GetExcitationEnergy
196
197 G4double T = std::max(Temperature,0.0012*MeV
198 G4double Ta = T;
199 G4double TotalEnergy = CalcEnergy(A,Z,aChann
200
201 G4double Da = (U - TotalEnergy)/U;
202 G4double Db = 0.0;
203
204 // bracketing the solution
205 if (Da == 0.0) {
206 Temperature = T;
207 return true;
208 } else if (Da < 0.0) {
209 do {
210 T *= 0.5;
211 if (T < 0.001*MeV) return false;
212
213 TotalEnergy = CalcEnergy(A,Z,aChannel,T)
214
215 Db = (U - TotalEnergy)/U;
216 // Loop checking, 05-Aug-2015, Vladimir
217 } while (Db < 0.0);
218
219 } else {
220 do {
221 T *= 1.5;
222
223 TotalEnergy = CalcEnergy(A,Z,aChannel,T)
224
225 Db = (U - TotalEnergy)/U;
226 // Loop checking, 05-Aug-2015, Vladimir
227 } while (Db > 0.0);
228 }
229
230 G4double eps = 1.0e-14 * std::abs(T-Ta);
231 //G4double eps = 1.0e-3 ;
232
233 // Start the bisection method
234 for (G4int j = 0; j < 1000; j++) {
235 G4double Tc = (Ta+T)*0.5;
236 if (std::abs(Ta-Tc) <= eps) {
237 Temperature = Tc;
238 return true;
239 }
240
241 T = Tc;
242 TotalEnergy = CalcEnergy(A,Z,aChannel,T);
243 G4double Dc = (U - TotalEnergy)/U;
244
245 if (Dc == 0.0) {
246 Temperature = Tc;
247 return true;
248 }
249 if (Da*Dc < 0.0) {
250 T = Tc;
251 Db = Dc;
252 } else {
253 Ta = Tc;
254 Da = Dc;
255 }
256 }
257
258 Temperature = (Ta+T)*0.5;
259 return false;
260 }
261
262 G4double G4StatMF::CalcEnergy(G4int A, G4int Z
263 G4double T)
264 {
265 G4double MassExcess0 = G4NucleiProperties::G
266 G4double ChannelEnergy = aChannel->GetFragme
267 return -MassExcess0 + G4StatMFParameters::Ge
268 }
269
270
271
272