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