| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // 23 //
>> 24 // $Id: G4StatMFMacroMultiplicity.cc,v 1.4 2005/06/04 13:27:49 jwellisc Exp $
>> 25 // GEANT4 tag $Name: geant4-07-01-patch-01 $
27 // 26 //
28 // Hadronic Process: Nuclear De-excitations 27 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara 28 // by V. Lara
30 // << 29
31 // Modified: <<
32 // 25.07.08 I.Pshenichnov (in collaboration wi <<
33 // Mishustin (FIAS, Frankfurt, INR, M <<
34 // Moscow, pshenich@fias.uni-frankfur <<
35 // solver of equation for the chemica <<
36 30
37 #include "G4StatMFMacroMultiplicity.hh" 31 #include "G4StatMFMacroMultiplicity.hh"
38 #include "G4PhysicalConstants.hh" <<
39 #include "G4Pow.hh" <<
40 32
41 // operators definitions 33 // operators definitions
42 G4StatMFMacroMultiplicity & 34 G4StatMFMacroMultiplicity &
43 G4StatMFMacroMultiplicity::operator=(const G4S 35 G4StatMFMacroMultiplicity::operator=(const G4StatMFMacroMultiplicity & )
44 { 36 {
45 throw G4HadronicException(__FILE__, __LINE << 37 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::operator= meant to not be accessable");
46 return *this; 38 return *this;
47 } 39 }
48 40
49 G4bool G4StatMFMacroMultiplicity::operator==(c 41 G4bool G4StatMFMacroMultiplicity::operator==(const G4StatMFMacroMultiplicity & ) const
50 { 42 {
51 throw G4HadronicException(__FILE__, __LINE << 43 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::operator== meant to not be accessable");
52 return false; 44 return false;
53 } 45 }
54 46
55 47
56 G4bool G4StatMFMacroMultiplicity::operator!=(c 48 G4bool G4StatMFMacroMultiplicity::operator!=(const G4StatMFMacroMultiplicity & ) const
57 { 49 {
58 throw G4HadronicException(__FILE__, __LINE << 50 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::operator!= meant to not be accessable");
59 return true; 51 return true;
60 } 52 }
61 53
>> 54
>> 55
>> 56
62 G4double G4StatMFMacroMultiplicity::CalcChemic 57 G4double G4StatMFMacroMultiplicity::CalcChemicalPotentialMu(void)
63 // Calculate Chemical potential \mu 58 // Calculate Chemical potential \mu
64 // For that is necesary to calculate mean 59 // For that is necesary to calculate mean multiplicities
65 { 60 {
66 G4Pow* g4calc = G4Pow::GetInstance(); << 61 G4double CP = ((3./5.)*elm_coupling/G4StatMFParameters::Getr0())*
67 G4double CP = G4StatMFParameters::GetCoulomb << 62 (1.0-1.0/std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1.0/3.0));
68 63
69 // starting value for chemical potential \mu << 64 // starting value for chemical potential \mu
70 // it is the derivative of F(T,V)-\nu*Z w.r. << 65 // it is the derivative of F(T,V)-\nu*Z w.r.t. Af in Af=5
71 G4double ZA5 = _theClusters->operator[](4)-> << 66 G4double ZA5 = _theClusters->operator[](4)->GetZARatio();
72 G4double ILD5 = _theClusters->operator[](4)- << 67 G4double ILD5 = _theClusters->operator[](4)->GetInvLevelDensity();
73 _ChemPotentialMu = -G4StatMFParameters::GetE << 68 _ChemPotentialMu = -G4StatMFParameters::GetE0()-
74 _MeanTemperature*_MeanTemperature/ILD5 - << 69 _MeanTemperature*_MeanTemperature/ILD5 -
75 _ChemPotentialNu*ZA5 + << 70 _ChemPotentialNu*ZA5 +
76 G4StatMFParameters::GetGamma0()*(1.0-2.0*Z << 71 G4StatMFParameters::GetGamma0()*(1.0-2.0*ZA5)*(1.0-2.0*ZA5) +
77 (2.0/3.0)*G4StatMFParameters::Beta(_MeanTe << 72 (2.0/3.0)*G4StatMFParameters::Beta(_MeanTemperature)/std::pow(5.,1./3.) +
78 (5.0/3.0)*CP*ZA5*ZA5*g4calc->Z23(5) - << 73 (5.0/3.0)*CP*ZA5*ZA5*std::pow(5.,2./3.) -
79 1.5*_MeanTemperature/5.0; << 74 1.5*_MeanTemperature/5.0;
80 75
81 G4double ChemPa = _ChemPotentialMu; << 76
82 if (ChemPa/_MeanTemperature > 10.0) ChemPa = << 77
83 G4double ChemPb = ChemPa - 0.5*std::abs(Chem << 78 G4double ChemPa = _ChemPotentialMu;
84 << 79 if (ChemPa/_MeanTemperature > 10.0) ChemPa = 10.0*_MeanTemperature;
85 G4double fChemPa = this->operator()(ChemPa); << 80 G4double ChemPb = ChemPa - 0.5*std::abs(ChemPa);
86 G4double fChemPb = this->operator()(ChemPb); << 81
87 << 82
88 // Set the precision level for locating the << 83 G4double fChemPa = this->operator()(ChemPa);
89 // If the root is inside this interval, then << 84 G4double fChemPb = this->operator()(ChemPb);
90 const G4double intervalWidth = 1.e-4; << 85
91 << 86 // bracketing the solution
92 // bracketing the solution << 87 G4int iterations = 0;
93 G4int iterations = 0; << 88 while (fChemPa*fChemPb > 0.0 && iterations < 10)
94 // Loop checking, 05-Aug-2015, Vladimir Ivan <<
95 while (fChemPa*fChemPb > 0.0 && iterations < <<
96 { 89 {
97 iterations++; << 90 if (std::abs(fChemPa) <= std::abs(fChemPb))
98 if (std::abs(fChemPa) <= std::abs(fChemP <<
99 { 91 {
100 ChemPa += 0.6*(ChemPa-ChemPb); << 92 ChemPa += 0.6*(ChemPa-ChemPb);
101 fChemPa = this->operator()(ChemPa); << 93 fChemPa = this->operator()(ChemPa);
102 } 94 }
103 else << 95 else
104 { 96 {
105 ChemPb += 0.6*(ChemPb-ChemPa); << 97 ChemPb += 0.6*(ChemPb-ChemPa);
106 fChemPb = this->operator()(ChemPb); << 98 fChemPb = this->operator()(ChemPb);
107 } 99 }
108 } 100 }
109 << 101 if (fChemPa*fChemPb > 0.0)
110 if (fChemPa*fChemPb > 0.0) // the bracketing <<
111 { 102 {
112 G4cout <<"G4StatMFMacroMultiplicity:"<<" << 103 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::CalcChemicalPotentialMu: I couldn't bracket the root.");
113 <<" ChemPb="<<ChemPb<< G4endl; <<
114 G4cout <<"G4StatMFMacroMultiplicity:"<<" <<
115 <<" fChemPb="<<fChemPb<< G4endl; <<
116 throw G4HadronicException(__FILE__, __LI <<
117 } 104 }
118 else if (fChemPa*fChemPb < 0.0 && std::abs(C << 105
119 { << 106
120 G4Solver<G4StatMFMacroMultiplicity> * theS << 107 G4Solver<G4StatMFMacroMultiplicity> * theSolver = new G4Solver<G4StatMFMacroMultiplicity>(100,1.e-4);
121 new G4Solver<G4StatMFMacroMultiplicity>( <<
122 theSolver->SetIntervalLimits(ChemPa,ChemPb 108 theSolver->SetIntervalLimits(ChemPa,ChemPb);
123 // if (!theSolver->Crenshaw(*this)) 109 // if (!theSolver->Crenshaw(*this))
124 if (!theSolver->Brent(*this)) 110 if (!theSolver->Brent(*this))
125 { 111 {
126 G4cout <<"G4StatMFMacroMultiplicity:"<<" << 112 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::CalcChemicalPotentialMu: I couldn't find the root.");
127 <<" ChemPb="<<ChemPb<< G4endl; <<
128 throw G4HadronicException(__FILE__, __LI <<
129 } 113 }
130 _ChemPotentialMu = theSolver->GetRoot(); 114 _ChemPotentialMu = theSolver->GetRoot();
131 delete theSolver; 115 delete theSolver;
132 } << 116 return _ChemPotentialMu;
133 else // the root is within the interval, whi <<
134 { <<
135 _ChemPotentialMu = ChemPa; <<
136 } <<
137 <<
138 return _ChemPotentialMu; <<
139 } 117 }
140 118
>> 119
>> 120
141 G4double G4StatMFMacroMultiplicity::CalcMeanA( 121 G4double G4StatMFMacroMultiplicity::CalcMeanA(const G4double mu)
142 { 122 {
143 G4double r0 = G4StatMFParameters::Getr0(); << 123 G4double r03 = G4StatMFParameters::Getr0(); r03 *= r03*r03;
144 G4double V0 = (4.0/3.0)*pi*theA*r0*r0*r0; << 124 G4double V0 = (4.0/3.0)*pi*theA*r03;
145 125
146 G4double MeanA = 0.0; 126 G4double MeanA = 0.0;
147 127
148 _MeanMultiplicity = 0.0; 128 _MeanMultiplicity = 0.0;
149 129
>> 130
150 G4int n = 1; 131 G4int n = 1;
151 for (std::vector<G4VStatMFMacroCluster*>::it << 132 for (std::vector<G4VStatMFMacroCluster*>::iterator i = _theClusters->begin();
152 i != _theClusters->end(); ++i) 133 i != _theClusters->end(); ++i)
153 { 134 {
154 G4double multip = (*i)->CalcMeanMultiplic << 135 G4double multip = (*i)->CalcMeanMultiplicity(V0*_Kappa,mu,_ChemPotentialNu,_MeanTemperature);
155 _MeanTemperature); << 136 MeanA += multip*static_cast<G4double>(n++);
156 MeanA += multip*(n++); <<
157 _MeanMultiplicity += multip; 137 _MeanMultiplicity += multip;
158 } 138 }
159 139
160 return MeanA; 140 return MeanA;
161 } 141 }
162 142