[819] | 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 | // * * |
---|
| 21 | // * Parts of this code which have been developed by QinetiQ Ltd * |
---|
| 22 | // * under contract to the European Space Agency (ESA) are the * |
---|
| 23 | // * intellectual property of ESA. Rights to use, copy, modify and * |
---|
| 24 | // * redistribute this software for general public use are granted * |
---|
| 25 | // * in compliance with any licensing, distribution and development * |
---|
| 26 | // * policy adopted by the Geant4 Collaboration. This code has been * |
---|
| 27 | // * written by QinetiQ Ltd for the European Space Agency, under ESA * |
---|
| 28 | // * contract 17191/03/NL/LvH (Aurora Programme). * |
---|
| 29 | // * * |
---|
| 30 | // * By using, copying, modifying or distributing the software (or * |
---|
| 31 | // * any work based on the software) you agree to acknowledge its * |
---|
| 32 | // * use in resulting scientific publications, and indicate your * |
---|
| 33 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 34 | // ******************************************************************** |
---|
| 35 | // |
---|
| 36 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 37 | // |
---|
| 38 | // MODULE: G4WilsonAbrasionModel.cc |
---|
| 39 | // |
---|
[1228] | 40 | // Version: 1.0 |
---|
| 41 | // Date: 08/12/2009 |
---|
[819] | 42 | // Author: P R Truscott |
---|
| 43 | // Organisation: QinetiQ Ltd, UK |
---|
| 44 | // Customer: ESA/ESTEC, NOORDWIJK |
---|
| 45 | // Contract: 17191/03/NL/LvH |
---|
| 46 | // |
---|
| 47 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 48 | // |
---|
| 49 | // CHANGE HISTORY |
---|
| 50 | // -------------- |
---|
| 51 | // |
---|
| 52 | // 6 October 2003, P R Truscott, QinetiQ Ltd, UK |
---|
| 53 | // Created. |
---|
| 54 | // |
---|
| 55 | // 15 March 2004, P R Truscott, QinetiQ Ltd, UK |
---|
| 56 | // Beta release |
---|
| 57 | // |
---|
| 58 | // 18 January 2005, M H Mendenhall, Vanderbilt University, US |
---|
| 59 | // Pointers to theAbrasionGeometry and products generated by the deexcitation |
---|
| 60 | // handler deleted to prevent memory leaks. Also particle change of projectile |
---|
| 61 | // fragment previously not properly defined. |
---|
| 62 | // |
---|
[1228] | 63 | // 08 December 2009, P R Truscott, QinetiQ Ltd, Ltd |
---|
| 64 | // ver 1.0 |
---|
| 65 | // There was originally a possibility of the minimum impact parameter AFTER |
---|
| 66 | // considering Couloumb repulsion, rm, being too large. Now if: |
---|
| 67 | // rm >= fradius * (rP + rT) |
---|
| 68 | // where fradius is currently 0.99, then it is assumed the primary track is |
---|
| 69 | // unchanged. Additional conditions to escape from while-loop over impact |
---|
| 70 | // parameter: if the loop counter evtcnt reaches 1000, the primary track is |
---|
| 71 | // continued. |
---|
| 72 | // Additional clauses have been included in |
---|
| 73 | // G4WilsonAbrasionModel::GetNucleonInducedExcitation |
---|
| 74 | // Previously it was possible to get sqrt of negative number as Wilson algorithm |
---|
| 75 | // not properly defined if either: |
---|
| 76 | // rT > rP && rsq < rTsq - rPsq) or (rP > rT && rsq < rPsq - rTsq) |
---|
| 77 | // |
---|
[819] | 78 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 79 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 80 | // |
---|
| 81 | #include "G4WilsonAbrasionModel.hh" |
---|
| 82 | #include "G4WilsonRadius.hh" |
---|
| 83 | #include "G4NuclearAbrasionGeometry.hh" |
---|
| 84 | #include "G4WilsonAblationModel.hh" |
---|
| 85 | |
---|
| 86 | #include "G4ExcitationHandler.hh" |
---|
| 87 | #include "G4Evaporation.hh" |
---|
| 88 | #include "G4FermiBreakUp.hh" |
---|
| 89 | #include "G4StatMF.hh" |
---|
| 90 | #include "G4ParticleDefinition.hh" |
---|
| 91 | #include "G4DynamicParticle.hh" |
---|
| 92 | #include "Randomize.hh" |
---|
| 93 | #include "G4Fragment.hh" |
---|
| 94 | #include "G4VNuclearDensity.hh" |
---|
| 95 | #include "G4NuclearShellModelDensity.hh" |
---|
| 96 | #include "G4NuclearFermiDensity.hh" |
---|
| 97 | #include "G4FermiMomentum.hh" |
---|
| 98 | #include "G4ReactionProductVector.hh" |
---|
| 99 | #include "G4LorentzVector.hh" |
---|
| 100 | #include "G4ParticleMomentum.hh" |
---|
| 101 | #include "G4Poisson.hh" |
---|
| 102 | #include "G4ParticleTable.hh" |
---|
| 103 | #include "G4IonTable.hh" |
---|
| 104 | #include "globals.hh" |
---|
| 105 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 106 | // |
---|
| 107 | G4WilsonAbrasionModel::G4WilsonAbrasionModel (G4bool useAblation1) |
---|
| 108 | :G4HadronicInteraction("G4WilsonAbrasion") |
---|
| 109 | { |
---|
| 110 | // |
---|
| 111 | // |
---|
| 112 | // Send message to stdout to advise that the G4Abrasion model is being used. |
---|
| 113 | // |
---|
| 114 | PrintWelcomeMessage(); |
---|
| 115 | // |
---|
| 116 | // |
---|
| 117 | // Set the default verbose level to 0 - no output. |
---|
| 118 | // |
---|
| 119 | verboseLevel = 0; |
---|
| 120 | useAblation = useAblation1; |
---|
| 121 | // |
---|
| 122 | // |
---|
| 123 | // No de-excitation handler has been supplied - define the default handler. |
---|
| 124 | // |
---|
| 125 | theExcitationHandler = new G4ExcitationHandler; |
---|
| 126 | theExcitationHandlerx = new G4ExcitationHandler; |
---|
| 127 | if (useAblation) |
---|
| 128 | { |
---|
| 129 | theAblation = new G4WilsonAblationModel; |
---|
| 130 | theAblation->SetVerboseLevel(verboseLevel); |
---|
| 131 | theExcitationHandler->SetEvaporation(theAblation); |
---|
| 132 | theExcitationHandlerx->SetEvaporation(theAblation); |
---|
| 133 | } |
---|
| 134 | else |
---|
| 135 | { |
---|
| 136 | theAblation = NULL; |
---|
| 137 | G4Evaporation * theEvaporation = new G4Evaporation; |
---|
| 138 | G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp; |
---|
| 139 | G4StatMF * theMF = new G4StatMF; |
---|
| 140 | theExcitationHandler->SetEvaporation(theEvaporation); |
---|
| 141 | theExcitationHandler->SetFermiModel(theFermiBreakUp); |
---|
| 142 | theExcitationHandler->SetMultiFragmentation(theMF); |
---|
| 143 | theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6); |
---|
| 144 | theExcitationHandler->SetMinEForMultiFrag(5.0*MeV); |
---|
| 145 | |
---|
| 146 | theEvaporation = new G4Evaporation; |
---|
| 147 | theFermiBreakUp = new G4FermiBreakUp; |
---|
| 148 | theExcitationHandlerx->SetEvaporation(theEvaporation); |
---|
| 149 | theExcitationHandlerx->SetFermiModel(theFermiBreakUp); |
---|
| 150 | theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6); |
---|
| 151 | } |
---|
| 152 | // |
---|
| 153 | // |
---|
| 154 | // Set the minimum and maximum range for the model (despite nomanclature, this |
---|
| 155 | // is in energy per nucleon number). |
---|
| 156 | // |
---|
| 157 | SetMinEnergy(70.0*MeV); |
---|
| 158 | SetMaxEnergy(10.1*GeV); |
---|
| 159 | isBlocked = false; |
---|
| 160 | // |
---|
| 161 | // |
---|
| 162 | // npK, when mutiplied by the nuclear Fermi momentum, determines the range of |
---|
| 163 | // momentum over which the secondary nucleon momentum is sampled. |
---|
| 164 | // |
---|
| 165 | npK = 5.0; |
---|
| 166 | B = 10.0 * MeV; |
---|
| 167 | third = 1.0 / 3.0; |
---|
[1228] | 168 | fradius = 0.99; |
---|
[819] | 169 | conserveEnergy = false; |
---|
| 170 | conserveMomentum = true; |
---|
| 171 | } |
---|
| 172 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 173 | // |
---|
| 174 | G4WilsonAbrasionModel::G4WilsonAbrasionModel (G4ExcitationHandler *aExcitationHandler) |
---|
| 175 | { |
---|
| 176 | // |
---|
| 177 | // |
---|
| 178 | // Send message to stdout to advise that the G4Abrasion model is being used. |
---|
| 179 | // |
---|
| 180 | PrintWelcomeMessage(); |
---|
| 181 | // |
---|
| 182 | // |
---|
| 183 | // Set the default verbose level to 0 - no output. |
---|
| 184 | // |
---|
| 185 | verboseLevel = 0; |
---|
| 186 | // |
---|
| 187 | // |
---|
| 188 | // The user is able to provide the excitation handler as well as an argument |
---|
| 189 | // which is provided in this instantiation is used to determine |
---|
| 190 | // whether the spectators of the interaction are free following the abrasion. |
---|
| 191 | // |
---|
| 192 | theExcitationHandler = aExcitationHandler; |
---|
| 193 | theExcitationHandlerx = new G4ExcitationHandler; |
---|
| 194 | G4Evaporation * theEvaporation = new G4Evaporation; |
---|
| 195 | G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp; |
---|
| 196 | theExcitationHandlerx->SetEvaporation(theEvaporation); |
---|
| 197 | theExcitationHandlerx->SetFermiModel(theFermiBreakUp); |
---|
| 198 | theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6); |
---|
| 199 | // |
---|
| 200 | // |
---|
| 201 | // Set the minimum and maximum range for the model (despite nomanclature, this |
---|
| 202 | // is in energy per nucleon number). |
---|
| 203 | // |
---|
| 204 | SetMinEnergy(70.0*MeV); |
---|
| 205 | SetMaxEnergy(10.1*GeV); |
---|
| 206 | isBlocked = false; |
---|
| 207 | // |
---|
| 208 | // |
---|
| 209 | // npK, when mutiplied by the nuclear Fermi momentum, determines the range of |
---|
| 210 | // momentum over which the secondary nucleon momentum is sampled. |
---|
| 211 | // |
---|
| 212 | npK = 5.0; |
---|
| 213 | B = 10.0 * MeV; |
---|
| 214 | third = 1.0 / 3.0; |
---|
[1228] | 215 | fradius = 0.99; |
---|
[819] | 216 | conserveEnergy = false; |
---|
| 217 | conserveMomentum = true; |
---|
| 218 | } |
---|
| 219 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 220 | // |
---|
| 221 | G4WilsonAbrasionModel::~G4WilsonAbrasionModel () |
---|
| 222 | { |
---|
| 223 | // |
---|
| 224 | // |
---|
| 225 | // The destructor doesn't have to do a great deal! |
---|
| 226 | // |
---|
[1228] | 227 | if (theExcitationHandler) delete theExcitationHandler; |
---|
| 228 | if (theExcitationHandlerx) delete theExcitationHandlerx; |
---|
| 229 | if (useAblation) delete theAblation; |
---|
| 230 | // delete theExcitationHandler; |
---|
| 231 | // delete theExcitationHandlerx; |
---|
[819] | 232 | } |
---|
| 233 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 234 | // |
---|
| 235 | G4HadFinalState *G4WilsonAbrasionModel::ApplyYourself ( |
---|
| 236 | const G4HadProjectile &theTrack, G4Nucleus &theTarget) |
---|
| 237 | { |
---|
| 238 | // |
---|
| 239 | // |
---|
| 240 | // The secondaries will be returned in G4HadFinalState &theParticleChange - |
---|
| 241 | // initialise this. The original track will always be discontinued and |
---|
| 242 | // secondaries followed. |
---|
| 243 | // |
---|
| 244 | theParticleChange.Clear(); |
---|
| 245 | theParticleChange.SetStatusChange(stopAndKill); |
---|
| 246 | // |
---|
| 247 | // |
---|
| 248 | // Get relevant information about the projectile and target (A, Z, energy/nuc, |
---|
| 249 | // momentum, etc). |
---|
| 250 | // |
---|
| 251 | const G4ParticleDefinition *definitionP = theTrack.GetDefinition(); |
---|
| 252 | const G4double AP = definitionP->GetBaryonNumber(); |
---|
| 253 | const G4double ZP = definitionP->GetPDGCharge(); |
---|
| 254 | G4LorentzVector pP = theTrack.Get4Momentum(); |
---|
| 255 | G4double E = theTrack.GetKineticEnergy()/AP; |
---|
| 256 | G4double AT = theTarget.GetN(); |
---|
| 257 | G4double ZT = theTarget.GetZ(); |
---|
| 258 | G4double TotalEPre = theTrack.GetTotalEnergy() + |
---|
| 259 | theTarget.AtomicMass(AT, ZT) + theTarget.GetEnergyDeposit(); |
---|
| 260 | G4double TotalEPost = 0.0; |
---|
| 261 | // |
---|
| 262 | // |
---|
| 263 | // Determine the radii of the projectile and target nuclei. |
---|
| 264 | // |
---|
| 265 | G4WilsonRadius aR; |
---|
| 266 | G4double rP = aR.GetWilsonRadius(AP); |
---|
| 267 | G4double rT = aR.GetWilsonRadius(AT); |
---|
| 268 | G4double rPsq = rP * rP; |
---|
| 269 | G4double rTsq = rT * rT; |
---|
| 270 | if (verboseLevel >= 2) |
---|
| 271 | { |
---|
| 272 | G4cout <<"########################################" |
---|
| 273 | <<"########################################" |
---|
| 274 | <<G4endl; |
---|
| 275 | G4cout.precision(6); |
---|
| 276 | G4cout <<"IN G4WilsonAbrasionModel" <<G4endl; |
---|
| 277 | G4cout <<"Initial projectile A=" <<AP |
---|
| 278 | <<", Z=" <<ZP |
---|
| 279 | <<", radius = " <<rP/fermi <<" fm" |
---|
| 280 | <<G4endl; |
---|
| 281 | G4cout <<"Initial target A=" <<AT |
---|
| 282 | <<", Z=" <<ZT |
---|
| 283 | <<", radius = " <<rT/fermi <<" fm" |
---|
| 284 | <<G4endl; |
---|
| 285 | G4cout <<"Projectile momentum and Energy/nuc = " <<pP <<" ," <<E <<G4endl; |
---|
| 286 | } |
---|
| 287 | // |
---|
| 288 | // |
---|
| 289 | // The following variables are used to determine the impact parameter in the |
---|
| 290 | // near-field (i.e. taking into consideration the electrostatic repulsion). |
---|
| 291 | // |
---|
| 292 | G4double rm = ZP * ZT * elm_coupling / (E * AP); |
---|
| 293 | G4double r = 0.0; |
---|
| 294 | G4double rsq = 0.0; |
---|
| 295 | // |
---|
| 296 | // |
---|
| 297 | // Initialise some of the variables which wll be used to calculate the chord- |
---|
| 298 | // length for nucleons in the projectile and target, and hence calculate the |
---|
| 299 | // number of abraded nucleons and the excitation energy. |
---|
| 300 | // |
---|
| 301 | G4NuclearAbrasionGeometry *theAbrasionGeometry = NULL; |
---|
| 302 | G4double CT = 0.0; |
---|
| 303 | G4double F = 0.0; |
---|
| 304 | G4int Dabr = 0; |
---|
| 305 | // |
---|
| 306 | // |
---|
| 307 | // The following loop is performed until the number of nucleons which are |
---|
| 308 | // abraded by the process is >1, i.e. an interaction MUST occur. |
---|
| 309 | // |
---|
| 310 | while (Dabr == 0) |
---|
| 311 | { |
---|
| 312 | // Added by MHM 20050119 to fix leaking memory on second pass through this loop |
---|
| 313 | if (theAbrasionGeometry) |
---|
| 314 | { |
---|
| 315 | delete theAbrasionGeometry; |
---|
| 316 | theAbrasionGeometry = NULL; |
---|
| 317 | } |
---|
| 318 | // |
---|
| 319 | // |
---|
| 320 | // Sample the impact parameter. For the moment, this class takes account of |
---|
| 321 | // electrostatic effects on the impact parameter, but (like HZETRN AND NUCFRG2) |
---|
| 322 | // does not make any correction for the effects of nuclear-nuclear repulsion. |
---|
| 323 | // |
---|
| 324 | G4double rPT = rP + rT; |
---|
| 325 | G4double rPTsq = rPT * rPT; |
---|
[1228] | 326 | // |
---|
| 327 | // |
---|
| 328 | // This is a "catch" to make sure we don't go into an infinite loop because the |
---|
| 329 | // energy is too low to overcome nuclear repulsion. PRT 20091023. If the |
---|
| 330 | // value of rm < fradius * rPT then we're unlikely to sample a small enough |
---|
| 331 | // impact parameter (energy of incident particle is too low). Return primary |
---|
| 332 | // and don't complete nuclear interaction analysis. |
---|
| 333 | // |
---|
| 334 | if (rm >= fradius * rPT) { |
---|
| 335 | theParticleChange.SetStatusChange(isAlive); |
---|
| 336 | theParticleChange.SetEnergyChange(theTrack.GetKineticEnergy()); |
---|
| 337 | theParticleChange.SetMomentumChange(theTrack.Get4Momentum().vect().unit()); |
---|
| 338 | if (verboseLevel >= 2) { |
---|
| 339 | G4cout <<"Particle energy too low to overcome repulsion." <<G4endl; |
---|
| 340 | G4cout <<"Event rejected and original track maintained" <<G4endl; |
---|
| 341 | G4cout <<"########################################" |
---|
| 342 | <<"########################################" |
---|
| 343 | <<G4endl; |
---|
| 344 | } |
---|
| 345 | return &theParticleChange; |
---|
| 346 | } |
---|
| 347 | // |
---|
| 348 | // |
---|
| 349 | // Now sample impact parameter until the criterion is met that projectile |
---|
| 350 | // and target overlap, but repulsion is taken into consideration. |
---|
| 351 | // |
---|
| 352 | G4int evtcnt = 0; |
---|
[819] | 353 | r = 1.1 * rPT; |
---|
[1228] | 354 | while (r > rPT && ++evtcnt < 1000) |
---|
[819] | 355 | { |
---|
| 356 | G4double bsq = rPTsq * G4UniformRand(); |
---|
| 357 | r = (rm + std::sqrt(rm*rm + 4.0*bsq)) / 2.0; |
---|
| 358 | } |
---|
[1228] | 359 | // |
---|
| 360 | // |
---|
| 361 | // We've tried to sample this 1000 times, but failed. Assume nuclei do not |
---|
| 362 | // collide. |
---|
| 363 | // |
---|
| 364 | if (evtcnt >= 1000) { |
---|
| 365 | theParticleChange.SetStatusChange(isAlive); |
---|
| 366 | theParticleChange.SetEnergyChange(theTrack.GetKineticEnergy()); |
---|
| 367 | theParticleChange.SetMomentumChange(theTrack.Get4Momentum().vect().unit()); |
---|
| 368 | if (verboseLevel >= 2) { |
---|
| 369 | G4cout <<"Particle energy too low to overcome repulsion." <<G4endl; |
---|
| 370 | G4cout <<"Event rejected and original track maintained" <<G4endl; |
---|
| 371 | G4cout <<"########################################" |
---|
| 372 | <<"########################################" |
---|
| 373 | <<G4endl; |
---|
| 374 | } |
---|
| 375 | return &theParticleChange; |
---|
| 376 | } |
---|
| 377 | |
---|
| 378 | |
---|
[819] | 379 | rsq = r * r; |
---|
| 380 | // |
---|
| 381 | // |
---|
| 382 | // Now determine the chord-length through the target nucleus. |
---|
| 383 | // |
---|
| 384 | if (rT > rP) |
---|
| 385 | { |
---|
| 386 | G4double x = (rPsq + rsq - rTsq) / 2.0 / r; |
---|
| 387 | if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x); |
---|
| 388 | else CT = 2.0 * std::sqrt(rTsq - rsq); |
---|
| 389 | } |
---|
| 390 | else |
---|
| 391 | { |
---|
| 392 | G4double x = (rTsq + rsq - rPsq) / 2.0 / r; |
---|
| 393 | if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x); |
---|
| 394 | else CT = 2.0 * rT; |
---|
| 395 | } |
---|
| 396 | // |
---|
| 397 | // |
---|
| 398 | // Determine the number of abraded nucleons. Note that the mean number of |
---|
| 399 | // abraded nucleons is used to sample the Poisson distribution. The Poisson |
---|
| 400 | // distribution is sampled only ten times with the current impact parameter, |
---|
| 401 | // and if it fails after this to find a case for which the number of abraded |
---|
| 402 | // nucleons >1, the impact parameter is re-sampled. |
---|
| 403 | // |
---|
| 404 | theAbrasionGeometry = new G4NuclearAbrasionGeometry(AP,AT,r); |
---|
| 405 | F = theAbrasionGeometry->F(); |
---|
| 406 | G4double lambda = 16.6*fermi / std::pow(E/MeV,0.26); |
---|
| 407 | G4double Mabr = F * AP * (1.0 - std::exp(-CT/lambda)); |
---|
| 408 | G4long n = 0; |
---|
| 409 | for (G4int i = 0; i<10; i++) |
---|
| 410 | { |
---|
| 411 | n = G4Poisson(Mabr); |
---|
| 412 | if (n > 0) |
---|
| 413 | { |
---|
| 414 | if (n>AP) Dabr = (G4int) AP; |
---|
| 415 | else Dabr = (G4int) n; |
---|
| 416 | break; |
---|
| 417 | } |
---|
| 418 | } |
---|
| 419 | } |
---|
| 420 | if (verboseLevel >= 2) |
---|
| 421 | { |
---|
| 422 | G4cout <<G4endl; |
---|
| 423 | G4cout <<"Impact parameter = " <<r/fermi <<" fm" <<G4endl; |
---|
| 424 | G4cout <<"# Abraded nucleons = " <<Dabr <<G4endl; |
---|
| 425 | } |
---|
| 426 | // |
---|
| 427 | // |
---|
| 428 | // The number of abraded nucleons must be no greater than the number of |
---|
| 429 | // nucleons in either the projectile or the target. If AP - Dabr < 2 or |
---|
| 430 | // AT - Dabr < 2 then either we have only a nucleon left behind in the |
---|
| 431 | // projectile/target or we've tried to abrade too many nucleons - and Dabr |
---|
| 432 | // should be limited. |
---|
| 433 | // |
---|
| 434 | if (AP - (G4double) Dabr < 2.0) Dabr = (G4int) AP; |
---|
| 435 | if (AT - (G4double) Dabr < 2.0) Dabr = (G4int) AT; |
---|
| 436 | // |
---|
| 437 | // |
---|
| 438 | // Determine the abraded secondary nucleons from the projectile. *fragmentP |
---|
| 439 | // is a pointer to the prefragment from the projectile and nSecP is the number |
---|
| 440 | // of nucleons in theParticleChange which have been abraded. The total energy |
---|
| 441 | // from these is determined. |
---|
| 442 | // |
---|
| 443 | G4ThreeVector boost = pP.findBoostToCM(); |
---|
| 444 | G4Fragment *fragmentP = GetAbradedNucleons (Dabr, AP, ZP, rP); |
---|
| 445 | G4int nSecP = theParticleChange.GetNumberOfSecondaries(); |
---|
| 446 | G4int i = 0; |
---|
| 447 | for (i=0; i<nSecP; i++) |
---|
| 448 | { |
---|
| 449 | TotalEPost += theParticleChange.GetSecondary(i)-> |
---|
| 450 | GetParticle()->GetTotalEnergy(); |
---|
| 451 | } |
---|
| 452 | // |
---|
| 453 | // |
---|
| 454 | // Determine the number of spectators in the interaction region for the |
---|
| 455 | // projectile. |
---|
| 456 | // |
---|
| 457 | G4int DspcP = (G4int) (AP*F) - Dabr; |
---|
| 458 | if (DspcP <= 0) DspcP = 0; |
---|
| 459 | else if (DspcP > AP-Dabr) DspcP = ((G4int) AP) - Dabr; |
---|
| 460 | // |
---|
| 461 | // |
---|
| 462 | // Determine excitation energy associated with excess surface area of the |
---|
| 463 | // projectile (EsP) and the excitation due to scattering of nucleons which are |
---|
| 464 | // retained within the projectile (ExP). Add the total energy from the excited |
---|
| 465 | // nucleus to the total energy of the secondaries. |
---|
| 466 | // |
---|
| 467 | G4bool excitationAbsorbedByProjectile = false; |
---|
| 468 | if (fragmentP != NULL) |
---|
| 469 | { |
---|
| 470 | G4double EsP = theAbrasionGeometry->GetExcitationEnergyOfProjectile(); |
---|
| 471 | G4double ExP = 0.0; |
---|
| 472 | if (Dabr < AT) |
---|
| 473 | excitationAbsorbedByProjectile = G4UniformRand() < 0.5; |
---|
| 474 | if (excitationAbsorbedByProjectile) |
---|
| 475 | ExP = GetNucleonInducedExcitation(rP, rT, r); |
---|
| 476 | G4double xP = EsP + ExP; |
---|
| 477 | if (xP > B*(AP-Dabr)) xP = B*(AP-Dabr); |
---|
| 478 | G4LorentzVector lorentzVector = fragmentP->GetMomentum(); |
---|
| 479 | lorentzVector.setE(lorentzVector.e()+xP); |
---|
| 480 | fragmentP->SetMomentum(lorentzVector); |
---|
| 481 | TotalEPost += lorentzVector.e(); |
---|
| 482 | } |
---|
| 483 | G4double EMassP = TotalEPost; |
---|
| 484 | // |
---|
| 485 | // |
---|
| 486 | // Determine the abraded secondary nucleons from the target. Note that it's |
---|
| 487 | // assumed that the same number of nucleons are abraded from the target as for |
---|
| 488 | // the projectile, and obviously no boost is applied to the products. *fragmentT |
---|
| 489 | // is a pointer to the prefragment from the target and nSec is the total number |
---|
| 490 | // of nucleons in theParticleChange which have been abraded. The total energy |
---|
| 491 | // from these is determined. |
---|
| 492 | // |
---|
| 493 | G4Fragment *fragmentT = GetAbradedNucleons (Dabr, AT, ZT, rT); |
---|
| 494 | G4int nSec = theParticleChange.GetNumberOfSecondaries(); |
---|
| 495 | for (i=nSecP; i<nSec; i++) |
---|
| 496 | { |
---|
| 497 | TotalEPost += theParticleChange.GetSecondary(i)-> |
---|
| 498 | GetParticle()->GetTotalEnergy(); |
---|
| 499 | } |
---|
| 500 | // |
---|
| 501 | // |
---|
| 502 | // Determine the number of spectators in the interaction region for the |
---|
| 503 | // target. |
---|
| 504 | // |
---|
| 505 | G4int DspcT = (G4int) (AT*F) - Dabr; |
---|
| 506 | if (DspcT <= 0) DspcT = 0; |
---|
| 507 | else if (DspcT > AP-Dabr) DspcT = ((G4int) AT) - Dabr; |
---|
| 508 | // |
---|
| 509 | // |
---|
| 510 | // Determine excitation energy associated with excess surface area of the |
---|
| 511 | // target (EsT) and the excitation due to scattering of nucleons which are |
---|
| 512 | // retained within the target (ExT). Add the total energy from the excited |
---|
| 513 | // nucleus to the total energy of the secondaries. |
---|
| 514 | // |
---|
| 515 | if (fragmentT != NULL) |
---|
| 516 | { |
---|
| 517 | G4double EsT = theAbrasionGeometry->GetExcitationEnergyOfTarget(); |
---|
| 518 | G4double ExT = 0.0; |
---|
| 519 | if (!excitationAbsorbedByProjectile) |
---|
| 520 | ExT = GetNucleonInducedExcitation(rT, rP, r); |
---|
| 521 | G4double xT = EsT + ExT; |
---|
| 522 | if (xT > B*(AT-Dabr)) xT = B*(AT-Dabr); |
---|
| 523 | G4LorentzVector lorentzVector = fragmentT->GetMomentum(); |
---|
| 524 | lorentzVector.setE(lorentzVector.e()+xT); |
---|
| 525 | fragmentT->SetMomentum(lorentzVector); |
---|
| 526 | TotalEPost += lorentzVector.e(); |
---|
| 527 | } |
---|
| 528 | // |
---|
| 529 | // |
---|
| 530 | // Now determine the difference between the pre and post interaction |
---|
| 531 | // energy - this will be used to determine the Lorentz boost if conservation |
---|
| 532 | // of energy is to be imposed/attempted. |
---|
| 533 | // |
---|
| 534 | G4double deltaE = TotalEPre - TotalEPost; |
---|
| 535 | if (deltaE > 0.0 && conserveEnergy) |
---|
| 536 | { |
---|
| 537 | G4double beta = std::sqrt(1.0 - EMassP*EMassP/std::pow(deltaE+EMassP,2.0)); |
---|
| 538 | boost = boost / boost.mag() * beta; |
---|
| 539 | } |
---|
| 540 | // |
---|
| 541 | // |
---|
| 542 | // Now boost the secondaries from the projectile. |
---|
| 543 | // |
---|
| 544 | G4ThreeVector pBalance = pP.vect(); |
---|
| 545 | for (i=0; i<nSecP; i++) |
---|
| 546 | { |
---|
| 547 | G4DynamicParticle *dynamicP = theParticleChange.GetSecondary(i)-> |
---|
| 548 | GetParticle(); |
---|
| 549 | G4LorentzVector lorentzVector = dynamicP->Get4Momentum(); |
---|
| 550 | lorentzVector.boost(-boost); |
---|
| 551 | dynamicP->Set4Momentum(lorentzVector); |
---|
| 552 | pBalance -= lorentzVector.vect(); |
---|
| 553 | } |
---|
| 554 | // |
---|
| 555 | // |
---|
| 556 | // Set the boost for the projectile prefragment. This is now based on the |
---|
| 557 | // conservation of momentum. However, if the user selected momentum of the |
---|
| 558 | // prefragment is not to be conserved this simply boosted to the velocity of the |
---|
| 559 | // original projectile times the ratio of the unexcited to the excited mass |
---|
| 560 | // of the prefragment (the excitation increases the effective mass of the |
---|
| 561 | // prefragment, and therefore modifying the boost is an attempt to prevent |
---|
| 562 | // the momentum of the prefragment being excessive). |
---|
| 563 | // |
---|
| 564 | if (fragmentP != NULL) |
---|
| 565 | { |
---|
| 566 | G4LorentzVector lorentzVector = fragmentP->GetMomentum(); |
---|
| 567 | G4double m = lorentzVector.m(); |
---|
| 568 | if (conserveMomentum) |
---|
| 569 | fragmentP->SetMomentum |
---|
| 570 | (G4LorentzVector(pBalance,std::sqrt(pBalance.mag2()+m*m+1.0*eV*eV))); |
---|
| 571 | else |
---|
| 572 | { |
---|
| 573 | G4double mg = fragmentP->GetGroundStateMass(); |
---|
| 574 | fragmentP->SetMomentum(lorentzVector.boost(-boost * mg/m)); |
---|
| 575 | } |
---|
| 576 | } |
---|
| 577 | // |
---|
| 578 | // |
---|
| 579 | // Output information to user if verbose information requested. |
---|
| 580 | // |
---|
| 581 | if (verboseLevel >= 2) |
---|
| 582 | { |
---|
| 583 | G4cout <<G4endl; |
---|
| 584 | G4cout <<"-----------------------------------" <<G4endl; |
---|
| 585 | G4cout <<"Secondary nucleons from projectile:" <<G4endl; |
---|
| 586 | G4cout <<"-----------------------------------" <<G4endl; |
---|
| 587 | G4cout.precision(7); |
---|
| 588 | for (i=0; i<nSecP; i++) |
---|
| 589 | { |
---|
| 590 | G4cout <<"Particle # " <<i <<G4endl; |
---|
| 591 | theParticleChange.GetSecondary(i)->GetParticle()->DumpInfo(); |
---|
| 592 | G4DynamicParticle *dyn = theParticleChange.GetSecondary(i)->GetParticle(); |
---|
| 593 | G4cout <<"New nucleon (P) " <<dyn->GetDefinition()->GetParticleName() |
---|
| 594 | <<" : " <<dyn->Get4Momentum() |
---|
| 595 | <<G4endl; |
---|
| 596 | } |
---|
| 597 | G4cout <<"---------------------------" <<G4endl; |
---|
| 598 | G4cout <<"The projectile prefragment:" <<G4endl; |
---|
| 599 | G4cout <<"---------------------------" <<G4endl; |
---|
| 600 | if (fragmentP != NULL) |
---|
| 601 | G4cout <<*fragmentP <<G4endl; |
---|
| 602 | else |
---|
| 603 | G4cout <<"(No residual prefragment)" <<G4endl; |
---|
| 604 | G4cout <<G4endl; |
---|
| 605 | G4cout <<"-------------------------------" <<G4endl; |
---|
| 606 | G4cout <<"Secondary nucleons from target:" <<G4endl; |
---|
| 607 | G4cout <<"-------------------------------" <<G4endl; |
---|
| 608 | G4cout.precision(7); |
---|
| 609 | for (i=nSecP; i<nSec; i++) |
---|
| 610 | { |
---|
| 611 | G4cout <<"Particle # " <<i <<G4endl; |
---|
| 612 | theParticleChange.GetSecondary(i)->GetParticle()->DumpInfo(); |
---|
| 613 | G4DynamicParticle *dyn = theParticleChange.GetSecondary(i)->GetParticle(); |
---|
| 614 | G4cout <<"New nucleon (T) " <<dyn->GetDefinition()->GetParticleName() |
---|
| 615 | <<" : " <<dyn->Get4Momentum() |
---|
| 616 | <<G4endl; |
---|
| 617 | } |
---|
| 618 | G4cout <<"-----------------------" <<G4endl; |
---|
| 619 | G4cout <<"The target prefragment:" <<G4endl; |
---|
| 620 | G4cout <<"-----------------------" <<G4endl; |
---|
| 621 | if (fragmentT != NULL) |
---|
| 622 | G4cout <<*fragmentT <<G4endl; |
---|
| 623 | else |
---|
| 624 | G4cout <<"(No residual prefragment)" <<G4endl; |
---|
| 625 | } |
---|
| 626 | // |
---|
| 627 | // |
---|
| 628 | // Now we can decay the nuclear fragments if present. The secondaries are |
---|
| 629 | // collected and boosted as well. This is performed first for the projectile... |
---|
| 630 | // |
---|
| 631 | if (fragmentP !=NULL) |
---|
| 632 | { |
---|
| 633 | G4ReactionProductVector *products = NULL; |
---|
| 634 | if (fragmentP->GetZ() != fragmentP->GetA()) |
---|
| 635 | products = theExcitationHandler->BreakItUp(*fragmentP); |
---|
| 636 | else |
---|
| 637 | products = theExcitationHandlerx->BreakItUp(*fragmentP); |
---|
| 638 | delete fragmentP; |
---|
| 639 | fragmentP = NULL; |
---|
| 640 | |
---|
| 641 | G4ReactionProductVector::iterator iter; |
---|
| 642 | for (iter = products->begin(); iter != products->end(); ++iter) |
---|
| 643 | { |
---|
| 644 | G4DynamicParticle *secondary = |
---|
| 645 | new G4DynamicParticle((*iter)->GetDefinition(), |
---|
| 646 | (*iter)->GetTotalEnergy(), (*iter)->GetMomentum()); |
---|
| 647 | theParticleChange.AddSecondary (secondary); // Added MHM 20050118 |
---|
| 648 | G4String particleName = (*iter)->GetDefinition()->GetParticleName(); |
---|
| 649 | delete (*iter); // get rid of leftover particle def! // Added MHM 20050118 |
---|
| 650 | if (verboseLevel >= 2 && particleName.find("[",0) < particleName.size()) |
---|
| 651 | { |
---|
| 652 | G4cout <<"------------------------" <<G4endl; |
---|
| 653 | G4cout <<"The projectile fragment:" <<G4endl; |
---|
| 654 | G4cout <<"------------------------" <<G4endl; |
---|
| 655 | G4cout <<" fragmentP = " <<particleName |
---|
| 656 | <<" Energy = " <<secondary->GetKineticEnergy() |
---|
| 657 | <<G4endl; |
---|
| 658 | } |
---|
| 659 | } |
---|
| 660 | delete products; // Added MHM 20050118 |
---|
| 661 | } |
---|
| 662 | // |
---|
| 663 | // |
---|
| 664 | // Now decay the target nucleus - no boost is applied since in this |
---|
| 665 | // approximation it is assumed that there is negligible momentum transfer from |
---|
| 666 | // the projectile. |
---|
| 667 | // |
---|
| 668 | if (fragmentT != NULL) |
---|
| 669 | { |
---|
| 670 | G4ReactionProductVector *products = NULL; |
---|
| 671 | if (fragmentT->GetZ() != fragmentT->GetA()) |
---|
| 672 | products = theExcitationHandler->BreakItUp(*fragmentT); |
---|
| 673 | else |
---|
| 674 | products = theExcitationHandlerx->BreakItUp(*fragmentT); |
---|
| 675 | delete fragmentT; |
---|
| 676 | fragmentT = NULL; |
---|
| 677 | |
---|
| 678 | G4ReactionProductVector::iterator iter; |
---|
| 679 | for (iter = products->begin(); iter != products->end(); ++iter) |
---|
| 680 | { |
---|
| 681 | G4DynamicParticle *secondary = |
---|
| 682 | new G4DynamicParticle((*iter)->GetDefinition(), |
---|
| 683 | (*iter)->GetTotalEnergy(), (*iter)->GetMomentum()); |
---|
| 684 | theParticleChange.AddSecondary (secondary); |
---|
| 685 | G4String particleName = (*iter)->GetDefinition()->GetParticleName(); |
---|
| 686 | delete (*iter); // get rid of leftover particle def! // Added MHM 20050118 |
---|
| 687 | if (verboseLevel >= 2 && particleName.find("[",0) < particleName.size()) |
---|
| 688 | { |
---|
| 689 | G4cout <<"--------------------" <<G4endl; |
---|
| 690 | G4cout <<"The target fragment:" <<G4endl; |
---|
| 691 | G4cout <<"--------------------" <<G4endl; |
---|
| 692 | G4cout <<" fragmentT = " <<particleName |
---|
| 693 | <<" Energy = " <<secondary->GetKineticEnergy() |
---|
| 694 | <<G4endl; |
---|
| 695 | } |
---|
| 696 | } |
---|
| 697 | delete products; // Added MHM 20050118 |
---|
| 698 | } |
---|
| 699 | |
---|
| 700 | if (verboseLevel >= 2) |
---|
| 701 | G4cout <<"########################################" |
---|
| 702 | <<"########################################" |
---|
| 703 | <<G4endl; |
---|
| 704 | |
---|
| 705 | delete theAbrasionGeometry; |
---|
| 706 | |
---|
| 707 | return &theParticleChange; |
---|
| 708 | } |
---|
| 709 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 710 | // |
---|
| 711 | G4Fragment *G4WilsonAbrasionModel::GetAbradedNucleons (G4int Dabr, G4double A, |
---|
| 712 | G4double Z, G4double r) |
---|
| 713 | { |
---|
| 714 | // |
---|
| 715 | // |
---|
| 716 | // Initialise variables. tau is the Fermi radius of the nucleus. The variables |
---|
| 717 | // p..., C... and g(amma) are used to help sample the secondary nucleon |
---|
| 718 | // spectrum. |
---|
| 719 | // |
---|
| 720 | |
---|
| 721 | G4double pK = hbarc * std::pow(9.0 * pi / 4.0 * A, third) / (1.29 * r); |
---|
| 722 | if (A <= 24.0) pK *= -0.229*std::pow(A,third) + 1.62; |
---|
| 723 | G4double pKsq = pK * pK; |
---|
| 724 | G4double p1sq = 2.0/5.0 * pKsq; |
---|
| 725 | G4double p2sq = 6.0/5.0 * pKsq; |
---|
| 726 | G4double p3sq = 500.0 * 500.0; |
---|
| 727 | G4double C1 = 1.0; |
---|
| 728 | G4double C2 = 0.03; |
---|
| 729 | G4double C3 = 0.0002; |
---|
| 730 | G4double g = 90.0 * MeV; |
---|
| 731 | G4double maxn = C1 + C2 + C3; |
---|
| 732 | // |
---|
| 733 | // |
---|
| 734 | // initialise the number of secondary nucleons abraded to zero, and initially set |
---|
| 735 | // the type of nucleon abraded to proton ... just for now. |
---|
| 736 | // |
---|
| 737 | G4double Aabr = 0.0; |
---|
| 738 | G4double Zabr = 0.0; |
---|
| 739 | G4ParticleDefinition *typeNucleon = G4Proton::ProtonDefinition(); |
---|
| 740 | G4DynamicParticle *dynamicNucleon = NULL; |
---|
| 741 | G4ParticleMomentum pabr(0.0, 0.0, 0.0); |
---|
| 742 | // |
---|
| 743 | // |
---|
| 744 | // Now go through each abraded nucleon and sample type, spectrum and angle. |
---|
| 745 | // |
---|
| 746 | for (G4int i=0; i<Dabr; i++) |
---|
| 747 | { |
---|
| 748 | // |
---|
| 749 | // |
---|
| 750 | // Sample the nucleon momentum distribution by simple rejection techniques. We |
---|
| 751 | // reject values of p == 0.0 since this causes bad behaviour in the sinh term. |
---|
| 752 | // |
---|
| 753 | G4double p = 0.0; |
---|
| 754 | G4bool found = false; |
---|
| 755 | while (!found) |
---|
| 756 | { |
---|
| 757 | while (p <= 0.0) p = npK * pK * G4UniformRand(); |
---|
| 758 | G4double psq = p * p; |
---|
| 759 | found = maxn * G4UniformRand() < C1*std::exp(-psq/p1sq/2.0) + |
---|
| 760 | C2*std::exp(-psq/p2sq/2.0) + C3*std::exp(-psq/p3sq/2.0) + p/g/std::sinh(p/g); |
---|
| 761 | } |
---|
| 762 | // |
---|
| 763 | // |
---|
| 764 | // Determine the type of particle abraded. Can only be proton or neutron, |
---|
| 765 | // and the probability is determine to be proportional to the ratio as found |
---|
| 766 | // in the nucleus at each stage. |
---|
| 767 | // |
---|
| 768 | G4double prob = (Z-Zabr)/(A-Aabr); |
---|
| 769 | if (G4UniformRand()<prob) |
---|
| 770 | { |
---|
| 771 | Zabr++; |
---|
| 772 | typeNucleon = G4Proton::ProtonDefinition(); |
---|
| 773 | } |
---|
| 774 | else |
---|
| 775 | typeNucleon = G4Neutron::NeutronDefinition(); |
---|
| 776 | Aabr++; |
---|
| 777 | // |
---|
| 778 | // |
---|
| 779 | // The angular distribution of the secondary nucleons is approximated to an |
---|
| 780 | // isotropic distribution in the rest frame of the nucleus (this will be Lorentz |
---|
| 781 | // boosted later. |
---|
| 782 | // |
---|
| 783 | G4double costheta = 2.*G4UniformRand()-1.0; |
---|
| 784 | G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta)); |
---|
| 785 | G4double phi = 2.0*pi*G4UniformRand()*rad; |
---|
| 786 | G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta); |
---|
| 787 | G4double nucleonMass = typeNucleon->GetPDGMass(); |
---|
| 788 | G4double E = std::sqrt(p*p + nucleonMass*nucleonMass)-nucleonMass; |
---|
| 789 | dynamicNucleon = new G4DynamicParticle(typeNucleon,direction,E); |
---|
| 790 | theParticleChange.AddSecondary (dynamicNucleon); |
---|
| 791 | pabr += p*direction; |
---|
| 792 | } |
---|
| 793 | // |
---|
| 794 | // |
---|
| 795 | // Next determine the details of the nuclear prefragment .. that is if there |
---|
| 796 | // is one or more protons in the residue. (Note that the 1 eV in the total |
---|
| 797 | // energy is a safety factor to avoid any possibility of negative rest mass |
---|
| 798 | // energy.) |
---|
| 799 | // |
---|
| 800 | G4Fragment *fragment = NULL; |
---|
| 801 | if (Z-Zabr>=1.0) |
---|
| 802 | { |
---|
| 803 | G4double ionMass = G4ParticleTable::GetParticleTable()->GetIonTable()-> |
---|
| 804 | GetIonMass(G4lrint(Z-Zabr),G4lrint(A-Aabr)); |
---|
| 805 | G4double E = std::sqrt(pabr.mag2() + ionMass*ionMass); |
---|
| 806 | G4LorentzVector lorentzVector = G4LorentzVector(-pabr, E + 1.0*eV); |
---|
| 807 | fragment = |
---|
| 808 | new G4Fragment((G4int) (A-Aabr), (G4int) (Z-Zabr), lorentzVector); |
---|
| 809 | } |
---|
| 810 | |
---|
| 811 | return fragment; |
---|
| 812 | } |
---|
| 813 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 814 | // |
---|
| 815 | G4double G4WilsonAbrasionModel::GetNucleonInducedExcitation |
---|
| 816 | (G4double rP, G4double rT, G4double r) |
---|
| 817 | { |
---|
| 818 | // |
---|
| 819 | // |
---|
| 820 | // Initialise variables. |
---|
| 821 | // |
---|
| 822 | G4double Cl = 0.0; |
---|
| 823 | G4double rPsq = rP * rP; |
---|
| 824 | G4double rTsq = rT * rT; |
---|
| 825 | G4double rsq = r * r; |
---|
| 826 | // |
---|
| 827 | // |
---|
| 828 | // Depending upon the impact parameter, a different form of the chord length is |
---|
| 829 | // is used. |
---|
| 830 | // |
---|
| 831 | if (r > rT) Cl = 2.0*std::sqrt(rPsq + 2.0*r*rT - rsq - rTsq); |
---|
| 832 | else Cl = 2.0*rP; |
---|
[1228] | 833 | // |
---|
| 834 | // |
---|
| 835 | // The next lines have been changed to include a "catch" to make sure if the |
---|
| 836 | // projectile and target are too close, Ct is set to twice rP or twice rT. |
---|
| 837 | // Otherwise the standard Wilson algorithm should work fine. |
---|
| 838 | // PRT 20091023. |
---|
| 839 | // |
---|
| 840 | G4double Ct = 0.0; |
---|
| 841 | if (rT > rP && rsq < rTsq - rPsq) Ct = 2.0 * rP; |
---|
| 842 | else if (rP > rT && rsq < rPsq - rTsq) Ct = 2.0 * rT; |
---|
| 843 | else { |
---|
| 844 | G4double bP = (rPsq+rsq-rTsq)/2.0/r; |
---|
| 845 | G4double x = rPsq - bP*bP; |
---|
| 846 | if (x < 0.0) { |
---|
| 847 | G4cerr <<"########################################" |
---|
| 848 | <<"########################################" |
---|
| 849 | <<G4endl; |
---|
| 850 | G4cerr <<"ERROR IN G4WilsonAbrasionModel::GetNucleonInducedExcitation" |
---|
| 851 | <<G4endl; |
---|
| 852 | G4cerr <<"rPsq - bP*bP < 0.0 and cannot be square-rooted" <<G4endl; |
---|
| 853 | G4cerr <<"Set to zero instead" <<G4endl; |
---|
| 854 | G4cerr <<"########################################" |
---|
| 855 | <<"########################################" |
---|
| 856 | <<G4endl; |
---|
| 857 | } |
---|
| 858 | Ct = 2.0*std::sqrt(x); |
---|
| 859 | } |
---|
[819] | 860 | |
---|
| 861 | G4double Ex = 13.0 * Cl / fermi; |
---|
| 862 | if (Ct > 1.5*fermi) |
---|
| 863 | Ex += 13.0 * Cl / fermi /3.0 * (Ct/fermi - 1.5); |
---|
| 864 | |
---|
| 865 | return Ex; |
---|
| 866 | } |
---|
| 867 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 868 | // |
---|
| 869 | void G4WilsonAbrasionModel::SetUseAblation (G4bool useAblation1) |
---|
| 870 | { |
---|
| 871 | if (useAblation != useAblation1) |
---|
| 872 | { |
---|
| 873 | useAblation = useAblation1; |
---|
| 874 | delete theExcitationHandler; |
---|
| 875 | delete theExcitationHandlerx; |
---|
| 876 | theExcitationHandler = new G4ExcitationHandler; |
---|
| 877 | theExcitationHandlerx = new G4ExcitationHandler; |
---|
| 878 | if (useAblation) |
---|
| 879 | { |
---|
| 880 | theAblation = new G4WilsonAblationModel; |
---|
| 881 | theAblation->SetVerboseLevel(verboseLevel); |
---|
| 882 | theExcitationHandler->SetEvaporation(theAblation); |
---|
| 883 | theExcitationHandlerx->SetEvaporation(theAblation); |
---|
| 884 | } |
---|
| 885 | else |
---|
| 886 | { |
---|
| 887 | theAblation = NULL; |
---|
| 888 | G4Evaporation * theEvaporation = new G4Evaporation; |
---|
| 889 | G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp; |
---|
| 890 | G4StatMF * theMF = new G4StatMF; |
---|
| 891 | theExcitationHandler->SetEvaporation(theEvaporation); |
---|
| 892 | theExcitationHandler->SetFermiModel(theFermiBreakUp); |
---|
| 893 | theExcitationHandler->SetMultiFragmentation(theMF); |
---|
| 894 | theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6); |
---|
| 895 | theExcitationHandler->SetMinEForMultiFrag(5.0*MeV); |
---|
| 896 | |
---|
| 897 | theEvaporation = new G4Evaporation; |
---|
| 898 | theFermiBreakUp = new G4FermiBreakUp; |
---|
| 899 | theExcitationHandlerx->SetEvaporation(theEvaporation); |
---|
| 900 | theExcitationHandlerx->SetFermiModel(theFermiBreakUp); |
---|
| 901 | theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6); |
---|
| 902 | } |
---|
| 903 | } |
---|
| 904 | return; |
---|
| 905 | } |
---|
| 906 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 907 | // |
---|
| 908 | void G4WilsonAbrasionModel::PrintWelcomeMessage () |
---|
| 909 | { |
---|
| 910 | G4cout <<G4endl; |
---|
| 911 | G4cout <<" *****************************************************************" |
---|
| 912 | <<G4endl; |
---|
| 913 | G4cout <<" Nuclear abrasion model for nuclear-nuclear interactions activated" |
---|
| 914 | <<G4endl; |
---|
| 915 | G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)" |
---|
| 916 | <<G4endl; |
---|
| 917 | G4cout <<" *****************************************************************" |
---|
| 918 | <<G4endl; |
---|
| 919 | G4cout << G4endl; |
---|
| 920 | |
---|
| 921 | return; |
---|
| 922 | } |
---|
| 923 | //////////////////////////////////////////////////////////////////////////////// |
---|
| 924 | // |
---|