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