| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 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
41 #include "G4StatMFMacroTemperature.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4Pow.hh"
45
46 G4StatMFMacroTemperature::G4StatMFMacroTempera
47 const G4double ExEnergy, const G4double Free
48 std::vector<G4VStatMFMacroCluster*> * Cluste
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
62 G4StatMFMacroTemperature::~G4StatMFMacroTemper
63 {}
64
65 G4double G4StatMFMacroTemperature::CalcTempera
66 {
67 // Inital guess for the interval of the ense
68 G4double Ta = 0.5;
69 G4double Tb = std::max(std::sqrt(_ExEnergy/(
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
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
133 _MeanTemperature = theSolverBrent->GetRo
134 FunctionValureAtRoot = this->operator()
135 delete theSolverBrent;
136 }
137 if (std::abs(FunctionValureAtRoot) > 5.e-2
138 G4cout << "Brent method failed; function
139 << " solution? = " << _MeanTemperature
140 throw G4HadronicException(__FILE__, __LI
141 }
142 }
143 //G4cout << "G4StatMFMacroTemperature::CalcT
144 //<< FunctionValureAtRoot
145 // << " T(MeV)= " << _MeanTemperature <<
146 return _MeanTemperature;
147 }
148
149 G4double G4StatMFMacroTemperature::FragsExcitE
150 // Calculates excitation energy per nucleon an
151 // multiplicity and entropy
152 {
153 // Model Parameters
154 G4Pow* g4calc = G4Pow::GetInstance();
155 G4double R0 = G4StatMFParameters::Getr0()*g4
156 G4double R = R0*g4calc->A13(1.0+G4StatMFPara
157 G4double FreeVol = _Kappa*(4.*pi/3.)*R0*R0*R
158
159 // Calculate Chemical potentials
160 CalcChemicalPotentialNu(T);
161
162
163 // Average total fragment energy
164 G4double AverageEnergy = 0.0;
165 std::vector<G4VStatMFMacroCluster*>::iterato
166 for (i = _theClusters->begin(); i != _theCl
167 {
168 AverageEnergy += (*i)->GetMeanMultiplici
169 }
170
171 // Add Coulomb energy
172 AverageEnergy += 0.6*elm_coupling*theZ*theZ/
173
174 // Calculate mean entropy
175 _MeanEntropy = 0.0;
176 for (i = _theClusters->begin(); i != _theClu
177 {
178 _MeanEntropy += (*i)->CalcEntropy(T,Free
179 }
180
181 // Excitation energy per nucleon
182 return AverageEnergy - _FreeInternalE0;
183 }
184
185 void G4StatMFMacroTemperature::CalcChemicalPot
186 // Calculates the chemical potential \nu
187 {
188 G4StatMFMacroChemicalPotential * theChemPot
189 G4StatMFMacroChemicalPotential(theA,theZ,_
190
191 _ChemPotentialNu = theChemPot->CalcChemicalP
192 _ChemPotentialMu = theChemPot->GetChemicalPo
193 _MeanMultiplicity = theChemPot->GetMeanMulti
194 delete theChemPot;
195
196 return;
197 }
198
199
200