Geant4_10
G4StatMFMacroChemicalPotential.cc
Go to the documentation of this file.
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 // $Id: G4StatMFMacroChemicalPotential.cc 67983 2013-03-13 10:42:03Z gcosmo $
28 //
29 // Hadronic Process: Nuclear De-excitations
30 // by V. Lara
31 
33 #include "G4PhysicalConstants.hh"
34 
35 // operators definitions
37 G4StatMFMacroChemicalPotential::operator=(const G4StatMFMacroChemicalPotential & )
38 {
39  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroChemicalPotential::operator= meant to not be accessable");
40  return *this;
41 }
42 
43 G4bool G4StatMFMacroChemicalPotential::operator==(const G4StatMFMacroChemicalPotential & ) const
44 {
45  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroChemicalPotential::operator== meant to not be accessable");
46  return false;
47 }
48 
49 
50 G4bool G4StatMFMacroChemicalPotential::operator!=(const G4StatMFMacroChemicalPotential & ) const
51 {
52  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroChemicalPotential::operator!= meant to not be accessable");
53  return true;
54 }
55 
56 
57 
58 
60  // Calculate Chemical potential \nu
61 {
63  (1.0-1.0/std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1.0/3.0));
64 
65  // Initial value for _ChemPotentialNu
66  _ChemPotentialNu = (theZ/theA)*(8.0*G4StatMFParameters::GetGamma0()+2.0*CP*std::pow(theA,2./3.)) -
68 
69 
70  G4double ChemPa = _ChemPotentialNu;
71  G4double ChemPb = 0.5*_ChemPotentialNu;
72 
73  G4double fChemPa = this->operator()(ChemPa);
74  G4double fChemPb = this->operator()(ChemPb);
75 
76  if (fChemPa*fChemPb > 0.0) {
77  // bracketing the solution
78  if (fChemPa < 0.0) {
79  do {
80  ChemPb -= 1.5*std::abs(ChemPb-ChemPa);
81  fChemPb = this->operator()(ChemPb);
82  } while (fChemPb < 0.0);
83  } else {
84  do {
85  ChemPb += 1.5*std::abs(ChemPb-ChemPa);
86  fChemPb = this->operator()(ChemPb);
87  } while (fChemPb > 0.0);
88  }
89  }
90 
93  theSolver->SetIntervalLimits(ChemPa,ChemPb);
94  // if (!theSolver->Crenshaw(*this))
95  if (!theSolver->Brent(*this)){
96  G4cerr <<"G4StatMFMacroChemicalPotential:"<<" ChemPa="<<ChemPa<<" ChemPb="<<ChemPb<< G4endl;
97  G4cerr <<"G4StatMFMacroChemicalPotential:"<<" fChemPa="<<fChemPa<<" fChemPb="<<fChemPb<< G4endl;
98  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroChemicalPotential::CalcChemicalPotentialNu: I couldn't find the root.");
99  }
100  _ChemPotentialNu = theSolver->GetRoot();
101  delete theSolver;
102  return _ChemPotentialNu;
103 }
104 
105 
106 
107 G4double G4StatMFMacroChemicalPotential::CalcMeanZ(const G4double nu)
108 {
109  std::vector<G4VStatMFMacroCluster*>::iterator i;
110  for (i= _theClusters->begin()+1; i != _theClusters->end(); ++i)
111  {
112  (*i)->CalcZARatio(nu);
113  }
114  CalcChemicalPotentialMu(nu);
115  // This is important, the Z over A ratio for proton and neutron depends on the
116  // chemical potential Mu, while for the first guess for Chemical potential mu
117  // some values of Z over A ratio. This is the reason for that.
118  (*_theClusters->begin())->CalcZARatio(nu);
119 
120  G4double MeanZ = 0.0;
121  G4int n = 1;
122  for (i = _theClusters->begin(); i != _theClusters->end(); ++i)
123  {
124  MeanZ += static_cast<G4double>(n++) *
125  (*i)->GetZARatio() *
126  (*i)->GetMeanMultiplicity();
127  }
128  return MeanZ;
129 }
130 
131 
132 void G4StatMFMacroChemicalPotential::CalcChemicalPotentialMu(const G4double nu)
133  // Calculate Chemical potential \mu
134  // For that is necesary to calculate mean multiplicities
135 {
136  G4StatMFMacroMultiplicity * theMultip = new
137  G4StatMFMacroMultiplicity(theA,_Kappa,_MeanTemperature,nu,_theClusters);
138 
139  _ChemPotentialMu = theMultip->CalcChemicalPotentialMu();
140  _MeanMultiplicity = theMultip->GetMeanMultiplicity();
141 
142  delete theMultip;
143 
144  return;
145 
146 }
static G4double GetGamma0()
static G4double GetKappaCoulomb()
tuple elm_coupling
Definition: hepunit.py:286
void SetIntervalLimits(const G4double Limit1, const G4double Limit2)
int G4int
Definition: G4Types.hh:78
G4bool Brent(Function &theFunction)
static G4double Getr0()
Char_t n[5]
bool G4bool
Definition: G4Types.hh:79
G4double GetMeanMultiplicity(void) const
G4double GetRoot(void) const
Definition: G4Solver.hh:77
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4GLOB_DLL std::ostream G4cerr