| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // 23 //
>> 24 // $Id: G4StatMFMacroTemperature.cc,v 1.2 2003/11/03 17:53:05 hpw Exp $
>> 25 // GEANT4 tag $Name: geant4-06-00-patch-01 $
27 // 26 //
28 // Hadronic Process: Nuclear De-excitations 27 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara 28 // by V. Lara
30 // << 29
31 // Modified: <<
32 // 25.07.08 I.Pshenichnov (in collaboration wi <<
33 // Mishustin (FIAS, Frankfurt, INR, M <<
34 // Moscow, pshenich@fias.uni-frankfur <<
35 // original MF model <<
36 // 16.04.10 V.Ivanchenko improved logic of sol <<
37 // to protect code from rare unwanted <<
38 // and destructor to source <<
39 // 28.10.10 V.Ivanchenko defined members in co <<
40 30
41 #include "G4StatMFMacroTemperature.hh" 31 #include "G4StatMFMacroTemperature.hh"
42 #include "G4PhysicalConstants.hh" << 32
43 #include "G4SystemOfUnits.hh" << 33 // operators definitions
44 #include "G4Pow.hh" << 34 G4StatMFMacroTemperature &
45 << 35 G4StatMFMacroTemperature::operator=(const G4StatMFMacroTemperature & )
46 G4StatMFMacroTemperature::G4StatMFMacroTempera << 36 {
47 const G4double ExEnergy, const G4double Free << 37 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::operator= meant to not be accessable");
48 std::vector<G4VStatMFMacroCluster*> * Cluste << 38 return *this;
49 theA(anA), << 39 }
50 theZ(aZ), << 40
51 _ExEnergy(ExEnergy), << 41 G4bool G4StatMFMacroTemperature::operator==(const G4StatMFMacroTemperature & ) const
52 _FreeInternalE0(FreeE0), << 42 {
53 _Kappa(kappa), << 43 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::operator== meant to not be accessable");
54 _MeanMultiplicity(0.0), << 44 return false;
55 _MeanTemperature(0.0), << 45 }
56 _ChemPotentialMu(0.0), << 46
57 _ChemPotentialNu(0.0), << 47
58 _MeanEntropy(0.0), << 48 G4bool G4StatMFMacroTemperature::operator!=(const G4StatMFMacroTemperature & ) const
59 _theClusters(ClusterVector) << 49 {
60 {} << 50 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::operator!= meant to not be accessable");
61 << 51 return true;
62 G4StatMFMacroTemperature::~G4StatMFMacroTemper << 52 }
63 {} << 53
>> 54
>> 55
64 56
65 G4double G4StatMFMacroTemperature::CalcTempera 57 G4double G4StatMFMacroTemperature::CalcTemperature(void)
>> 58 // Calculate Chemical potential \nu
66 { 59 {
67 // Inital guess for the interval of the ense << 60 // Temperature
68 G4double Ta = 0.5; << 61 G4double Ta = 0.00012;
69 G4double Tb = std::max(std::sqrt(_ExEnergy/( << 62 G4double Tb = std::max(sqrt(_ExEnergy/(theA*0.12)),0.01*MeV);
70 63
71 G4double fTa = this->operator()(Ta); << 64 G4double fTa = this->operator()(Ta);
72 G4double fTb = this->operator()(Tb); << 65 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 <<
79 // the case, we decrease Ta. But carefully, <<
80 // fTa growes very fast when Ta is near 0 an <<
81 // an overflow. <<
82 <<
83 G4int iterations = 0; <<
84 // Loop checking, 05-Aug-2015, Vladimir Ivan <<
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 <<
90 iterations = 0; <<
91 // Loop checking, 05-Aug-2015, Vladimir Ivan <<
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 <<
99 G4cerr <<"G4StatMFMacroTemperature:"<<" fT <<
100 throw G4HadronicException(__FILE__, __LINE <<
101 } <<
102 <<
103 G4Solver<G4StatMFMacroTemperature> * theSolv <<
104 new G4Solver<G4StatMFMacroTemperature>(100 <<
105 theSolver->SetIntervalLimits(Ta,Tb); <<
106 if (!theSolver->Crenshaw(*this)){ <<
107 G4cout <<"G4StatMFMacroTemperature, Crensh <<
108 <<Ta<<" Tb="<<Tb<< G4endl; <<
109 G4cout <<"G4StatMFMacroTemperature, Crensh <<
110 <<fTa<<" fTb="<<fTb<< G4endl; <<
111 } <<
112 _MeanTemperature = theSolver->GetRoot(); <<
113 G4double FunctionValureAtRoot = this->opera <<
114 delete theSolver; <<
115 <<
116 // Verify if the root is found and it is ind <<
117 // say, between 1 and 50 MeV, otherwise try <<
118 if (std::fabs(FunctionValureAtRoot) > 5.e-2) <<
119 if (_MeanTemperature < 1. || _MeanTemperat <<
120 G4cout << "Crenshaw method failed; funct <<
121 << " solution? = " << _MeanTemperature <<
122 G4Solver<G4StatMFMacroTemperature> * the <<
123 new G4Solver<G4StatMFMacroTemperature>(200,1 <<
124 theSolverBrent->SetIntervalLimits(Ta,Tb) <<
125 if (!theSolverBrent->Brent(*this)){ <<
126 G4cout <<"G4StatMFMacroTemperature, Brent me <<
127 <<" Ta="<<Ta<<" Tb="<<Tb<< G4endl; <<
128 G4cout <<"G4StatMFMacroTemperature, Brent me <<
129 <<" fTa="<<fTa<<" fTb="<<fTb<< G4endl <<
130 throw G4HadronicException(__FILE__, __LINE__ <<
131 } <<
132 66
133 _MeanTemperature = theSolverBrent->GetRo << 67 // Bracketing the solution
134 FunctionValureAtRoot = this->operator() << 68 // T should be greater than 0.
135 delete theSolverBrent; << 69 // The interval is [Ta,Tb]
>> 70 // We start with a low value for Ta = 0.0012 K
>> 71 // it should be enough to have fTa > 0 If it isn't
>> 72 // the case, we decrease Ta. But carefully, because
>> 73 // fTa growes very fast when Ta is near 0 and we could have
>> 74 // an overflow.
>> 75
>> 76 G4int iterations = 0;
>> 77 while (fTa < 0.0 && iterations++ < 10) {
>> 78 Ta -= 0.5*Ta;
>> 79 fTa = this->operator()(Ta);
136 } 80 }
137 if (std::abs(FunctionValureAtRoot) > 5.e-2 << 81 // Usually, fTb will be less than 0, but if it is not the case:
138 G4cout << "Brent method failed; function << 82 iterations = 0;
139 << " solution? = " << _MeanTemperature << 83 while (fTa*fTb > 0.0 && iterations++ < 10) {
140 throw G4HadronicException(__FILE__, __LI << 84 Tb += 1.5*abs(Tb-Ta);
>> 85 fTb = this->operator()(Tb);
141 } 86 }
142 } << 87
143 //G4cout << "G4StatMFMacroTemperature::CalcT << 88 if (fTa*fTb > 0.0) {
144 //<< FunctionValureAtRoot << 89 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't bracket the solution.");
145 // << " T(MeV)= " << _MeanTemperature << << 90 }
146 return _MeanTemperature; << 91
>> 92 G4Solver<G4StatMFMacroTemperature> * theSolver = new G4Solver<G4StatMFMacroTemperature>(100,1.e-4);
>> 93 theSolver->SetIntervalLimits(Ta,Tb);
>> 94 // if (!theSolver->Crenshaw(*this))
>> 95 if (!theSolver->Brent(*this))
>> 96 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root.");
>> 97 _MeanTemperature = theSolver->GetRoot();
>> 98 delete theSolver;
>> 99 return _MeanTemperature;
147 } 100 }
148 101
>> 102
>> 103
149 G4double G4StatMFMacroTemperature::FragsExcitE 104 G4double G4StatMFMacroTemperature::FragsExcitEnergy(const G4double T)
150 // Calculates excitation energy per nucleon an << 105 // Calculates excitation energy per nucleon and summed fragment multiplicity and entropy
151 // multiplicity and entropy <<
152 { 106 {
153 // Model Parameters << 107
154 G4Pow* g4calc = G4Pow::GetInstance(); << 108 // Model Parameters
155 G4double R0 = G4StatMFParameters::Getr0()*g4 << 109 G4double R0 = G4StatMFParameters::Getr0()*pow(theA,1./3.);
156 G4double R = R0*g4calc->A13(1.0+G4StatMFPara << 110 G4double R = R0*pow(1.0+G4StatMFParameters::GetKappaCoulomb(), 1./3.);
157 G4double FreeVol = _Kappa*(4.*pi/3.)*R0*R0*R << 111 G4double FreeVol = _Kappa*(4.*pi/3.)*R0*R0*R0;
>> 112
158 113
159 // Calculate Chemical potentials << 114 // Calculate Chemical potentials
160 CalcChemicalPotentialNu(T); << 115 CalcChemicalPotentialNu(T);
161 116
162 117
163 // Average total fragment energy << 118 // Average total fragment energy
164 G4double AverageEnergy = 0.0; << 119 G4double AverageEnergy = 0.0;
165 std::vector<G4VStatMFMacroCluster*>::iterato << 120 std::vector<G4VStatMFMacroCluster*>::iterator i;
166 for (i = _theClusters->begin(); i != _theCl << 121 for (i = _theClusters->begin(); i != _theClusters->end(); ++i)
167 { << 122 {
168 AverageEnergy += (*i)->GetMeanMultiplici << 123 AverageEnergy += (*i)->GetMeanMultiplicity() * (*i)->CalcEnergy(T);
169 } << 124 }
170 125
171 // Add Coulomb energy << 126 // Add Coulomb energy
172 AverageEnergy += 0.6*elm_coupling*theZ*theZ/ << 127 AverageEnergy += (3./5.)*elm_coupling*theZ*theZ/R;
173 128
174 // Calculate mean entropy << 129 // Calculate mean entropy
175 _MeanEntropy = 0.0; << 130 _MeanEntropy = 0.0;
176 for (i = _theClusters->begin(); i != _theClu << 131 for (i = _theClusters->begin(); i != _theClusters->end(); ++i)
177 { << 132 {
178 _MeanEntropy += (*i)->CalcEntropy(T,Free << 133 _MeanEntropy += (*i)->CalcEntropy(T,FreeVol);
179 } << 134 }
>> 135
>> 136 // Excitation energy per nucleon
>> 137 G4double FragsExcitEnergy = AverageEnergy - _FreeInternalE0;
>> 138
>> 139 return FragsExcitEnergy;
180 140
181 // Excitation energy per nucleon <<
182 return AverageEnergy - _FreeInternalE0; <<
183 } 141 }
184 142
>> 143
185 void G4StatMFMacroTemperature::CalcChemicalPot 144 void G4StatMFMacroTemperature::CalcChemicalPotentialNu(const G4double T)
186 // Calculates the chemical potential \nu << 145 // Calculates the chemical potential \nu
>> 146
187 { 147 {
188 G4StatMFMacroChemicalPotential * theChemPot << 148 G4StatMFMacroChemicalPotential * theChemPot = new
189 G4StatMFMacroChemicalPotential(theA,theZ,_ << 149 G4StatMFMacroChemicalPotential(theA,theZ,_Kappa,T,_theClusters);
>> 150
190 151
191 _ChemPotentialNu = theChemPot->CalcChemicalP << 152 _ChemPotentialNu = theChemPot->CalcChemicalPotentialNu();
192 _ChemPotentialMu = theChemPot->GetChemicalPo << 153 _ChemPotentialMu = theChemPot->GetChemicalPotentialMu();
193 _MeanMultiplicity = theChemPot->GetMeanMulti << 154 _MeanMultiplicity = theChemPot->GetMeanMultiplicity();
194 delete theChemPot; << 155
>> 156 delete theChemPot;
195 157
196 return; << 158 return;
>> 159
197 } 160 }
198 161
199 162
200 163