| Geant4 Cross Reference |
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) additional checks in
35 // solver of equation for the chemical potential
36
37 #include "G4StatMFMacroMultiplicity.hh"
38 #include "G4PhysicalConstants.hh"
39 #include "G4Pow.hh"
40
41 // operators definitions
42 G4StatMFMacroMultiplicity &
43 G4StatMFMacroMultiplicity::operator=(const G4StatMFMacroMultiplicity & )
44 {
45 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::operator= meant to not be accessible");
46 return *this;
47 }
48
49 G4bool G4StatMFMacroMultiplicity::operator==(const G4StatMFMacroMultiplicity & ) const
50 {
51 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::operator== meant to not be accessible");
52 return false;
53 }
54
55
56 G4bool G4StatMFMacroMultiplicity::operator!=(const G4StatMFMacroMultiplicity & ) const
57 {
58 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::operator!= meant to not be accessible");
59 return true;
60 }
61
62 G4double G4StatMFMacroMultiplicity::CalcChemicalPotentialMu(void)
63 // Calculate Chemical potential \mu
64 // For that is necesary to calculate mean multiplicities
65 {
66 G4Pow* g4calc = G4Pow::GetInstance();
67 G4double CP = G4StatMFParameters::GetCoulomb();
68
69 // starting value for chemical potential \mu
70 // it is the derivative of F(T,V)-\nu*Z w.r.t. Af in Af=5
71 G4double ZA5 = _theClusters->operator[](4)->GetZARatio();
72 G4double ILD5 = _theClusters->operator[](4)->GetInvLevelDensity();
73 _ChemPotentialMu = -G4StatMFParameters::GetE0()-
74 _MeanTemperature*_MeanTemperature/ILD5 -
75 _ChemPotentialNu*ZA5 +
76 G4StatMFParameters::GetGamma0()*(1.0-2.0*ZA5)*(1.0-2.0*ZA5) +
77 (2.0/3.0)*G4StatMFParameters::Beta(_MeanTemperature)/g4calc->Z13(5) +
78 (5.0/3.0)*CP*ZA5*ZA5*g4calc->Z23(5) -
79 1.5*_MeanTemperature/5.0;
80
81 G4double ChemPa = _ChemPotentialMu;
82 if (ChemPa/_MeanTemperature > 10.0) ChemPa = 10.0*_MeanTemperature;
83 G4double ChemPb = ChemPa - 0.5*std::abs(ChemPa);
84
85 G4double fChemPa = this->operator()(ChemPa);
86 G4double fChemPb = this->operator()(ChemPb);
87
88 // Set the precision level for locating the root.
89 // If the root is inside this interval, then it's done!
90 const G4double intervalWidth = 1.e-4;
91
92 // bracketing the solution
93 G4int iterations = 0;
94 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
95 while (fChemPa*fChemPb > 0.0 && iterations < 100)
96 {
97 iterations++;
98 if (std::abs(fChemPa) <= std::abs(fChemPb))
99 {
100 ChemPa += 0.6*(ChemPa-ChemPb);
101 fChemPa = this->operator()(ChemPa);
102 }
103 else
104 {
105 ChemPb += 0.6*(ChemPb-ChemPa);
106 fChemPb = this->operator()(ChemPb);
107 }
108 }
109
110 if (fChemPa*fChemPb > 0.0) // the bracketing failed, complain
111 {
112 G4cout <<"G4StatMFMacroMultiplicity:"<<" ChemPa="<<ChemPa
113 <<" ChemPb="<<ChemPb<< G4endl;
114 G4cout <<"G4StatMFMacroMultiplicity:"<<" fChemPa="<<fChemPa
115 <<" fChemPb="<<fChemPb<< G4endl;
116 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::CalcChemicalPotentialMu: I couldn't bracket the root.");
117 }
118 else if (fChemPa*fChemPb < 0.0 && std::abs(ChemPa-ChemPb) > intervalWidth)
119 {
120 G4Solver<G4StatMFMacroMultiplicity> * theSolver =
121 new G4Solver<G4StatMFMacroMultiplicity>(100,intervalWidth);
122 theSolver->SetIntervalLimits(ChemPa,ChemPb);
123 // if (!theSolver->Crenshaw(*this))
124 if (!theSolver->Brent(*this))
125 {
126 G4cout <<"G4StatMFMacroMultiplicity:"<<" ChemPa="<<ChemPa
127 <<" ChemPb="<<ChemPb<< G4endl;
128 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::CalcChemicalPotentialMu: I couldn't find the root.");
129 }
130 _ChemPotentialMu = theSolver->GetRoot();
131 delete theSolver;
132 }
133 else // the root is within the interval, which is shorter then the precision level - all done
134 {
135 _ChemPotentialMu = ChemPa;
136 }
137
138 return _ChemPotentialMu;
139 }
140
141 G4double G4StatMFMacroMultiplicity::CalcMeanA(const G4double mu)
142 {
143 G4double r0 = G4StatMFParameters::Getr0();
144 G4double V0 = (4.0/3.0)*pi*theA*r0*r0*r0;
145
146 G4double MeanA = 0.0;
147
148 _MeanMultiplicity = 0.0;
149
150 G4int n = 1;
151 for (std::vector<G4VStatMFMacroCluster*>::iterator i = _theClusters->begin();
152 i != _theClusters->end(); ++i)
153 {
154 G4double multip = (*i)->CalcMeanMultiplicity(V0*_Kappa,mu,_ChemPotentialNu,
155 _MeanTemperature);
156 MeanA += multip*(n++);
157 _MeanMultiplicity += multip;
158 }
159
160 return MeanA;
161 }
162