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