| 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 // $Id$
27 // 28 //
28 // Hadronic Process: Nuclear De-excitations 29 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara 30 // by V. Lara
30 // 31 //
31 // Modified: 32 // Modified:
32 // 25.07.08 I.Pshenichnov (in collaboration wi 33 // 25.07.08 I.Pshenichnov (in collaboration with Alexander Botvina and Igor
33 // Mishustin (FIAS, Frankfurt, INR, M 34 // Mishustin (FIAS, Frankfurt, INR, Moscow and Kurchatov Institute,
34 // Moscow, pshenich@fias.uni-frankfur 35 // Moscow, pshenich@fias.uni-frankfurt.de) additional checks in
35 // solver of equation for the chemica 36 // solver of equation for the chemical potential
36 37
37 #include "G4StatMFMacroMultiplicity.hh" 38 #include "G4StatMFMacroMultiplicity.hh"
38 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.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 accessable");
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 accessable");
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 accessable");
59 return true; 59 return true;
60 } 60 }
61 61
>> 62
>> 63
>> 64
62 G4double G4StatMFMacroMultiplicity::CalcChemic 65 G4double G4StatMFMacroMultiplicity::CalcChemicalPotentialMu(void)
63 // Calculate Chemical potential \mu 66 // Calculate Chemical potential \mu
64 // For that is necesary to calculate mean 67 // For that is necesary to calculate mean multiplicities
65 { 68 {
66 G4Pow* g4calc = G4Pow::GetInstance(); << 69 G4double CP = ((3./5.)*elm_coupling/G4StatMFParameters::Getr0())*
67 G4double CP = G4StatMFParameters::GetCoulomb << 70 (1.0-1.0/std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1.0/3.0));
68 71
69 // starting value for chemical potential \mu << 72 // starting value for chemical potential \mu
70 // it is the derivative of F(T,V)-\nu*Z w.r. << 73 // it is the derivative of F(T,V)-\nu*Z w.r.t. Af in Af=5
71 G4double ZA5 = _theClusters->operator[](4)-> << 74 G4double ZA5 = _theClusters->operator[](4)->GetZARatio();
72 G4double ILD5 = _theClusters->operator[](4)- << 75 G4double ILD5 = _theClusters->operator[](4)->GetInvLevelDensity();
73 _ChemPotentialMu = -G4StatMFParameters::GetE << 76 _ChemPotentialMu = -G4StatMFParameters::GetE0()-
74 _MeanTemperature*_MeanTemperature/ILD5 - << 77 _MeanTemperature*_MeanTemperature/ILD5 -
75 _ChemPotentialNu*ZA5 + << 78 _ChemPotentialNu*ZA5 +
76 G4StatMFParameters::GetGamma0()*(1.0-2.0*Z << 79 G4StatMFParameters::GetGamma0()*(1.0-2.0*ZA5)*(1.0-2.0*ZA5) +
77 (2.0/3.0)*G4StatMFParameters::Beta(_MeanTe << 80 (2.0/3.0)*G4StatMFParameters::Beta(_MeanTemperature)/std::pow(5.,1./3.) +
78 (5.0/3.0)*CP*ZA5*ZA5*g4calc->Z23(5) - << 81 (5.0/3.0)*CP*ZA5*ZA5*std::pow(5.,2./3.) -
79 1.5*_MeanTemperature/5.0; << 82 1.5*_MeanTemperature/5.0;
80 83
81 G4double ChemPa = _ChemPotentialMu; << 84
82 if (ChemPa/_MeanTemperature > 10.0) ChemPa = << 85
83 G4double ChemPb = ChemPa - 0.5*std::abs(Chem << 86 G4double ChemPa = _ChemPotentialMu;
84 << 87 if (ChemPa/_MeanTemperature > 10.0) ChemPa = 10.0*_MeanTemperature;
85 G4double fChemPa = this->operator()(ChemPa); << 88 G4double ChemPb = ChemPa - 0.5*std::abs(ChemPa);
86 G4double fChemPb = this->operator()(ChemPb); << 89
87 << 90
88 // Set the precision level for locating the << 91 G4double fChemPa = this->operator()(ChemPa);
89 // If the root is inside this interval, then << 92 G4double fChemPb = this->operator()(ChemPb);
90 const G4double intervalWidth = 1.e-4; << 93
91 << 94
92 // bracketing the solution << 95 // Set the precision level for locating the root.
93 G4int iterations = 0; << 96 // If the root is inside this interval, then it's done!
94 // Loop checking, 05-Aug-2015, Vladimir Ivan << 97 G4double intervalWidth = 1.e-4;
95 while (fChemPa*fChemPb > 0.0 && iterations < << 98
>> 99 // bracketing the solution
>> 100 G4int iterations = 0;
>> 101 while (fChemPa*fChemPb > 0.0 && iterations < 100)
96 { 102 {
97 iterations++; << 103 if (std::abs(fChemPa) <= std::abs(fChemPb))
98 if (std::abs(fChemPa) <= std::abs(fChemP <<
99 { 104 {
100 ChemPa += 0.6*(ChemPa-ChemPb); << 105 ChemPa += 0.6*(ChemPa-ChemPb);
101 fChemPa = this->operator()(ChemPa); << 106 fChemPa = this->operator()(ChemPa);
>> 107 iterations++;
102 } 108 }
103 else << 109 else
104 { 110 {
105 ChemPb += 0.6*(ChemPb-ChemPa); << 111 ChemPb += 0.6*(ChemPb-ChemPa);
106 fChemPb = this->operator()(ChemPb); << 112 fChemPb = this->operator()(ChemPb);
>> 113 iterations++;
107 } 114 }
108 } 115 }
109 116
110 if (fChemPa*fChemPb > 0.0) // the bracketing << 117 if (fChemPa*fChemPb > 0.0) // the bracketing failed, complain
111 { 118 {
112 G4cout <<"G4StatMFMacroMultiplicity:"<<" << 119 G4cerr <<"G4StatMFMacroMultiplicity:"<<" ChemPa="<<ChemPa<<" ChemPb="<<ChemPb<< G4endl;
113 <<" ChemPb="<<ChemPb<< G4endl; << 120 G4cerr <<"G4StatMFMacroMultiplicity:"<<" fChemPa="<<fChemPa<<" fChemPb="<<fChemPb<< G4endl;
114 G4cout <<"G4StatMFMacroMultiplicity:"<<" << 121 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::CalcChemicalPotentialMu: I couldn't bracket the root.");
115 <<" fChemPb="<<fChemPb<< G4endl; <<
116 throw G4HadronicException(__FILE__, __LI <<
117 } 122 }
118 else if (fChemPa*fChemPb < 0.0 && std::abs(C << 123 else if (fChemPa*fChemPb < 0.0 && std::abs(ChemPa-ChemPb) > intervalWidth) // the bracketing was OK, try to locate the root
119 { 124 {
120 G4Solver<G4StatMFMacroMultiplicity> * theS << 125 G4Solver<G4StatMFMacroMultiplicity> * theSolver = new G4Solver<G4StatMFMacroMultiplicity>(100,intervalWidth);
121 new G4Solver<G4StatMFMacroMultiplicity>( <<
122 theSolver->SetIntervalLimits(ChemPa,ChemPb 126 theSolver->SetIntervalLimits(ChemPa,ChemPb);
123 // if (!theSolver->Crenshaw(*this)) 127 // if (!theSolver->Crenshaw(*this))
124 if (!theSolver->Brent(*this)) 128 if (!theSolver->Brent(*this))
125 { 129 {
126 G4cout <<"G4StatMFMacroMultiplicity:"<<" << 130 G4cerr <<"G4StatMFMacroMultiplicity:"<<" ChemPa="<<ChemPa<<" ChemPb="<<ChemPb<< G4endl;
127 <<" ChemPb="<<ChemPb<< G4endl; << 131 G4cerr <<"G4StatMFMacroMultiplicity:"<<" fChemPa="<<fChemPa<<" fChemPb="<<fChemPb<< G4endl;
128 throw G4HadronicException(__FILE__, __LI << 132 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::CalcChemicalPotentialMu: I couldn't find the root.");
129 } 133 }
130 _ChemPotentialMu = theSolver->GetRoot(); 134 _ChemPotentialMu = theSolver->GetRoot();
131 delete theSolver; 135 delete theSolver;
132 } 136 }
133 else // the root is within the interval, whi << 137 else // the root is within the interval, which is shorter then the precision level - all done
134 { 138 {
135 _ChemPotentialMu = ChemPa; 139 _ChemPotentialMu = ChemPa;
136 } 140 }
137 141
138 return _ChemPotentialMu; << 142 return _ChemPotentialMu;
139 } 143 }
140 144
>> 145
>> 146
141 G4double G4StatMFMacroMultiplicity::CalcMeanA( 147 G4double G4StatMFMacroMultiplicity::CalcMeanA(const G4double mu)
142 { 148 {
143 G4double r0 = G4StatMFParameters::Getr0(); << 149 G4double r03 = G4StatMFParameters::Getr0(); r03 *= r03*r03;
144 G4double V0 = (4.0/3.0)*pi*theA*r0*r0*r0; << 150 G4double V0 = (4.0/3.0)*pi*theA*r03;
145 151
146 G4double MeanA = 0.0; 152 G4double MeanA = 0.0;
147 153
148 _MeanMultiplicity = 0.0; 154 _MeanMultiplicity = 0.0;
149 155
>> 156
150 G4int n = 1; 157 G4int n = 1;
151 for (std::vector<G4VStatMFMacroCluster*>::it << 158 for (std::vector<G4VStatMFMacroCluster*>::iterator i = _theClusters->begin();
152 i != _theClusters->end(); ++i) 159 i != _theClusters->end(); ++i)
153 { 160 {
154 G4double multip = (*i)->CalcMeanMultiplic << 161 G4double multip = (*i)->CalcMeanMultiplicity(V0*_Kappa,mu,_ChemPotentialNu,_MeanTemperature);
155 _MeanTemperature); << 162 MeanA += multip*static_cast<G4double>(n++);
156 MeanA += multip*(n++); <<
157 _MeanMultiplicity += multip; 163 _MeanMultiplicity += multip;
158 } 164 }
159 165
160 return MeanA; 166 return MeanA;
161 } 167 }
162 168