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