ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4StatMFMacroMultiplicity.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4StatMFMacroMultiplicity.cc
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 //
28 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara
30 //
31 // Modified:
32 // 25.07.08 I.Pshenichnov (in collaboration with Alexander Botvina and Igor
33 // Mishustin (FIAS, Frankfurt, INR, Moscow and Kurchatov Institute,
34 // Moscow, pshenich@fias.uni-frankfurt.de) additional checks in
35 // solver of equation for the chemical potential
36 
38 #include "G4PhysicalConstants.hh"
39 #include "G4Pow.hh"
40 
41 // operators definitions
44 {
45  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::operator= meant to not be accessible");
46  return *this;
47 }
48 
50 {
51  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::operator== meant to not be accessible");
52  return false;
53 }
54 
55 
57 {
58  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::operator!= meant to not be accessible");
59  return true;
60 }
61 
63  // Calculate Chemical potential \mu
64  // For that is necesary to calculate mean multiplicities
65 {
66  G4Pow* g4calc = G4Pow::GetInstance();
68 
69  // starting value for chemical potential \mu
70  // it is the derivative of F(T,V)-\nu*Z w.r.t. Af in Af=5
71  G4double ZA5 = _theClusters->operator[](4)->GetZARatio();
72  G4double ILD5 = _theClusters->operator[](4)->GetInvLevelDensity();
75  _ChemPotentialNu*ZA5 +
76  G4StatMFParameters::GetGamma0()*(1.0-2.0*ZA5)*(1.0-2.0*ZA5) +
77  (2.0/3.0)*G4StatMFParameters::Beta(_MeanTemperature)/g4calc->Z13(5) +
78  (5.0/3.0)*CP*ZA5*ZA5*g4calc->Z23(5) -
79  1.5*_MeanTemperature/5.0;
80 
81  G4double ChemPa = _ChemPotentialMu;
82  if (ChemPa/_MeanTemperature > 10.0) ChemPa = 10.0*_MeanTemperature;
83  G4double ChemPb = ChemPa - 0.5*std::abs(ChemPa);
84 
85  G4double fChemPa = this->operator()(ChemPa);
86  G4double fChemPb = this->operator()(ChemPb);
87 
88  // Set the precision level for locating the root.
89  // If the root is inside this interval, then it's done!
90  const G4double intervalWidth = 1.e-4;
91 
92  // bracketing the solution
93  G4int iterations = 0;
94  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
95  while (fChemPa*fChemPb > 0.0 && iterations < 100)
96  {
97  iterations++;
98  if (std::abs(fChemPa) <= std::abs(fChemPb))
99  {
100  ChemPa += 0.6*(ChemPa-ChemPb);
101  fChemPa = this->operator()(ChemPa);
102  }
103  else
104  {
105  ChemPb += 0.6*(ChemPb-ChemPa);
106  fChemPb = this->operator()(ChemPb);
107  }
108  }
109 
110  if (fChemPa*fChemPb > 0.0) // the bracketing failed, complain
111  {
112  G4cout <<"G4StatMFMacroMultiplicity:"<<" ChemPa="<<ChemPa
113  <<" ChemPb="<<ChemPb<< G4endl;
114  G4cout <<"G4StatMFMacroMultiplicity:"<<" fChemPa="<<fChemPa
115  <<" fChemPb="<<fChemPb<< G4endl;
116  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::CalcChemicalPotentialMu: I couldn't bracket the root.");
117  }
118  else if (fChemPa*fChemPb < 0.0 && std::abs(ChemPa-ChemPb) > intervalWidth)
119  {
121  new G4Solver<G4StatMFMacroMultiplicity>(100,intervalWidth);
122  theSolver->SetIntervalLimits(ChemPa,ChemPb);
123  // if (!theSolver->Crenshaw(*this))
124  if (!theSolver->Brent(*this))
125  {
126  G4cout <<"G4StatMFMacroMultiplicity:"<<" ChemPa="<<ChemPa
127  <<" ChemPb="<<ChemPb<< G4endl;
128  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::CalcChemicalPotentialMu: I couldn't find the root.");
129  }
130  _ChemPotentialMu = theSolver->GetRoot();
131  delete theSolver;
132  }
133  else // the root is within the interval, which is shorter then the precision level - all done
134  {
135  _ChemPotentialMu = ChemPa;
136  }
137 
138  return _ChemPotentialMu;
139 }
140 
142 {
144  G4double V0 = (4.0/3.0)*pi*theA*r0*r0*r0;
145 
146  G4double MeanA = 0.0;
147 
148  _MeanMultiplicity = 0.0;
149 
150  G4int n = 1;
151  for (std::vector<G4VStatMFMacroCluster*>::iterator i = _theClusters->begin();
152  i != _theClusters->end(); ++i)
153  {
154  G4double multip = (*i)->CalcMeanMultiplicity(V0*_Kappa,mu,_ChemPotentialNu,
156  MeanA += multip*(n++);
157  _MeanMultiplicity += multip;
158  }
159 
160  return MeanA;
161 }