Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMacroTemperature.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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) make algorithm closer to
 35 //          original MF model
 36 // 16.04.10 V.Ivanchenko improved logic of solving equation for temperature
 37 //          to protect code from rare unwanted exception; moved constructor 
 38 //          and destructor to source  
 39 // 28.10.10 V.Ivanchenko defined members in constructor and cleaned up
 40 
 41 #include "G4StatMFMacroTemperature.hh"
 42 #include "G4PhysicalConstants.hh"
 43 #include "G4SystemOfUnits.hh"
 44 #include "G4Pow.hh"
 45 
 46 G4StatMFMacroTemperature::G4StatMFMacroTemperature(const G4double anA, const G4double aZ, 
 47   const G4double ExEnergy, const G4double FreeE0, const G4double kappa, 
 48   std::vector<G4VStatMFMacroCluster*> * ClusterVector) :
 49   theA(anA),
 50   theZ(aZ),
 51   _ExEnergy(ExEnergy),
 52   _FreeInternalE0(FreeE0),
 53   _Kappa(kappa),
 54   _MeanMultiplicity(0.0),
 55   _MeanTemperature(0.0),
 56   _ChemPotentialMu(0.0),
 57   _ChemPotentialNu(0.0),
 58   _MeanEntropy(0.0),
 59   _theClusters(ClusterVector) 
 60 {}
 61   
 62 G4StatMFMacroTemperature::~G4StatMFMacroTemperature() 
 63 {}
 64 
 65 G4double G4StatMFMacroTemperature::CalcTemperature(void) 
 66 {
 67   // Inital guess for the interval of the ensemble temperature values
 68   G4double Ta = 0.5; 
 69   G4double Tb = std::max(std::sqrt(_ExEnergy/(theA*0.12)),0.01*MeV);
 70     
 71   G4double fTa = this->operator()(Ta); 
 72   G4double fTb = this->operator()(Tb); 
 73 
 74   // Bracketing the solution
 75   // T should be greater than 0.
 76   // The interval is [Ta,Tb]
 77   // We start with a value for Ta = 0.5 MeV
 78   // it should be enough to have fTa > 0 If it isn't 
 79   // the case, we decrease Ta. But carefully, because 
 80   // fTa growes very fast when Ta is near 0 and we could have
 81   // an overflow.
 82 
 83   G4int iterations = 0;  
 84   // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
 85   while (fTa < 0.0 && ++iterations < 10) {
 86     Ta -= 0.5*Ta;
 87     fTa = this->operator()(Ta);
 88   }
 89   // Usually, fTb will be less than 0, but if it is not the case: 
 90   iterations = 0;  
 91   // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
 92   while (fTa*fTb > 0.0 && iterations++ < 10) {
 93     Tb += 2.*std::fabs(Tb-Ta);
 94     fTb = this->operator()(Tb);
 95   }
 96   
 97   if (fTa*fTb > 0.0) {
 98     G4cerr <<"G4StatMFMacroTemperature:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
 99     G4cerr <<"G4StatMFMacroTemperature:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl;
100     throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't bracket the solution.");
101   }
102 
103   G4Solver<G4StatMFMacroTemperature> * theSolver = 
104     new G4Solver<G4StatMFMacroTemperature>(100,1.e-4);
105   theSolver->SetIntervalLimits(Ta,Tb);
106   if (!theSolver->Crenshaw(*this)){ 
107     G4cout <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" Ta="
108      <<Ta<<" Tb="<<Tb<< G4endl;
109     G4cout <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" fTa="
110      <<fTa<<" fTb="<<fTb<< G4endl;
111   }
112   _MeanTemperature = theSolver->GetRoot();
113   G4double FunctionValureAtRoot =  this->operator()(_MeanTemperature);
114   delete  theSolver;
115 
116   // Verify if the root is found and it is indeed within the physical domain, 
117   // say, between 1 and 50 MeV, otherwise try Brent method:
118   if (std::fabs(FunctionValureAtRoot) > 5.e-2) {
119     if (_MeanTemperature < 1. || _MeanTemperature > 50.) {
120       G4cout << "Crenshaw method failed; function = " << FunctionValureAtRoot 
121        << " solution? = " << _MeanTemperature << " MeV " << G4endl;
122       G4Solver<G4StatMFMacroTemperature> * theSolverBrent = 
123   new G4Solver<G4StatMFMacroTemperature>(200,1.e-3);
124       theSolverBrent->SetIntervalLimits(Ta,Tb);
125       if (!theSolverBrent->Brent(*this)){
126   G4cout <<"G4StatMFMacroTemperature, Brent method failed:"
127          <<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
128   G4cout <<"G4StatMFMacroTemperature, Brent method failed:"
129          <<" fTa="<<fTa<<" fTb="<<fTb<< G4endl; 
130   throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root with any method.");
131       }
132 
133       _MeanTemperature = theSolverBrent->GetRoot();
134       FunctionValureAtRoot =  this->operator()(_MeanTemperature);
135       delete theSolverBrent;
136     }
137     if (std::abs(FunctionValureAtRoot) > 5.e-2) {
138       G4cout << "Brent method failed; function = " << FunctionValureAtRoot 
139        << " solution? = " << _MeanTemperature << " MeV " << G4endl;
140       throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root with any method.");
141     }
142   }
143   //G4cout << "G4StatMFMacroTemperature::CalcTemperature: function = " 
144   //<< FunctionValureAtRoot 
145   //     << " T(MeV)= " << _MeanTemperature << G4endl;
146   return _MeanTemperature;
147 }
148 
149 G4double G4StatMFMacroTemperature::FragsExcitEnergy(const G4double T)
150 // Calculates excitation energy per nucleon and summed fragment 
151 // multiplicity and entropy
152 {
153   // Model Parameters
154   G4Pow* g4calc = G4Pow::GetInstance();
155   G4double R0 = G4StatMFParameters::Getr0()*g4calc->Z13(theA);
156   G4double R = R0*g4calc->A13(1.0+G4StatMFParameters::GetKappaCoulomb());
157   G4double FreeVol = _Kappa*(4.*pi/3.)*R0*R0*R0; 
158  
159   // Calculate Chemical potentials
160   CalcChemicalPotentialNu(T);
161 
162 
163   // Average total fragment energy
164   G4double AverageEnergy = 0.0;
165   std::vector<G4VStatMFMacroCluster*>::iterator i;
166   for (i =  _theClusters->begin(); i != _theClusters->end(); ++i) 
167     {
168       AverageEnergy += (*i)->GetMeanMultiplicity() * (*i)->CalcEnergy(T);
169     }
170     
171   // Add Coulomb energy     
172   AverageEnergy += 0.6*elm_coupling*theZ*theZ/R;    
173     
174   // Calculate mean entropy
175   _MeanEntropy = 0.0;
176   for (i = _theClusters->begin(); i != _theClusters->end(); ++i) 
177     {
178       _MeanEntropy += (*i)->CalcEntropy(T,FreeVol); 
179     }
180 
181   // Excitation energy per nucleon
182   return AverageEnergy - _FreeInternalE0;
183 }
184 
185 void G4StatMFMacroTemperature::CalcChemicalPotentialNu(const G4double T)
186 // Calculates the chemical potential \nu 
187 {
188   G4StatMFMacroChemicalPotential * theChemPot = new
189     G4StatMFMacroChemicalPotential(theA,theZ,_Kappa,T,_theClusters);
190 
191   _ChemPotentialNu = theChemPot->CalcChemicalPotentialNu();
192   _ChemPotentialMu = theChemPot->GetChemicalPotentialMu();
193   _MeanMultiplicity = theChemPot->GetMeanMultiplicity();    
194   delete theChemPot;
195         
196   return;
197 }
198 
199 
200