Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4StatMFMacroTemperature.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: G4StatMFMacroTemperature.cc 100379 2016-10-19 15:05:35Z gcosmo $
28 //
29 // Hadronic Process: Nuclear De-excitations
30 // by V. Lara
31 //
32 // Modified:
33 // 25.07.08 I.Pshenichnov (in collaboration with Alexander Botvina and Igor
34 // Mishustin (FIAS, Frankfurt, INR, Moscow and Kurchatov Institute,
35 // Moscow, pshenich@fias.uni-frankfurt.de) make algorithm closer to
36 // original MF model
37 // 16.04.10 V.Ivanchenko improved logic of solving equation for tempetature
38 // to protect code from rare unwanted exception; moved constructor
39 // and destructor to source
40 // 28.10.10 V.Ivanchenko defined members in constructor and cleaned up
41 
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4Pow.hh"
46 
48  const G4double ExEnergy, const G4double FreeE0, const G4double kappa,
49  std::vector<G4VStatMFMacroCluster*> * ClusterVector) :
50  theA(anA),
51  theZ(aZ),
52  _ExEnergy(ExEnergy),
53  _FreeInternalE0(FreeE0),
54  _Kappa(kappa),
55  _MeanMultiplicity(0.0),
56  _MeanTemperature(0.0),
57  _ChemPotentialMu(0.0),
58  _ChemPotentialNu(0.0),
59  _MeanEntropy(0.0),
60  _theClusters(ClusterVector)
61 {}
62 
64 {}
65 
67 {
68  // Inital guess for the interval of the ensemble temperature values
69  G4double Ta = 0.5;
70  G4double Tb = std::max(std::sqrt(_ExEnergy/(theA*0.12)),0.01*MeV);
71 
72  G4double fTa = this->operator()(Ta);
73  G4double fTb = this->operator()(Tb);
74 
75  // Bracketing the solution
76  // T should be greater than 0.
77  // The interval is [Ta,Tb]
78  // We start with a value for Ta = 0.5 MeV
79  // it should be enough to have fTa > 0 If it isn't
80  // the case, we decrease Ta. But carefully, because
81  // fTa growes very fast when Ta is near 0 and we could have
82  // an overflow.
83 
84  G4int iterations = 0;
85  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
86  while (fTa < 0.0 && ++iterations < 10) {
87  Ta -= 0.5*Ta;
88  fTa = this->operator()(Ta);
89  }
90  // Usually, fTb will be less than 0, but if it is not the case:
91  iterations = 0;
92  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
93  while (fTa*fTb > 0.0 && iterations++ < 10) {
94  Tb += 2.*std::fabs(Tb-Ta);
95  fTb = this->operator()(Tb);
96  }
97 
98  if (fTa*fTb > 0.0) {
99  G4cerr <<"G4StatMFMacroTemperature:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
100  G4cerr <<"G4StatMFMacroTemperature:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl;
101  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't bracket the solution.");
102  }
103 
105  new G4Solver<G4StatMFMacroTemperature>(100,1.e-4);
106  theSolver->SetIntervalLimits(Ta,Tb);
107  if (!theSolver->Crenshaw(*this)){
108  G4cout <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" Ta="
109  <<Ta<<" Tb="<<Tb<< G4endl;
110  G4cout <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" fTa="
111  <<fTa<<" fTb="<<fTb<< G4endl;
112  }
113  _MeanTemperature = theSolver->GetRoot();
114  G4double FunctionValureAtRoot = this->operator()(_MeanTemperature);
115  delete theSolver;
116 
117  // Verify if the root is found and it is indeed within the physical domain,
118  // say, between 1 and 50 MeV, otherwise try Brent method:
119  if (std::fabs(FunctionValureAtRoot) > 5.e-2) {
120  if (_MeanTemperature < 1. || _MeanTemperature > 50.) {
121  G4cout << "Crenshaw method failed; function = " << FunctionValureAtRoot
122  << " solution? = " << _MeanTemperature << " MeV " << G4endl;
123  G4Solver<G4StatMFMacroTemperature> * theSolverBrent =
124  new G4Solver<G4StatMFMacroTemperature>(200,1.e-3);
125  theSolverBrent->SetIntervalLimits(Ta,Tb);
126  if (!theSolverBrent->Brent(*this)){
127  G4cout <<"G4StatMFMacroTemperature, Brent method failed:"
128  <<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
129  G4cout <<"G4StatMFMacroTemperature, Brent method failed:"
130  <<" fTa="<<fTa<<" fTb="<<fTb<< G4endl;
131  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root with any method.");
132  }
133 
134  _MeanTemperature = theSolverBrent->GetRoot();
135  FunctionValureAtRoot = this->operator()(_MeanTemperature);
136  delete theSolverBrent;
137  }
138  if (std::abs(FunctionValureAtRoot) > 5.e-2) {
139  G4cout << "Brent method failed; function = " << FunctionValureAtRoot
140  << " solution? = " << _MeanTemperature << " MeV " << G4endl;
141  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root with any method.");
142  }
143  }
144  //G4cout << "G4StatMFMacroTemperature::CalcTemperature: function = "
145  //<< FunctionValureAtRoot
146  // << " T(MeV)= " << _MeanTemperature << G4endl;
147  return _MeanTemperature;
148 }
149 
150 G4double G4StatMFMacroTemperature::FragsExcitEnergy(const G4double T)
151 // Calculates excitation energy per nucleon and summed fragment
152 // multiplicity and entropy
153 {
154  // Model Parameters
155  G4Pow* g4calc = G4Pow::GetInstance();
156  G4double R0 = G4StatMFParameters::Getr0()*g4calc->Z13(theA);
157  G4double R = R0*g4calc->A13(1.0+G4StatMFParameters::GetKappaCoulomb());
158  G4double FreeVol = _Kappa*(4.*pi/3.)*R0*R0*R0;
159 
160  // Calculate Chemical potentials
161  CalcChemicalPotentialNu(T);
162 
163 
164  // Average total fragment energy
165  G4double AverageEnergy = 0.0;
166  std::vector<G4VStatMFMacroCluster*>::iterator i;
167  for (i = _theClusters->begin(); i != _theClusters->end(); ++i)
168  {
169  AverageEnergy += (*i)->GetMeanMultiplicity() * (*i)->CalcEnergy(T);
170  }
171 
172  // Add Coulomb energy
173  AverageEnergy += 0.6*elm_coupling*theZ*theZ/R;
174 
175  // Calculate mean entropy
176  _MeanEntropy = 0.0;
177  for (i = _theClusters->begin(); i != _theClusters->end(); ++i)
178  {
179  _MeanEntropy += (*i)->CalcEntropy(T,FreeVol);
180  }
181 
182  // Excitation energy per nucleon
183  return AverageEnergy - _FreeInternalE0;
184 }
185 
186 void G4StatMFMacroTemperature::CalcChemicalPotentialNu(const G4double T)
187 // Calculates the chemical potential \nu
188 {
189  G4StatMFMacroChemicalPotential * theChemPot = new
190  G4StatMFMacroChemicalPotential(theA,theZ,_Kappa,T,_theClusters);
191 
192  _ChemPotentialNu = theChemPot->CalcChemicalPotentialNu();
193  _ChemPotentialMu = theChemPot->GetChemicalPotentialMu();
194  _MeanMultiplicity = theChemPot->GetMeanMultiplicity();
195  delete theChemPot;
196 
197  return;
198 }
199 
200 
G4StatMFMacroTemperature(const G4double anA, const G4double aZ, const G4double ExEnergy, const G4double FreeE0, const G4double kappa, std::vector< G4VStatMFMacroCluster * > *ClusterVector)
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
static G4double GetKappaCoulomb()
Definition: G4Pow.hh:56
tuple elm_coupling
Definition: hepunit.py:286
void SetIntervalLimits(const G4double Limit1, const G4double Limit2)
int G4int
Definition: G4Types.hh:78
G4bool Brent(Function &theFunction)
static G4double Getr0()
G4GLOB_DLL std::ostream G4cout
G4double Z13(G4int Z) const
Definition: G4Pow.hh:127
G4double GetRoot(void) const
Definition: G4Solver.hh:77
G4double A13(G4double A) const
Definition: G4Pow.hh:132
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double operator()(const G4double T)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
static constexpr double pi
Definition: G4SIunits.hh:75
G4bool Crenshaw(Function &theFunction)
double G4double
Definition: G4Types.hh:76
G4GLOB_DLL std::ostream G4cerr