Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // 23 // >> 24 // $Id: G4StatMF.cc,v 1.2 2003/11/03 17:53:05 hpw Exp $ >> 25 // GEANT4 tag $Name: geant4-06-00-patch-01 $ 27 // 26 // 28 // Hadronic Process: Nuclear De-excitations 27 // Hadronic Process: Nuclear De-excitations 29 // by V. Lara 28 // by V. Lara 30 29 31 #include "G4StatMF.hh" 30 #include "G4StatMF.hh" 32 #include "G4PhysicalConstants.hh" << 33 #include "G4SystemOfUnits.hh" << 34 #include "G4Pow.hh" << 35 #include "G4PhysicsModelCatalog.hh" << 36 31 37 G4StatMF::G4StatMF() << 32 >> 33 >> 34 // Default constructor >> 35 G4StatMF::G4StatMF() : _theEnsemble(0) {} >> 36 >> 37 >> 38 // Destructor >> 39 G4StatMF::~G4StatMF() {} //{if (_theEnsemble != 0) delete _theEnsemble;} >> 40 >> 41 >> 42 // Copy constructor >> 43 G4StatMF::G4StatMF(const G4StatMF & ) : G4VMultiFragmentation() 38 { 44 { 39 _secID = G4PhysicsModelCatalog::GetModelID(" << 45 throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::copy_constructor meant to not be accessable"); 40 } 46 } 41 47 42 G4StatMF::~G4StatMF() {} << 43 48 44 G4FragmentVector* G4StatMF::BreakItUp(const G4 << 49 // Operators >> 50 >> 51 G4StatMF & G4StatMF::operator=(const G4StatMF & ) >> 52 { >> 53 throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::operator= meant to not be accessable"); >> 54 return *this; >> 55 } >> 56 >> 57 >> 58 G4bool G4StatMF::operator==(const G4StatMF & ) >> 59 { >> 60 throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::operator== meant to not be accessable"); >> 61 return false; >> 62 } >> 63 >> 64 >> 65 G4bool G4StatMF::operator!=(const G4StatMF & ) 45 { 66 { >> 67 throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::operator!= meant to not be accessable"); >> 68 return true; >> 69 } >> 70 >> 71 >> 72 >> 73 >> 74 >> 75 G4FragmentVector * G4StatMF::BreakItUp(const G4Fragment &theFragment) >> 76 { >> 77 // G4FragmentVector * theResult = new G4FragmentVector; >> 78 46 if (theFragment.GetExcitationEnergy() <= 0.0 79 if (theFragment.GetExcitationEnergy() <= 0.0) { 47 return nullptr; << 80 G4FragmentVector * theResult = new G4FragmentVector; >> 81 theResult->push_back(new G4Fragment(theFragment)); >> 82 return 0; 48 } 83 } 49 84 >> 85 50 // Maximun average multiplicity: M_0 = 2.6 f 86 // Maximun average multiplicity: M_0 = 2.6 for A ~ 200 51 // and M_0 = 3.3 for A <= 110 87 // and M_0 = 3.3 for A <= 110 52 G4double MaxAverageMultiplicity = 88 G4double MaxAverageMultiplicity = 53 G4StatMFParameters::GetMaxAverageMultiplic << 89 G4StatMFParameters::GetMaxAverageMultiplicity(static_cast<G4int>(theFragment.GetA())); 54 90 55 91 56 // We'll use two kinds of ensembles 92 // We'll use two kinds of ensembles 57 G4StatMFMicroCanonical * theMicrocanonicalEn 93 G4StatMFMicroCanonical * theMicrocanonicalEnsemble = 0; 58 G4StatMFMacroCanonical * theMacrocanonicalEn 94 G4StatMFMacroCanonical * theMacrocanonicalEnsemble = 0; >> 95 59 96 60 //------------------------------------------ << 97 //------------------------------------------------------- 61 // Direct simulation part (Microcanonical en << 98 // Direct simulation part (Microcanonical ensemble) 62 //------------------------------------------ << 99 //------------------------------------------------------- 63 100 64 // Microcanonical ensemble initialization << 101 // Microcanonical ensemble initialization 65 theMicrocanonicalEnsemble = new G4StatMFMicr 102 theMicrocanonicalEnsemble = new G4StatMFMicroCanonical(theFragment); 66 103 67 G4int Iterations = 0; 104 G4int Iterations = 0; 68 G4int IterationsLimit = 100000; << 69 G4double Temperature = 0.0; 105 G4double Temperature = 0.0; 70 106 71 G4bool FirstTime = true; 107 G4bool FirstTime = true; 72 G4StatMFChannel * theChannel = 0; 108 G4StatMFChannel * theChannel = 0; 73 << 109 74 G4bool ChannelOk; 110 G4bool ChannelOk; 75 do { // Try to de-excite as much as Iterati << 111 do { // Try to de-excite as much as 10 times 76 do { 112 do { 77 113 78 G4double theMeanMult = theMicrocanonical 114 G4double theMeanMult = theMicrocanonicalEnsemble->GetMeanMultiplicity(); 79 if (theMeanMult <= MaxAverageMultiplicit 115 if (theMeanMult <= MaxAverageMultiplicity) { 80 // G4cout << "MICROCANONICAL" << G4endl; 116 // G4cout << "MICROCANONICAL" << G4endl; 81 // Choose fragments atomic numbers and charg 117 // Choose fragments atomic numbers and charges from direct simulation 82 theChannel = theMicrocanonicalEnsemble->Choo 118 theChannel = theMicrocanonicalEnsemble->ChooseAandZ(theFragment); 83 _theEnsemble = theMicrocanonicalEnsemble; 119 _theEnsemble = theMicrocanonicalEnsemble; 84 } else { 120 } else { 85 //------------------------------------------ 121 //----------------------------------------------------- 86 // Non direct simulation part (Macrocanonica 122 // Non direct simulation part (Macrocanonical Ensemble) 87 //------------------------------------------ 123 //----------------------------------------------------- 88 if (FirstTime) { 124 if (FirstTime) { 89 // Macrocanonical ensemble initialization 125 // Macrocanonical ensemble initialization 90 theMacrocanonicalEnsemble = new G4StatMFMa 126 theMacrocanonicalEnsemble = new G4StatMFMacroCanonical(theFragment); 91 _theEnsemble = theMacrocanonicalEnsemble; 127 _theEnsemble = theMacrocanonicalEnsemble; 92 FirstTime = false; 128 FirstTime = false; 93 } 129 } 94 // G4cout << "MACROCANONICAL" << G4endl; 130 // G4cout << "MACROCANONICAL" << G4endl; 95 // Select calculated fragment total multipli 131 // Select calculated fragment total multiplicity, 96 // fragment atomic numbers and fragment char 132 // fragment atomic numbers and fragment charges. 97 theChannel = theMacrocanonicalEnsemble->Choo 133 theChannel = theMacrocanonicalEnsemble->ChooseAandZ(theFragment); 98 } 134 } 99 135 100 ChannelOk = theChannel->CheckFragments() << 136 if (!(ChannelOk = theChannel->CheckFragments())) delete theChannel; 101 if (!ChannelOk) delete theChannel; << 102 137 103 // Loop checking, 05-Aug-2015, Vladimir << 104 } while (!ChannelOk); 138 } while (!ChannelOk); 105 139 106 140 107 if (theChannel->GetMultiplicity() <= 1) { 141 if (theChannel->GetMultiplicity() <= 1) { 108 G4FragmentVector * theResult = new G4Fra 142 G4FragmentVector * theResult = new G4FragmentVector; 109 theResult->push_back(new G4Fragment(theF 143 theResult->push_back(new G4Fragment(theFragment)); 110 delete theMicrocanonicalEnsemble; 144 delete theMicrocanonicalEnsemble; 111 if (theMacrocanonicalEnsemble != 0) dele 145 if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble; 112 delete theChannel; 146 delete theChannel; 113 return theResult; 147 return theResult; 114 } 148 } 115 149 116 //-------------------------------------- 150 //-------------------------------------- 117 // Second part of simulation procedure. 151 // Second part of simulation procedure. 118 //-------------------------------------- 152 //-------------------------------------- 119 153 120 // Find temperature of breaking channel. 154 // Find temperature of breaking channel. 121 Temperature = _theEnsemble->GetMeanTempera << 155 Temperature = _theEnsemble->GetMeanTemperature(); // Initial value for Temperature 122 << 156 123 if (FindTemperatureOfBreakingChannel(theFr 157 if (FindTemperatureOfBreakingChannel(theFragment,theChannel,Temperature)) break; 124 << 158 125 // Do not forget to delete this unusable c << 159 } while (Iterations++ < 10); 126 // otherwise for very proton-reach nuclei << 160 127 // number of iterations. N.B. "theChannel" << 161 128 << 162 // If Iterations >= 10 means that we couldn't solve for temperature 129 // G4cout << " Iteration # " << Iterations << 163 if (Iterations >= 10) 130 << 131 delete theChannel; << 132 << 133 // Loop checking, 05-Aug-2015, Vladimir Iv << 134 } while (Iterations++ < IterationsLimit ); << 135 << 136 // If Iterations >= IterationsLimit means th << 137 if (Iterations >= IterationsLimit) << 138 throw G4HadronicException(__FILE__, __LINE 164 throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::BreakItUp: Was not possible to solve for temperature of breaking channel"); 139 << 165 >> 166 140 G4FragmentVector * theResult = theChannel-> 167 G4FragmentVector * theResult = theChannel-> 141 GetFragments(theFragment.GetA_asInt(),theF << 168 GetFragments(theFragment.GetA(),theFragment.GetZ(),Temperature); 142 << 169 >> 170 >> 171 143 // ~~~~~~ Energy conservation Patch !!!!!!!! 172 // ~~~~~~ Energy conservation Patch !!!!!!!!!!!!!!!!!!!!!! 144 // Original nucleus 4-momentum in CM system 173 // Original nucleus 4-momentum in CM system 145 G4LorentzVector InitialMomentum(theFragment. 174 G4LorentzVector InitialMomentum(theFragment.GetMomentum()); 146 InitialMomentum.boost(-InitialMomentum.boost 175 InitialMomentum.boost(-InitialMomentum.boostVector()); 147 G4double ScaleFactor = 0.0; 176 G4double ScaleFactor = 0.0; 148 G4double SavedScaleFactor = 0.0; 177 G4double SavedScaleFactor = 0.0; 149 do { 178 do { 150 G4double FragmentsEnergy = 0.0; 179 G4double FragmentsEnergy = 0.0; 151 for (auto const & ptr : *theResult) { << 180 G4FragmentVector::iterator j; 152 FragmentsEnergy += ptr->GetMomentum().e( << 181 for (j = theResult->begin(); j != theResult->end(); j++) 153 } << 182 FragmentsEnergy += (*j)->GetMomentum().e(); 154 if (0.0 == FragmentsEnergy) { break; } << 155 SavedScaleFactor = ScaleFactor; 183 SavedScaleFactor = ScaleFactor; 156 ScaleFactor = InitialMomentum.e()/Fragment 184 ScaleFactor = InitialMomentum.e()/FragmentsEnergy; 157 G4ThreeVector ScaledMomentum(0.0,0.0,0.0); 185 G4ThreeVector ScaledMomentum(0.0,0.0,0.0); 158 for (auto const & ptr : *theResult) { << 186 for (j = theResult->begin(); j != theResult->end(); j++) { 159 ScaledMomentum = ScaleFactor * ptr->GetM << 187 ScaledMomentum = ScaleFactor * (*j)->GetMomentum().vect(); 160 G4double Mass = ptr->GetMomentum().mag() << 188 G4double Mass = (*j)->GetMomentum().m(); 161 G4LorentzVector NewMomentum; 189 G4LorentzVector NewMomentum; 162 NewMomentum.setVect(ScaledMomentum); 190 NewMomentum.setVect(ScaledMomentum); 163 NewMomentum.setE(std::sqrt(ScaledMomentu << 191 NewMomentum.setE(sqrt(ScaledMomentum.mag2()+Mass*Mass)); 164 ptr->SetMomentum(NewMomentum); << 192 (*j)->SetMomentum(NewMomentum); 165 } 193 } 166 // Loop checking, 05-Aug-2015, Vladimir Iv << 194 } while (ScaleFactor > 1.0+1.e-5 && abs(ScaleFactor-SavedScaleFactor)/ScaleFactor > 1.e-10); 167 } while (ScaleFactor > 1.0+1.e-5 && std::abs << 168 // ~~~~~~ End of patch !!!!!!!!!!!!!!!!!!!!! 195 // ~~~~~~ End of patch !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 169 196 170 // Perform Lorentz boost 197 // Perform Lorentz boost 171 G4FragmentVector::iterator i; 198 G4FragmentVector::iterator i; 172 for (i = theResult->begin(); i != theResult- 199 for (i = theResult->begin(); i != theResult->end(); i++) { 173 G4LorentzVector FourMom = (*i)->GetMomentu 200 G4LorentzVector FourMom = (*i)->GetMomentum(); 174 FourMom.boost(theFragment.GetMomentum().bo 201 FourMom.boost(theFragment.GetMomentum().boostVector()); 175 (*i)->SetMomentum(FourMom); 202 (*i)->SetMomentum(FourMom); 176 (*i)->SetCreatorModelID(_secID); << 203 #ifdef PRECOMPOUND_TEST >> 204 (*i)->SetCreatorModel(G4String("G4StatMF")); >> 205 #endif 177 } 206 } 178 207 179 // garbage collection 208 // garbage collection 180 delete theMicrocanonicalEnsemble; 209 delete theMicrocanonicalEnsemble; 181 if (theMacrocanonicalEnsemble != 0) delete t 210 if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble; 182 delete theChannel; 211 delete theChannel; 183 212 184 return theResult; 213 return theResult; 185 } 214 } 186 215 187 216 188 G4bool G4StatMF::FindTemperatureOfBreakingChan 217 G4bool G4StatMF::FindTemperatureOfBreakingChannel(const G4Fragment & theFragment, 189 const G4StatMFChannel * aChannel 218 const G4StatMFChannel * aChannel, 190 G4double & Temperature) 219 G4double & Temperature) 191 // This finds temperature of breaking channe 220 // This finds temperature of breaking channel. 192 { 221 { 193 G4int A = theFragment.GetA_asInt(); << 222 G4double A = theFragment.GetA(); 194 G4int Z = theFragment.GetZ_asInt(); << 223 G4double Z = theFragment.GetZ(); 195 G4double U = theFragment.GetExcitationEnergy 224 G4double U = theFragment.GetExcitationEnergy(); 196 225 197 G4double T = std::max(Temperature,0.0012*MeV << 226 G4double T = std::max(Temperature,0.0012*MeV); >> 227 198 G4double Ta = T; 228 G4double Ta = T; >> 229 G4double Tb = T; >> 230 >> 231 199 G4double TotalEnergy = CalcEnergy(A,Z,aChann 232 G4double TotalEnergy = CalcEnergy(A,Z,aChannel,T); 200 233 201 G4double Da = (U - TotalEnergy)/U; 234 G4double Da = (U - TotalEnergy)/U; 202 G4double Db = 0.0; 235 G4double Db = 0.0; 203 236 204 // bracketing the solution 237 // bracketing the solution 205 if (Da == 0.0) { 238 if (Da == 0.0) { 206 Temperature = T; 239 Temperature = T; 207 return true; 240 return true; 208 } else if (Da < 0.0) { 241 } else if (Da < 0.0) { 209 do { 242 do { 210 T *= 0.5; << 243 Tb -= 0.5 * abs(Tb); 211 if (T < 0.001*MeV) return false; << 244 T = Tb; >> 245 if (Tb < 0.001*MeV) return false; 212 246 213 TotalEnergy = CalcEnergy(A,Z,aChannel,T) 247 TotalEnergy = CalcEnergy(A,Z,aChannel,T); 214 248 215 Db = (U - TotalEnergy)/U; 249 Db = (U - TotalEnergy)/U; 216 // Loop checking, 05-Aug-2015, Vladimir << 217 } while (Db < 0.0); 250 } while (Db < 0.0); 218 251 219 } else { 252 } else { 220 do { 253 do { 221 T *= 1.5; << 254 Tb += 0.5 * abs(Tb); >> 255 T = Tb; 222 256 223 TotalEnergy = CalcEnergy(A,Z,aChannel,T) 257 TotalEnergy = CalcEnergy(A,Z,aChannel,T); 224 258 225 Db = (U - TotalEnergy)/U; 259 Db = (U - TotalEnergy)/U; 226 // Loop checking, 05-Aug-2015, Vladimir << 227 } while (Db > 0.0); 260 } while (Db > 0.0); 228 } 261 } 229 262 230 G4double eps = 1.0e-14 * std::abs(T-Ta); << 263 G4double eps = 1.0e-14 * abs(Tb-Ta); 231 //G4double eps = 1.0e-3 ; 264 //G4double eps = 1.0e-3 ; 232 265 233 // Start the bisection method 266 // Start the bisection method 234 for (G4int j = 0; j < 1000; j++) { 267 for (G4int j = 0; j < 1000; j++) { 235 G4double Tc = (Ta+T)*0.5; << 268 G4double Tc = (Ta+Tb)/2.0; 236 if (std::abs(Ta-Tc) <= eps) { << 269 if (abs(Ta-Tc) <= eps) { 237 Temperature = Tc; 270 Temperature = Tc; 238 return true; 271 return true; 239 } 272 } 240 273 241 T = Tc; << 274 T = Tc; >> 275 242 TotalEnergy = CalcEnergy(A,Z,aChannel,T); 276 TotalEnergy = CalcEnergy(A,Z,aChannel,T); >> 277 243 G4double Dc = (U - TotalEnergy)/U; 278 G4double Dc = (U - TotalEnergy)/U; 244 279 245 if (Dc == 0.0) { 280 if (Dc == 0.0) { 246 Temperature = Tc; 281 Temperature = Tc; 247 return true; 282 return true; 248 } 283 } >> 284 249 if (Da*Dc < 0.0) { 285 if (Da*Dc < 0.0) { 250 T = Tc; << 286 Tb = Tc; 251 Db = Dc; 287 Db = Dc; 252 } else { 288 } else { 253 Ta = Tc; 289 Ta = Tc; 254 Da = Dc; 290 Da = Dc; 255 } 291 } 256 } 292 } 257 293 258 Temperature = (Ta+T)*0.5; << 294 Temperature = (Ta+Tb)/2.0; 259 return false; 295 return false; 260 } 296 } 261 297 262 G4double G4StatMF::CalcEnergy(G4int A, G4int Z << 298 263 G4double T) << 299 >> 300 G4double G4StatMF::CalcEnergy(const G4double A, const G4double Z, const G4StatMFChannel * aChannel, >> 301 const G4double T) 264 { 302 { 265 G4double MassExcess0 = G4NucleiProperties::G << 303 G4double MassExcess0 = G4NucleiProperties::GetMassExcess(static_cast<G4int>(A),static_cast<G4int>(Z)); 266 G4double ChannelEnergy = aChannel->GetFragme << 304 267 return -MassExcess0 + G4StatMFParameters::Ge << 305 G4double Coulomb = (3./5.)*(elm_coupling*Z*Z)*pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1./3.)/ >> 306 (G4StatMFParameters::Getr0()*pow(A,1./3.)); >> 307 >> 308 G4double ChannelEnergy = aChannel->GetFragmentsEnergy(T); >> 309 >> 310 return -MassExcess0 + Coulomb + ChannelEnergy; >> 311 268 } 312 } 269 313 270 314 271 315 272 316