G4StatMFMacroChemicalPotential.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // Hadronic Process: Nuclear De-excitations
00030 // by V. Lara
00031 
00032 #include "G4StatMFMacroChemicalPotential.hh"
00033 #include "G4PhysicalConstants.hh"
00034 
00035 // operators definitions
00036 G4StatMFMacroChemicalPotential & 
00037 G4StatMFMacroChemicalPotential::operator=(const G4StatMFMacroChemicalPotential & ) 
00038 {
00039     throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroChemicalPotential::operator= meant to not be accessable");
00040     return *this;
00041 }
00042 
00043 G4bool G4StatMFMacroChemicalPotential::operator==(const G4StatMFMacroChemicalPotential & ) const 
00044 {
00045     throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroChemicalPotential::operator== meant to not be accessable");
00046     return false;
00047 }
00048 
00049 
00050 G4bool G4StatMFMacroChemicalPotential::operator!=(const G4StatMFMacroChemicalPotential & ) const 
00051 {
00052     throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroChemicalPotential::operator!= meant to not be accessable");
00053     return true;
00054 }
00055 
00056 
00057 
00058 
00059 G4double G4StatMFMacroChemicalPotential::CalcChemicalPotentialNu(void) 
00060     //  Calculate Chemical potential \nu
00061 {
00062     G4double CP = ((3./5.)*elm_coupling/G4StatMFParameters::Getr0())*
00063         (1.0-1.0/std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1.0/3.0));
00064 
00065     // Initial value for _ChemPotentialNu       
00066     _ChemPotentialNu = (theZ/theA)*(8.0*G4StatMFParameters::GetGamma0()+2.0*CP*std::pow(theA,2./3.)) -
00067         4.0*G4StatMFParameters::GetGamma0();
00068                 
00069 
00070     G4double ChemPa = _ChemPotentialNu;
00071     G4double ChemPb = 0.5*_ChemPotentialNu;
00072     
00073     G4double fChemPa = this->operator()(ChemPa); 
00074     G4double fChemPb = this->operator()(ChemPb); 
00075 
00076     if (fChemPa*fChemPb > 0.0) {    
00077         // bracketing the solution
00078         if (fChemPa < 0.0) {
00079             do {
00080                 ChemPb -= 1.5*std::abs(ChemPb-ChemPa);
00081                 fChemPb = this->operator()(ChemPb);   
00082             } while (fChemPb < 0.0);
00083         } else {
00084             do {
00085                 ChemPb += 1.5*std::abs(ChemPb-ChemPa);
00086                 fChemPb = this->operator()(ChemPb);
00087             } while (fChemPb > 0.0);
00088         }
00089     }
00090 
00091     G4Solver<G4StatMFMacroChemicalPotential> * theSolver =
00092       new G4Solver<G4StatMFMacroChemicalPotential>(100,1.e-4);
00093     theSolver->SetIntervalLimits(ChemPa,ChemPb);
00094     //    if (!theSolver->Crenshaw(*this)) 
00095     if (!theSolver->Brent(*this)){
00096       G4cerr <<"G4StatMFMacroChemicalPotential:"<<" ChemPa="<<ChemPa<<" ChemPb="<<ChemPb<< G4endl;
00097       G4cerr <<"G4StatMFMacroChemicalPotential:"<<" fChemPa="<<fChemPa<<" fChemPb="<<fChemPb<< G4endl;
00098       throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroChemicalPotential::CalcChemicalPotentialNu: I couldn't find the root.");
00099     }
00100     _ChemPotentialNu = theSolver->GetRoot();
00101     delete theSolver;
00102     return _ChemPotentialNu;
00103 }
00104 
00105 
00106 
00107 G4double G4StatMFMacroChemicalPotential::CalcMeanZ(const G4double nu)
00108 {
00109   std::vector<G4VStatMFMacroCluster*>::iterator i;
00110   for (i= _theClusters->begin()+1; i != _theClusters->end(); ++i) 
00111     { 
00112       (*i)->CalcZARatio(nu);
00113     }
00114   CalcChemicalPotentialMu(nu);
00115   // This is important, the Z over A ratio for proton and neutron depends on the 
00116   // chemical potential Mu, while for the first guess for Chemical potential mu 
00117   // some values of Z over A ratio. This is the reason for that.
00118   (*_theClusters->begin())->CalcZARatio(nu);
00119   
00120   G4double MeanZ = 0.0;
00121   G4int n = 1;
00122   for (i = _theClusters->begin(); i != _theClusters->end(); ++i)
00123     {
00124       MeanZ += static_cast<G4double>(n++) *
00125         (*i)->GetZARatio() *
00126         (*i)->GetMeanMultiplicity(); 
00127     }
00128   return MeanZ;
00129 }
00130 
00131 
00132 void G4StatMFMacroChemicalPotential::CalcChemicalPotentialMu(const G4double nu)
00133   //    Calculate Chemical potential \mu
00134   // For that is necesary to calculate mean multiplicities
00135 {
00136   G4StatMFMacroMultiplicity * theMultip = new
00137     G4StatMFMacroMultiplicity(theA,_Kappa,_MeanTemperature,nu,_theClusters);
00138   
00139   _ChemPotentialMu = theMultip->CalcChemicalPotentialMu();
00140   _MeanMultiplicity = theMultip->GetMeanMultiplicity();
00141   
00142   delete theMultip;
00143   
00144   return;
00145   
00146 }

Generated on Mon May 27 17:49:53 2013 for Geant4 by  doxygen 1.4.7