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