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$ 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 #include "G4StatMF.hh" 32 #include "G4StatMF.hh" 32 #include "G4PhysicalConstants.hh" 33 #include "G4PhysicalConstants.hh" 33 #include "G4SystemOfUnits.hh" 34 #include "G4SystemOfUnits.hh" 34 #include "G4Pow.hh" 35 #include "G4Pow.hh" 35 #include "G4PhysicsModelCatalog.hh" << 36 36 37 G4StatMF::G4StatMF() << 38 { << 39 _secID = G4PhysicsModelCatalog::GetModelID(" << 40 } << 41 37 42 G4StatMF::~G4StatMF() {} << 38 // Default constructor >> 39 G4StatMF::G4StatMF() : _theEnsemble(0) {} >> 40 >> 41 >> 42 // Destructor >> 43 G4StatMF::~G4StatMF() {} //{if (_theEnsemble != 0) delete _theEnsemble;} 43 44 44 G4FragmentVector* G4StatMF::BreakItUp(const G4 << 45 >> 46 G4FragmentVector * G4StatMF::BreakItUp(const G4Fragment &theFragment) 45 { 47 { >> 48 // G4FragmentVector * theResult = new G4FragmentVector; >> 49 46 if (theFragment.GetExcitationEnergy() <= 0.0 50 if (theFragment.GetExcitationEnergy() <= 0.0) { 47 return nullptr; << 51 //G4FragmentVector * theResult = new G4FragmentVector; >> 52 //theResult->push_back(new G4Fragment(theFragment)); >> 53 return 0; 48 } 54 } 49 55 >> 56 50 // Maximun average multiplicity: M_0 = 2.6 f 57 // Maximun average multiplicity: M_0 = 2.6 for A ~ 200 51 // and M_0 = 3.3 for A <= 110 58 // and M_0 = 3.3 for A <= 110 52 G4double MaxAverageMultiplicity = 59 G4double MaxAverageMultiplicity = 53 G4StatMFParameters::GetMaxAverageMultiplic << 60 G4StatMFParameters::GetMaxAverageMultiplicity(static_cast<G4int>(theFragment.GetA())); 54 61 55 62 56 // We'll use two kinds of ensembles 63 // We'll use two kinds of ensembles 57 G4StatMFMicroCanonical * theMicrocanonicalEn 64 G4StatMFMicroCanonical * theMicrocanonicalEnsemble = 0; 58 G4StatMFMacroCanonical * theMacrocanonicalEn 65 G4StatMFMacroCanonical * theMacrocanonicalEnsemble = 0; >> 66 59 67 60 //------------------------------------------ << 68 //------------------------------------------------------- 61 // Direct simulation part (Microcanonical en << 69 // Direct simulation part (Microcanonical ensemble) 62 //------------------------------------------ << 70 //------------------------------------------------------- 63 71 64 // Microcanonical ensemble initialization << 72 // Microcanonical ensemble initialization 65 theMicrocanonicalEnsemble = new G4StatMFMicr 73 theMicrocanonicalEnsemble = new G4StatMFMicroCanonical(theFragment); 66 74 67 G4int Iterations = 0; 75 G4int Iterations = 0; 68 G4int IterationsLimit = 100000; 76 G4int IterationsLimit = 100000; 69 G4double Temperature = 0.0; 77 G4double Temperature = 0.0; 70 78 71 G4bool FirstTime = true; 79 G4bool FirstTime = true; 72 G4StatMFChannel * theChannel = 0; 80 G4StatMFChannel * theChannel = 0; 73 81 74 G4bool ChannelOk; 82 G4bool ChannelOk; 75 do { // Try to de-excite as much as Iterati 83 do { // Try to de-excite as much as IterationLimit permits 76 do { 84 do { 77 85 78 G4double theMeanMult = theMicrocanonical 86 G4double theMeanMult = theMicrocanonicalEnsemble->GetMeanMultiplicity(); 79 if (theMeanMult <= MaxAverageMultiplicit 87 if (theMeanMult <= MaxAverageMultiplicity) { 80 // G4cout << "MICROCANONICAL" << G4endl; 88 // G4cout << "MICROCANONICAL" << G4endl; 81 // Choose fragments atomic numbers and charg 89 // Choose fragments atomic numbers and charges from direct simulation 82 theChannel = theMicrocanonicalEnsemble->Choo 90 theChannel = theMicrocanonicalEnsemble->ChooseAandZ(theFragment); 83 _theEnsemble = theMicrocanonicalEnsemble; 91 _theEnsemble = theMicrocanonicalEnsemble; 84 } else { 92 } else { 85 //------------------------------------------ 93 //----------------------------------------------------- 86 // Non direct simulation part (Macrocanonica 94 // Non direct simulation part (Macrocanonical Ensemble) 87 //------------------------------------------ 95 //----------------------------------------------------- 88 if (FirstTime) { 96 if (FirstTime) { 89 // Macrocanonical ensemble initialization 97 // Macrocanonical ensemble initialization 90 theMacrocanonicalEnsemble = new G4StatMFMa 98 theMacrocanonicalEnsemble = new G4StatMFMacroCanonical(theFragment); 91 _theEnsemble = theMacrocanonicalEnsemble; 99 _theEnsemble = theMacrocanonicalEnsemble; 92 FirstTime = false; 100 FirstTime = false; 93 } 101 } 94 // G4cout << "MACROCANONICAL" << G4endl; 102 // G4cout << "MACROCANONICAL" << G4endl; 95 // Select calculated fragment total multipli 103 // Select calculated fragment total multiplicity, 96 // fragment atomic numbers and fragment char 104 // fragment atomic numbers and fragment charges. 97 theChannel = theMacrocanonicalEnsemble->Choo 105 theChannel = theMacrocanonicalEnsemble->ChooseAandZ(theFragment); 98 } 106 } 99 107 100 ChannelOk = theChannel->CheckFragments() 108 ChannelOk = theChannel->CheckFragments(); 101 if (!ChannelOk) delete theChannel; 109 if (!ChannelOk) delete theChannel; 102 110 103 // Loop checking, 05-Aug-2015, Vladimir << 104 } while (!ChannelOk); 111 } while (!ChannelOk); 105 112 106 113 107 if (theChannel->GetMultiplicity() <= 1) { 114 if (theChannel->GetMultiplicity() <= 1) { 108 G4FragmentVector * theResult = new G4Fra 115 G4FragmentVector * theResult = new G4FragmentVector; 109 theResult->push_back(new G4Fragment(theF 116 theResult->push_back(new G4Fragment(theFragment)); 110 delete theMicrocanonicalEnsemble; 117 delete theMicrocanonicalEnsemble; 111 if (theMacrocanonicalEnsemble != 0) dele 118 if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble; 112 delete theChannel; 119 delete theChannel; 113 return theResult; 120 return theResult; 114 } 121 } 115 122 116 //-------------------------------------- 123 //-------------------------------------- 117 // Second part of simulation procedure. 124 // Second part of simulation procedure. 118 //-------------------------------------- 125 //-------------------------------------- 119 126 120 // Find temperature of breaking channel. 127 // Find temperature of breaking channel. 121 Temperature = _theEnsemble->GetMeanTempera 128 Temperature = _theEnsemble->GetMeanTemperature(); // Initial guess for Temperature 122 129 123 if (FindTemperatureOfBreakingChannel(theFr 130 if (FindTemperatureOfBreakingChannel(theFragment,theChannel,Temperature)) break; 124 131 125 // Do not forget to delete this unusable c 132 // Do not forget to delete this unusable channel, for which we failed to find the temperature, 126 // otherwise for very proton-reach nuclei 133 // otherwise for very proton-reach nuclei it would lead to memory leak due to large 127 // number of iterations. N.B. "theChannel" 134 // number of iterations. N.B. "theChannel" is created in G4StatMFMacroCanonical::ChooseZ() 128 135 129 // G4cout << " Iteration # " << Iterations 136 // G4cout << " Iteration # " << Iterations << " Mean Temperature = " << Temperature << G4endl; 130 137 131 delete theChannel; 138 delete theChannel; 132 139 133 // Loop checking, 05-Aug-2015, Vladimir Iv << 134 } while (Iterations++ < IterationsLimit ); 140 } while (Iterations++ < IterationsLimit ); >> 141 >> 142 135 143 136 // If Iterations >= IterationsLimit means th 144 // If Iterations >= IterationsLimit means that we couldn't solve for temperature 137 if (Iterations >= IterationsLimit) 145 if (Iterations >= IterationsLimit) 138 throw G4HadronicException(__FILE__, __LINE 146 throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::BreakItUp: Was not possible to solve for temperature of breaking channel"); 139 << 147 >> 148 140 G4FragmentVector * theResult = theChannel-> 149 G4FragmentVector * theResult = theChannel-> 141 GetFragments(theFragment.GetA_asInt(),theF 150 GetFragments(theFragment.GetA_asInt(),theFragment.GetZ_asInt(),Temperature); 142 << 151 >> 152 >> 153 143 // ~~~~~~ Energy conservation Patch !!!!!!!! 154 // ~~~~~~ Energy conservation Patch !!!!!!!!!!!!!!!!!!!!!! 144 // Original nucleus 4-momentum in CM system 155 // Original nucleus 4-momentum in CM system 145 G4LorentzVector InitialMomentum(theFragment. 156 G4LorentzVector InitialMomentum(theFragment.GetMomentum()); 146 InitialMomentum.boost(-InitialMomentum.boost 157 InitialMomentum.boost(-InitialMomentum.boostVector()); 147 G4double ScaleFactor = 0.0; 158 G4double ScaleFactor = 0.0; 148 G4double SavedScaleFactor = 0.0; 159 G4double SavedScaleFactor = 0.0; 149 do { 160 do { 150 G4double FragmentsEnergy = 0.0; 161 G4double FragmentsEnergy = 0.0; 151 for (auto const & ptr : *theResult) { << 162 G4FragmentVector::iterator j; 152 FragmentsEnergy += ptr->GetMomentum().e( << 163 for (j = theResult->begin(); j != theResult->end(); j++) 153 } << 164 FragmentsEnergy += (*j)->GetMomentum().e(); 154 if (0.0 == FragmentsEnergy) { break; } << 155 SavedScaleFactor = ScaleFactor; 165 SavedScaleFactor = ScaleFactor; 156 ScaleFactor = InitialMomentum.e()/Fragment 166 ScaleFactor = InitialMomentum.e()/FragmentsEnergy; 157 G4ThreeVector ScaledMomentum(0.0,0.0,0.0); 167 G4ThreeVector ScaledMomentum(0.0,0.0,0.0); 158 for (auto const & ptr : *theResult) { << 168 for (j = theResult->begin(); j != theResult->end(); j++) { 159 ScaledMomentum = ScaleFactor * ptr->GetM << 169 ScaledMomentum = ScaleFactor * (*j)->GetMomentum().vect(); 160 G4double Mass = ptr->GetMomentum().mag() << 170 G4double Mass = (*j)->GetMomentum().m(); 161 G4LorentzVector NewMomentum; 171 G4LorentzVector NewMomentum; 162 NewMomentum.setVect(ScaledMomentum); 172 NewMomentum.setVect(ScaledMomentum); 163 NewMomentum.setE(std::sqrt(ScaledMomentu 173 NewMomentum.setE(std::sqrt(ScaledMomentum.mag2()+Mass*Mass)); 164 ptr->SetMomentum(NewMomentum); << 174 (*j)->SetMomentum(NewMomentum); 165 } 175 } 166 // Loop checking, 05-Aug-2015, Vladimir Iv << 167 } while (ScaleFactor > 1.0+1.e-5 && std::abs 176 } while (ScaleFactor > 1.0+1.e-5 && std::abs(ScaleFactor-SavedScaleFactor)/ScaleFactor > 1.e-10); 168 // ~~~~~~ End of patch !!!!!!!!!!!!!!!!!!!!! 177 // ~~~~~~ End of patch !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 169 178 170 // Perform Lorentz boost 179 // Perform Lorentz boost 171 G4FragmentVector::iterator i; 180 G4FragmentVector::iterator i; 172 for (i = theResult->begin(); i != theResult- 181 for (i = theResult->begin(); i != theResult->end(); i++) { 173 G4LorentzVector FourMom = (*i)->GetMomentu 182 G4LorentzVector FourMom = (*i)->GetMomentum(); 174 FourMom.boost(theFragment.GetMomentum().bo 183 FourMom.boost(theFragment.GetMomentum().boostVector()); 175 (*i)->SetMomentum(FourMom); 184 (*i)->SetMomentum(FourMom); 176 (*i)->SetCreatorModelID(_secID); << 185 #ifdef PRECOMPOUND_TEST >> 186 (*i)->SetCreatorModel(G4String("G4StatMF")); >> 187 #endif 177 } 188 } 178 189 179 // garbage collection 190 // garbage collection 180 delete theMicrocanonicalEnsemble; 191 delete theMicrocanonicalEnsemble; 181 if (theMacrocanonicalEnsemble != 0) delete t 192 if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble; 182 delete theChannel; 193 delete theChannel; 183 194 184 return theResult; 195 return theResult; 185 } 196 } 186 197 187 198 188 G4bool G4StatMF::FindTemperatureOfBreakingChan 199 G4bool G4StatMF::FindTemperatureOfBreakingChannel(const G4Fragment & theFragment, 189 const G4StatMFChannel * aChannel 200 const G4StatMFChannel * aChannel, 190 G4double & Temperature) 201 G4double & Temperature) 191 // This finds temperature of breaking channe 202 // This finds temperature of breaking channel. 192 { 203 { 193 G4int A = theFragment.GetA_asInt(); 204 G4int A = theFragment.GetA_asInt(); 194 G4int Z = theFragment.GetZ_asInt(); 205 G4int Z = theFragment.GetZ_asInt(); 195 G4double U = theFragment.GetExcitationEnergy 206 G4double U = theFragment.GetExcitationEnergy(); 196 207 197 G4double T = std::max(Temperature,0.0012*MeV << 208 G4double T = std::max(Temperature,0.0012*MeV); >> 209 198 G4double Ta = T; 210 G4double Ta = T; >> 211 G4double Tb = T; >> 212 >> 213 199 G4double TotalEnergy = CalcEnergy(A,Z,aChann 214 G4double TotalEnergy = CalcEnergy(A,Z,aChannel,T); 200 215 201 G4double Da = (U - TotalEnergy)/U; 216 G4double Da = (U - TotalEnergy)/U; 202 G4double Db = 0.0; 217 G4double Db = 0.0; 203 218 204 // bracketing the solution 219 // bracketing the solution 205 if (Da == 0.0) { 220 if (Da == 0.0) { 206 Temperature = T; 221 Temperature = T; 207 return true; 222 return true; 208 } else if (Da < 0.0) { 223 } else if (Da < 0.0) { 209 do { 224 do { 210 T *= 0.5; << 225 Tb -= 0.5 * std::fabs(Tb); 211 if (T < 0.001*MeV) return false; << 226 T = Tb; >> 227 if (Tb < 0.001*MeV) return false; 212 228 213 TotalEnergy = CalcEnergy(A,Z,aChannel,T) 229 TotalEnergy = CalcEnergy(A,Z,aChannel,T); 214 230 215 Db = (U - TotalEnergy)/U; 231 Db = (U - TotalEnergy)/U; 216 // Loop checking, 05-Aug-2015, Vladimir << 217 } while (Db < 0.0); 232 } while (Db < 0.0); 218 233 219 } else { 234 } else { 220 do { 235 do { 221 T *= 1.5; << 236 Tb += 0.5 * std::fabs(Tb); >> 237 T = Tb; 222 238 223 TotalEnergy = CalcEnergy(A,Z,aChannel,T) 239 TotalEnergy = CalcEnergy(A,Z,aChannel,T); 224 240 225 Db = (U - TotalEnergy)/U; 241 Db = (U - TotalEnergy)/U; 226 // Loop checking, 05-Aug-2015, Vladimir << 227 } while (Db > 0.0); 242 } while (Db > 0.0); 228 } 243 } 229 244 230 G4double eps = 1.0e-14 * std::abs(T-Ta); << 245 G4double eps = 1.0e-14 * std::abs(Tb-Ta); 231 //G4double eps = 1.0e-3 ; 246 //G4double eps = 1.0e-3 ; 232 247 233 // Start the bisection method 248 // Start the bisection method 234 for (G4int j = 0; j < 1000; j++) { 249 for (G4int j = 0; j < 1000; j++) { 235 G4double Tc = (Ta+T)*0.5; << 250 G4double Tc = (Ta+Tb)/2.0; 236 if (std::abs(Ta-Tc) <= eps) { 251 if (std::abs(Ta-Tc) <= eps) { 237 Temperature = Tc; 252 Temperature = Tc; 238 return true; 253 return true; 239 } 254 } 240 255 241 T = Tc; << 256 T = Tc; >> 257 242 TotalEnergy = CalcEnergy(A,Z,aChannel,T); 258 TotalEnergy = CalcEnergy(A,Z,aChannel,T); >> 259 243 G4double Dc = (U - TotalEnergy)/U; 260 G4double Dc = (U - TotalEnergy)/U; 244 261 245 if (Dc == 0.0) { 262 if (Dc == 0.0) { 246 Temperature = Tc; 263 Temperature = Tc; 247 return true; 264 return true; 248 } 265 } >> 266 249 if (Da*Dc < 0.0) { 267 if (Da*Dc < 0.0) { 250 T = Tc; << 268 Tb = Tc; 251 Db = Dc; 269 Db = Dc; 252 } else { 270 } else { 253 Ta = Tc; 271 Ta = Tc; 254 Da = Dc; 272 Da = Dc; 255 } 273 } 256 } 274 } 257 275 258 Temperature = (Ta+T)*0.5; << 276 Temperature = (Ta+Tb)/2.0; 259 return false; 277 return false; 260 } 278 } 261 279 262 G4double G4StatMF::CalcEnergy(G4int A, G4int Z 280 G4double G4StatMF::CalcEnergy(G4int A, G4int Z, const G4StatMFChannel * aChannel, 263 G4double T) 281 G4double T) 264 { 282 { 265 G4double MassExcess0 = G4NucleiProperties::G << 283 G4double MassExcess0 = G4NucleiProperties::GetMassExcess(A,Z); 266 G4double ChannelEnergy = aChannel->GetFragme << 284 267 return -MassExcess0 + G4StatMFParameters::Ge << 285 G4double Coulomb = (3./5.)*(elm_coupling*Z*Z) >> 286 *std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1./3.)/ >> 287 (G4StatMFParameters::Getr0()*G4Pow::GetInstance()->Z13(A)); >> 288 >> 289 G4double ChannelEnergy = aChannel->GetFragmentsEnergy(T); >> 290 >> 291 return -MassExcess0 + Coulomb + ChannelEnergy; >> 292 268 } 293 } 269 294 270 295 271 296 272 297