Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // >> 27 // $Id: G4StatMFMacroTemperature.cc 67983 2013-03-13 10:42:03Z gcosmo $ 27 // 28 // 28 // Hadronic Process: Nuclear De-excitations 29 // Hadronic Process: Nuclear De-excitations 29 // by V. Lara 30 // by V. Lara 30 // 31 // 31 // Modified: 32 // Modified: 32 // 25.07.08 I.Pshenichnov (in collaboration wi 33 // 25.07.08 I.Pshenichnov (in collaboration with Alexander Botvina and Igor 33 // Mishustin (FIAS, Frankfurt, INR, M 34 // Mishustin (FIAS, Frankfurt, INR, Moscow and Kurchatov Institute, 34 // Moscow, pshenich@fias.uni-frankfur 35 // Moscow, pshenich@fias.uni-frankfurt.de) make algorithm closer to 35 // original MF model 36 // original MF model 36 // 16.04.10 V.Ivanchenko improved logic of sol << 37 // 16.04.10 V.Ivanchenko improved logic of solving equation for tempetature 37 // to protect code from rare unwanted 38 // to protect code from rare unwanted exception; moved constructor 38 // and destructor to source 39 // and destructor to source 39 // 28.10.10 V.Ivanchenko defined members in co 40 // 28.10.10 V.Ivanchenko defined members in constructor and cleaned up 40 41 41 #include "G4StatMFMacroTemperature.hh" 42 #include "G4StatMFMacroTemperature.hh" 42 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicalConstants.hh" 43 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh" 44 #include "G4Pow.hh" << 45 45 46 G4StatMFMacroTemperature::G4StatMFMacroTempera 46 G4StatMFMacroTemperature::G4StatMFMacroTemperature(const G4double anA, const G4double aZ, 47 const G4double ExEnergy, const G4double Free 47 const G4double ExEnergy, const G4double FreeE0, const G4double kappa, 48 std::vector<G4VStatMFMacroCluster*> * Cluste 48 std::vector<G4VStatMFMacroCluster*> * ClusterVector) : 49 theA(anA), 49 theA(anA), 50 theZ(aZ), 50 theZ(aZ), 51 _ExEnergy(ExEnergy), 51 _ExEnergy(ExEnergy), 52 _FreeInternalE0(FreeE0), 52 _FreeInternalE0(FreeE0), 53 _Kappa(kappa), 53 _Kappa(kappa), 54 _MeanMultiplicity(0.0), 54 _MeanMultiplicity(0.0), 55 _MeanTemperature(0.0), 55 _MeanTemperature(0.0), 56 _ChemPotentialMu(0.0), 56 _ChemPotentialMu(0.0), 57 _ChemPotentialNu(0.0), 57 _ChemPotentialNu(0.0), 58 _MeanEntropy(0.0), 58 _MeanEntropy(0.0), 59 _theClusters(ClusterVector) 59 _theClusters(ClusterVector) 60 {} 60 {} 61 61 62 G4StatMFMacroTemperature::~G4StatMFMacroTemper 62 G4StatMFMacroTemperature::~G4StatMFMacroTemperature() 63 {} 63 {} 64 64 >> 65 65 G4double G4StatMFMacroTemperature::CalcTempera 66 G4double G4StatMFMacroTemperature::CalcTemperature(void) 66 { 67 { 67 // Inital guess for the interval of the ense << 68 // Inital guess for the interval of the ensemble temperature values 68 G4double Ta = 0.5; << 69 G4double Ta = 0.5; 69 G4double Tb = std::max(std::sqrt(_ExEnergy/( << 70 G4double Tb = std::max(std::sqrt(_ExEnergy/(theA*0.12)),0.01*MeV); 70 71 71 G4double fTa = this->operator()(Ta); << 72 G4double fTa = this->operator()(Ta); 72 G4double fTb = this->operator()(Tb); << 73 G4double fTb = this->operator()(Tb); 73 74 74 // Bracketing the solution << 75 // Bracketing the solution 75 // T should be greater than 0. << 76 // T should be greater than 0. 76 // The interval is [Ta,Tb] << 77 // The interval is [Ta,Tb] 77 // We start with a value for Ta = 0.5 MeV << 78 // We start with a value for Ta = 0.5 MeV 78 // it should be enough to have fTa > 0 If it << 79 // it should be enough to have fTa > 0 If it isn't 79 // the case, we decrease Ta. But carefully, << 80 // the case, we decrease Ta. But carefully, because 80 // fTa growes very fast when Ta is near 0 an << 81 // fTa growes very fast when Ta is near 0 and we could have 81 // an overflow. << 82 // an overflow. 82 << 83 83 G4int iterations = 0; << 84 G4int iterations = 0; 84 // Loop checking, 05-Aug-2015, Vladimir Ivan << 85 while (fTa < 0.0 && ++iterations < 10) { 85 while (fTa < 0.0 && ++iterations < 10) { << 86 Ta -= 0.5*Ta; 86 Ta -= 0.5*Ta; << 87 fTa = this->operator()(Ta); 87 fTa = this->operator()(Ta); << 88 } 88 } << 89 // Usually, fTb will be less than 0, but if it is not the case: 89 // Usually, fTb will be less than 0, but if << 90 iterations = 0; 90 iterations = 0; << 91 while (fTa*fTb > 0.0 && iterations++ < 10) { 91 // Loop checking, 05-Aug-2015, Vladimir Ivan << 92 Tb += 2.*std::fabs(Tb-Ta); 92 while (fTa*fTb > 0.0 && iterations++ < 10) { << 93 fTb = this->operator()(Tb); 93 Tb += 2.*std::fabs(Tb-Ta); << 94 } 94 fTb = this->operator()(Tb); << 95 } << 96 95 97 if (fTa*fTb > 0.0) { << 96 if (fTa*fTb > 0.0) { 98 G4cerr <<"G4StatMFMacroTemperature:"<<" Ta << 97 G4cerr <<"G4StatMFMacroTemperature:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl; 99 G4cerr <<"G4StatMFMacroTemperature:"<<" fT << 98 G4cerr <<"G4StatMFMacroTemperature:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl; 100 throw G4HadronicException(__FILE__, __LINE << 99 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't bracket the solution."); 101 } << 100 } 102 << 103 G4Solver<G4StatMFMacroTemperature> * theSolv << 104 new G4Solver<G4StatMFMacroTemperature>(100 << 105 theSolver->SetIntervalLimits(Ta,Tb); << 106 if (!theSolver->Crenshaw(*this)){ << 107 G4cout <<"G4StatMFMacroTemperature, Crensh << 108 <<Ta<<" Tb="<<Tb<< G4endl; << 109 G4cout <<"G4StatMFMacroTemperature, Crensh << 110 <<fTa<<" fTb="<<fTb<< G4endl; << 111 } << 112 _MeanTemperature = theSolver->GetRoot(); << 113 G4double FunctionValureAtRoot = this->opera << 114 delete theSolver; << 115 << 116 // Verify if the root is found and it is ind << 117 // say, between 1 and 50 MeV, otherwise try << 118 if (std::fabs(FunctionValureAtRoot) > 5.e-2) << 119 if (_MeanTemperature < 1. || _MeanTemperat << 120 G4cout << "Crenshaw method failed; funct << 121 << " solution? = " << _MeanTemperature << 122 G4Solver<G4StatMFMacroTemperature> * the << 123 new G4Solver<G4StatMFMacroTemperature>(200,1 << 124 theSolverBrent->SetIntervalLimits(Ta,Tb) << 125 if (!theSolverBrent->Brent(*this)){ << 126 G4cout <<"G4StatMFMacroTemperature, Brent me << 127 <<" Ta="<<Ta<<" Tb="<<Tb<< G4endl; << 128 G4cout <<"G4StatMFMacroTemperature, Brent me << 129 <<" fTa="<<fTa<<" fTb="<<fTb<< G4endl << 130 throw G4HadronicException(__FILE__, __LINE__ << 131 } << 132 101 133 _MeanTemperature = theSolverBrent->GetRo << 102 G4Solver<G4StatMFMacroTemperature> * theSolver = new G4Solver<G4StatMFMacroTemperature>(100,1.e-4); 134 FunctionValureAtRoot = this->operator() << 103 theSolver->SetIntervalLimits(Ta,Tb); 135 delete theSolverBrent; << 104 if (!theSolver->Crenshaw(*this)){ >> 105 G4cout <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl; >> 106 G4cout <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl; 136 } 107 } 137 if (std::abs(FunctionValureAtRoot) > 5.e-2 << 108 _MeanTemperature = theSolver->GetRoot(); 138 G4cout << "Brent method failed; function << 109 G4double FunctionValureAtRoot = this->operator()(_MeanTemperature); 139 << " solution? = " << _MeanTemperature << 110 delete theSolver; 140 throw G4HadronicException(__FILE__, __LI << 111 >> 112 // Verify if the root is found and it is indeed within the physical domain, >> 113 // say, between 1 and 50 MeV, otherwise try Brent method: >> 114 if (std::fabs(FunctionValureAtRoot) > 5.e-2) { >> 115 if (_MeanTemperature < 1. || _MeanTemperature > 50.) { >> 116 G4cout << "Crenshaw method failed; function = " << FunctionValureAtRoot >> 117 << " solution? = " << _MeanTemperature << " MeV " << G4endl; >> 118 G4Solver<G4StatMFMacroTemperature> * theSolverBrent = new G4Solver<G4StatMFMacroTemperature>(200,1.e-3); >> 119 theSolverBrent->SetIntervalLimits(Ta,Tb); >> 120 if (!theSolverBrent->Brent(*this)){ >> 121 G4cout <<"G4StatMFMacroTemperature, Brent method failed:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl; >> 122 G4cout <<"G4StatMFMacroTemperature, Brent method failed:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl; >> 123 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root with any method."); >> 124 } >> 125 >> 126 _MeanTemperature = theSolverBrent->GetRoot(); >> 127 FunctionValureAtRoot = this->operator()(_MeanTemperature); >> 128 delete theSolverBrent; >> 129 } >> 130 if (std::abs(FunctionValureAtRoot) > 5.e-2) { >> 131 //if (_MeanTemperature < 1. || _MeanTemperature > 50. || std::abs(FunctionValureAtRoot) > 5.e-2) { >> 132 G4cout << "Brent method failed; function = " << FunctionValureAtRoot << " solution? = " << _MeanTemperature << " MeV " << G4endl; >> 133 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root with any method."); >> 134 } 141 } 135 } 142 } << 136 //G4cout << "G4StatMFMacroTemperature::CalcTemperature: function = " << FunctionValureAtRoot 143 //G4cout << "G4StatMFMacroTemperature::CalcT << 137 // << " T(MeV)= " << _MeanTemperature << G4endl; 144 //<< FunctionValureAtRoot << 138 return _MeanTemperature; 145 // << " T(MeV)= " << _MeanTemperature << << 146 return _MeanTemperature; << 147 } 139 } 148 140 >> 141 >> 142 149 G4double G4StatMFMacroTemperature::FragsExcitE 143 G4double G4StatMFMacroTemperature::FragsExcitEnergy(const G4double T) 150 // Calculates excitation energy per nucleon an << 144 // Calculates excitation energy per nucleon and summed fragment multiplicity and entropy 151 // multiplicity and entropy << 152 { 145 { 153 // Model Parameters << 146 154 G4Pow* g4calc = G4Pow::GetInstance(); << 147 // Model Parameters 155 G4double R0 = G4StatMFParameters::Getr0()*g4 << 148 G4double R0 = G4StatMFParameters::Getr0()*std::pow(theA,1./3.); 156 G4double R = R0*g4calc->A13(1.0+G4StatMFPara << 149 G4double R = R0*std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(), 1./3.); 157 G4double FreeVol = _Kappa*(4.*pi/3.)*R0*R0*R << 150 G4double FreeVol = _Kappa*(4.*pi/3.)*R0*R0*R0; 158 151 159 // Calculate Chemical potentials << 152 160 CalcChemicalPotentialNu(T); << 153 // Calculate Chemical potentials >> 154 CalcChemicalPotentialNu(T); 161 155 162 156 163 // Average total fragment energy << 157 // Average total fragment energy 164 G4double AverageEnergy = 0.0; << 158 G4double AverageEnergy = 0.0; 165 std::vector<G4VStatMFMacroCluster*>::iterato << 159 std::vector<G4VStatMFMacroCluster*>::iterator i; 166 for (i = _theClusters->begin(); i != _theCl << 160 for (i = _theClusters->begin(); i != _theClusters->end(); ++i) 167 { << 161 { 168 AverageEnergy += (*i)->GetMeanMultiplici << 162 AverageEnergy += (*i)->GetMeanMultiplicity() * (*i)->CalcEnergy(T); 169 } << 163 } 170 164 171 // Add Coulomb energy << 165 // Add Coulomb energy 172 AverageEnergy += 0.6*elm_coupling*theZ*theZ/ << 166 AverageEnergy += (3./5.)*elm_coupling*theZ*theZ/R; 173 167 174 // Calculate mean entropy << 168 // Calculate mean entropy 175 _MeanEntropy = 0.0; << 169 _MeanEntropy = 0.0; 176 for (i = _theClusters->begin(); i != _theClu << 170 for (i = _theClusters->begin(); i != _theClusters->end(); ++i) 177 { << 171 { 178 _MeanEntropy += (*i)->CalcEntropy(T,Free << 172 _MeanEntropy += (*i)->CalcEntropy(T,FreeVol); 179 } << 173 } >> 174 >> 175 // Excitation energy per nucleon >> 176 return AverageEnergy - _FreeInternalE0; 180 177 181 // Excitation energy per nucleon << 182 return AverageEnergy - _FreeInternalE0; << 183 } 178 } 184 179 >> 180 185 void G4StatMFMacroTemperature::CalcChemicalPot 181 void G4StatMFMacroTemperature::CalcChemicalPotentialNu(const G4double T) 186 // Calculates the chemical potential \nu << 182 // Calculates the chemical potential \nu >> 183 187 { 184 { 188 G4StatMFMacroChemicalPotential * theChemPot << 185 G4StatMFMacroChemicalPotential * theChemPot = new 189 G4StatMFMacroChemicalPotential(theA,theZ,_ << 186 G4StatMFMacroChemicalPotential(theA,theZ,_Kappa,T,_theClusters); >> 187 190 188 191 _ChemPotentialNu = theChemPot->CalcChemicalP << 189 _ChemPotentialNu = theChemPot->CalcChemicalPotentialNu(); 192 _ChemPotentialMu = theChemPot->GetChemicalPo << 190 _ChemPotentialMu = theChemPot->GetChemicalPotentialMu(); 193 _MeanMultiplicity = theChemPot->GetMeanMulti << 191 _MeanMultiplicity = theChemPot->GetMeanMultiplicity(); 194 delete theChemPot; << 192 >> 193 delete theChemPot; 195 194 196 return; << 195 return; >> 196 197 } 197 } 198 198 199 199 200 200