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