00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00083
00084 #include "G4WilsonAbrasionModel.hh"
00085 #include "G4WilsonRadius.hh"
00086 #include "G4NuclearAbrasionGeometry.hh"
00087 #include "G4WilsonAblationModel.hh"
00088
00089 #include "G4PhysicalConstants.hh"
00090 #include "G4SystemOfUnits.hh"
00091 #include "G4ExcitationHandler.hh"
00092 #include "G4Evaporation.hh"
00093 #include "G4FermiBreakUp.hh"
00094 #include "G4StatMF.hh"
00095 #include "G4ParticleDefinition.hh"
00096 #include "G4DynamicParticle.hh"
00097 #include "Randomize.hh"
00098 #include "G4Fragment.hh"
00099 #include "G4ReactionProductVector.hh"
00100 #include "G4LorentzVector.hh"
00101 #include "G4ParticleMomentum.hh"
00102 #include "G4Poisson.hh"
00103 #include "G4ParticleTable.hh"
00104 #include "G4IonTable.hh"
00105 #include "globals.hh"
00106
00107
00108 G4WilsonAbrasionModel::G4WilsonAbrasionModel(G4bool useAblation1)
00109 :G4HadronicInteraction("G4WilsonAbrasion")
00110 {
00111
00112 PrintWelcomeMessage();
00113
00114
00115 verboseLevel = 0;
00116 useAblation = useAblation1;
00117
00118
00119
00120 theExcitationHandler = new G4ExcitationHandler;
00121 theExcitationHandlerx = new G4ExcitationHandler;
00122 if (useAblation)
00123 {
00124 theAblation = new G4WilsonAblationModel;
00125 theAblation->SetVerboseLevel(verboseLevel);
00126 theExcitationHandler->SetEvaporation(theAblation);
00127 theExcitationHandlerx->SetEvaporation(theAblation);
00128 }
00129 else
00130 {
00131 theAblation = NULL;
00132 G4Evaporation * theEvaporation = new G4Evaporation;
00133 G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp;
00134 G4StatMF * theMF = new G4StatMF;
00135 theExcitationHandler->SetEvaporation(theEvaporation);
00136 theExcitationHandler->SetFermiModel(theFermiBreakUp);
00137 theExcitationHandler->SetMultiFragmentation(theMF);
00138 theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6);
00139 theExcitationHandler->SetMinEForMultiFrag(5.0*MeV);
00140
00141 theEvaporation = new G4Evaporation;
00142 theFermiBreakUp = new G4FermiBreakUp;
00143 theExcitationHandlerx->SetEvaporation(theEvaporation);
00144 theExcitationHandlerx->SetFermiModel(theFermiBreakUp);
00145 theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6);
00146 }
00147
00148
00149
00150
00151 SetMinEnergy(70.0*MeV);
00152 SetMaxEnergy(10.1*GeV);
00153 isBlocked = false;
00154
00155
00156
00157
00158 r0sq = 0.0;
00159 npK = 5.0;
00160 B = 10.0 * MeV;
00161 third = 1.0 / 3.0;
00162 fradius = 0.99;
00163 conserveEnergy = false;
00164 conserveMomentum = true;
00165 }
00166
00167 void G4WilsonAbrasionModel::ModelDescription(std::ostream& outFile) const
00168 {
00169 outFile << "G4WilsonAbrasionModel is a macroscopic treatment of\n"
00170 << "nucleus-nucleus collisions using simple geometric arguments.\n"
00171 << "The smaller projectile nucleus gouges out a part of the larger\n"
00172 << "target nucleus, leaving a residual nucleus and a fireball\n"
00173 << "region where the projectile and target intersect. The fireball"
00174 << "is then treated as a highly excited nuclear fragment. This\n"
00175 << "model is based on the NUCFRG2 model and is valid for all\n"
00176 << "projectile energies between 70 MeV/n and 10.1 GeV/n. \n";
00177 }
00178
00179 G4WilsonAbrasionModel::G4WilsonAbrasionModel(G4ExcitationHandler* aExcitationHandler)
00180 {
00181
00182
00183 PrintWelcomeMessage();
00184
00185
00186
00187 verboseLevel = 0;
00188
00189 theAblation = NULL;
00190 useAblation = false;
00191
00192
00193
00194
00195
00196
00197 theExcitationHandler = aExcitationHandler;
00198 theExcitationHandlerx = new G4ExcitationHandler;
00199 G4Evaporation * theEvaporation = new G4Evaporation;
00200 G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp;
00201 theExcitationHandlerx->SetEvaporation(theEvaporation);
00202 theExcitationHandlerx->SetFermiModel(theFermiBreakUp);
00203 theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6);
00204
00205
00206
00207
00208
00209 SetMinEnergy(70.0*MeV);
00210 SetMaxEnergy(10.1*GeV);
00211 isBlocked = false;
00212
00213
00214
00215
00216
00217 r0sq = 0.0;
00218 npK = 5.0;
00219 B = 10.0 * MeV;
00220 third = 1.0 / 3.0;
00221 fradius = 0.99;
00222 conserveEnergy = false;
00223 conserveMomentum = true;
00224 }
00226
00227 G4WilsonAbrasionModel::~G4WilsonAbrasionModel ()
00228 {
00229
00230
00231
00232
00233 if (theExcitationHandler) delete theExcitationHandler;
00234 if (theExcitationHandlerx) delete theExcitationHandlerx;
00235 if (useAblation) delete theAblation;
00236
00237
00238 }
00240
00241 G4HadFinalState *G4WilsonAbrasionModel::ApplyYourself (
00242 const G4HadProjectile &theTrack, G4Nucleus &theTarget)
00243 {
00244
00245
00246
00247
00248
00249
00250 theParticleChange.Clear();
00251 theParticleChange.SetStatusChange(stopAndKill);
00252
00253
00254
00255
00256
00257 const G4ParticleDefinition *definitionP = theTrack.GetDefinition();
00258 const G4double AP = definitionP->GetBaryonNumber();
00259 const G4double ZP = definitionP->GetPDGCharge();
00260 G4LorentzVector pP = theTrack.Get4Momentum();
00261 G4double E = theTrack.GetKineticEnergy()/AP;
00262 G4double AT = theTarget.GetA_asInt();
00263 G4double ZT = theTarget.GetZ_asInt();
00264 G4double TotalEPre = theTrack.GetTotalEnergy() +
00265 theTarget.AtomicMass(AT, ZT) + theTarget.GetEnergyDeposit();
00266 G4double TotalEPost = 0.0;
00267
00268
00269
00270
00271 G4WilsonRadius aR;
00272 G4double rP = aR.GetWilsonRadius(AP);
00273 G4double rT = aR.GetWilsonRadius(AT);
00274 G4double rPsq = rP * rP;
00275 G4double rTsq = rT * rT;
00276 if (verboseLevel >= 2)
00277 {
00278 G4cout <<"########################################"
00279 <<"########################################"
00280 <<G4endl;
00281 G4cout.precision(6);
00282 G4cout <<"IN G4WilsonAbrasionModel" <<G4endl;
00283 G4cout <<"Initial projectile A=" <<AP
00284 <<", Z=" <<ZP
00285 <<", radius = " <<rP/fermi <<" fm"
00286 <<G4endl;
00287 G4cout <<"Initial target A=" <<AT
00288 <<", Z=" <<ZT
00289 <<", radius = " <<rT/fermi <<" fm"
00290 <<G4endl;
00291 G4cout <<"Projectile momentum and Energy/nuc = " <<pP <<" ," <<E <<G4endl;
00292 }
00293
00294
00295
00296
00297
00298 G4double rm = ZP * ZT * elm_coupling / (E * AP);
00299 G4double r = 0.0;
00300 G4double rsq = 0.0;
00301
00302
00303
00304
00305
00306
00307 G4NuclearAbrasionGeometry *theAbrasionGeometry = NULL;
00308 G4double CT = 0.0;
00309 G4double F = 0.0;
00310 G4int Dabr = 0;
00311
00312
00313
00314
00315
00316 while (Dabr == 0)
00317 {
00318
00319 if (theAbrasionGeometry)
00320 {
00321 delete theAbrasionGeometry;
00322 theAbrasionGeometry = NULL;
00323 }
00324
00325
00326
00327
00328
00329
00330 G4double rPT = rP + rT;
00331 G4double rPTsq = rPT * rPT;
00332
00333
00334
00335
00336
00337
00338
00339
00340 if (rm >= fradius * rPT) {
00341 theParticleChange.SetStatusChange(isAlive);
00342 theParticleChange.SetEnergyChange(theTrack.GetKineticEnergy());
00343 theParticleChange.SetMomentumChange(theTrack.Get4Momentum().vect().unit());
00344 if (verboseLevel >= 2) {
00345 G4cout <<"Particle energy too low to overcome repulsion." <<G4endl;
00346 G4cout <<"Event rejected and original track maintained" <<G4endl;
00347 G4cout <<"########################################"
00348 <<"########################################"
00349 <<G4endl;
00350 }
00351 return &theParticleChange;
00352 }
00353
00354
00355
00356
00357
00358 G4int evtcnt = 0;
00359 r = 1.1 * rPT;
00360 while (r > rPT && ++evtcnt < 1000)
00361 {
00362 G4double bsq = rPTsq * G4UniformRand();
00363 r = (rm + std::sqrt(rm*rm + 4.0*bsq)) / 2.0;
00364 }
00365
00366
00367
00368
00369
00370 if (evtcnt >= 1000) {
00371 theParticleChange.SetStatusChange(isAlive);
00372 theParticleChange.SetEnergyChange(theTrack.GetKineticEnergy());
00373 theParticleChange.SetMomentumChange(theTrack.Get4Momentum().vect().unit());
00374 if (verboseLevel >= 2) {
00375 G4cout <<"Particle energy too low to overcome repulsion." <<G4endl;
00376 G4cout <<"Event rejected and original track maintained" <<G4endl;
00377 G4cout <<"########################################"
00378 <<"########################################"
00379 <<G4endl;
00380 }
00381 return &theParticleChange;
00382 }
00383
00384
00385 rsq = r * r;
00386
00387
00388
00389
00390 if (rT > rP)
00391 {
00392 G4double x = (rPsq + rsq - rTsq) / 2.0 / r;
00393 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
00394 else CT = 2.0 * std::sqrt(rTsq - rsq);
00395 }
00396 else
00397 {
00398 G4double x = (rTsq + rsq - rPsq) / 2.0 / r;
00399 if (x > 0.0) CT = 2.0 * std::sqrt(rTsq - x*x);
00400 else CT = 2.0 * rT;
00401 }
00402
00403
00404
00405
00406
00407
00408
00409
00410 theAbrasionGeometry = new G4NuclearAbrasionGeometry(AP,AT,r);
00411 F = theAbrasionGeometry->F();
00412 G4double lambda = 16.6*fermi / std::pow(E/MeV,0.26);
00413 G4double Mabr = F * AP * (1.0 - std::exp(-CT/lambda));
00414 G4long n = 0;
00415 for (G4int i = 0; i<10; i++)
00416 {
00417 n = G4Poisson(Mabr);
00418 if (n > 0)
00419 {
00420 if (n>AP) Dabr = (G4int) AP;
00421 else Dabr = (G4int) n;
00422 break;
00423 }
00424 }
00425 }
00426 if (verboseLevel >= 2)
00427 {
00428 G4cout <<G4endl;
00429 G4cout <<"Impact parameter = " <<r/fermi <<" fm" <<G4endl;
00430 G4cout <<"# Abraded nucleons = " <<Dabr <<G4endl;
00431 }
00432
00433
00434
00435
00436
00437
00438
00439
00440 if (AP - (G4double) Dabr < 2.0) Dabr = (G4int) AP;
00441 if (AT - (G4double) Dabr < 2.0) Dabr = (G4int) AT;
00442
00443
00444
00445
00446
00447
00448
00449 G4ThreeVector boost = pP.findBoostToCM();
00450 G4Fragment *fragmentP = GetAbradedNucleons (Dabr, AP, ZP, rP);
00451 G4int nSecP = theParticleChange.GetNumberOfSecondaries();
00452 G4int i = 0;
00453 for (i=0; i<nSecP; i++)
00454 {
00455 TotalEPost += theParticleChange.GetSecondary(i)->
00456 GetParticle()->GetTotalEnergy();
00457 }
00458
00459
00460
00461
00462
00463 G4int DspcP = (G4int) (AP*F) - Dabr;
00464 if (DspcP <= 0) DspcP = 0;
00465 else if (DspcP > AP-Dabr) DspcP = ((G4int) AP) - Dabr;
00466
00467
00468
00469
00470
00471
00472
00473 G4bool excitationAbsorbedByProjectile = false;
00474 if (fragmentP != NULL)
00475 {
00476 G4double EsP = theAbrasionGeometry->GetExcitationEnergyOfProjectile();
00477 G4double ExP = 0.0;
00478 if (Dabr < AT)
00479 excitationAbsorbedByProjectile = G4UniformRand() < 0.5;
00480 if (excitationAbsorbedByProjectile)
00481 ExP = GetNucleonInducedExcitation(rP, rT, r);
00482 G4double xP = EsP + ExP;
00483 if (xP > B*(AP-Dabr)) xP = B*(AP-Dabr);
00484 G4LorentzVector lorentzVector = fragmentP->GetMomentum();
00485 lorentzVector.setE(lorentzVector.e()+xP);
00486 fragmentP->SetMomentum(lorentzVector);
00487 TotalEPost += lorentzVector.e();
00488 }
00489 G4double EMassP = TotalEPost;
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499 G4Fragment *fragmentT = GetAbradedNucleons (Dabr, AT, ZT, rT);
00500 G4int nSec = theParticleChange.GetNumberOfSecondaries();
00501 for (i=nSecP; i<nSec; i++)
00502 {
00503 TotalEPost += theParticleChange.GetSecondary(i)->
00504 GetParticle()->GetTotalEnergy();
00505 }
00506
00507
00508
00509
00510
00511 G4int DspcT = (G4int) (AT*F) - Dabr;
00512 if (DspcT <= 0) DspcT = 0;
00513 else if (DspcT > AP-Dabr) DspcT = ((G4int) AT) - Dabr;
00514
00515
00516
00517
00518
00519
00520
00521 if (fragmentT != NULL)
00522 {
00523 G4double EsT = theAbrasionGeometry->GetExcitationEnergyOfTarget();
00524 G4double ExT = 0.0;
00525 if (!excitationAbsorbedByProjectile)
00526 ExT = GetNucleonInducedExcitation(rT, rP, r);
00527 G4double xT = EsT + ExT;
00528 if (xT > B*(AT-Dabr)) xT = B*(AT-Dabr);
00529 G4LorentzVector lorentzVector = fragmentT->GetMomentum();
00530 lorentzVector.setE(lorentzVector.e()+xT);
00531 fragmentT->SetMomentum(lorentzVector);
00532 TotalEPost += lorentzVector.e();
00533 }
00534
00535
00536
00537
00538
00539
00540 G4double deltaE = TotalEPre - TotalEPost;
00541 if (deltaE > 0.0 && conserveEnergy)
00542 {
00543 G4double beta = std::sqrt(1.0 - EMassP*EMassP/std::pow(deltaE+EMassP,2.0));
00544 boost = boost / boost.mag() * beta;
00545 }
00546
00547
00548
00549
00550 G4ThreeVector pBalance = pP.vect();
00551 for (i=0; i<nSecP; i++)
00552 {
00553 G4DynamicParticle *dynamicP = theParticleChange.GetSecondary(i)->
00554 GetParticle();
00555 G4LorentzVector lorentzVector = dynamicP->Get4Momentum();
00556 lorentzVector.boost(-boost);
00557 dynamicP->Set4Momentum(lorentzVector);
00558 pBalance -= lorentzVector.vect();
00559 }
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570 if (fragmentP != NULL)
00571 {
00572 G4LorentzVector lorentzVector = fragmentP->GetMomentum();
00573 G4double fragmentM = lorentzVector.m();
00574 if (conserveMomentum)
00575 fragmentP->SetMomentum
00576 (G4LorentzVector(pBalance,std::sqrt(pBalance.mag2()+fragmentM*fragmentM+1.0*eV*eV)));
00577 else
00578 {
00579 G4double fragmentGroundStateM = fragmentP->GetGroundStateMass();
00580 fragmentP->SetMomentum(lorentzVector.boost(-boost * fragmentGroundStateM/fragmentM));
00581 }
00582 }
00583
00584
00585
00586
00587 if (verboseLevel >= 2)
00588 {
00589 G4cout <<G4endl;
00590 G4cout <<"-----------------------------------" <<G4endl;
00591 G4cout <<"Secondary nucleons from projectile:" <<G4endl;
00592 G4cout <<"-----------------------------------" <<G4endl;
00593 G4cout.precision(7);
00594 for (i=0; i<nSecP; i++)
00595 {
00596 G4cout <<"Particle # " <<i <<G4endl;
00597 theParticleChange.GetSecondary(i)->GetParticle()->DumpInfo();
00598 G4DynamicParticle *dyn = theParticleChange.GetSecondary(i)->GetParticle();
00599 G4cout <<"New nucleon (P) " <<dyn->GetDefinition()->GetParticleName()
00600 <<" : " <<dyn->Get4Momentum()
00601 <<G4endl;
00602 }
00603 G4cout <<"---------------------------" <<G4endl;
00604 G4cout <<"The projectile prefragment:" <<G4endl;
00605 G4cout <<"---------------------------" <<G4endl;
00606 if (fragmentP != NULL)
00607 G4cout <<*fragmentP <<G4endl;
00608 else
00609 G4cout <<"(No residual prefragment)" <<G4endl;
00610 G4cout <<G4endl;
00611 G4cout <<"-------------------------------" <<G4endl;
00612 G4cout <<"Secondary nucleons from target:" <<G4endl;
00613 G4cout <<"-------------------------------" <<G4endl;
00614 G4cout.precision(7);
00615 for (i=nSecP; i<nSec; i++)
00616 {
00617 G4cout <<"Particle # " <<i <<G4endl;
00618 theParticleChange.GetSecondary(i)->GetParticle()->DumpInfo();
00619 G4DynamicParticle *dyn = theParticleChange.GetSecondary(i)->GetParticle();
00620 G4cout <<"New nucleon (T) " <<dyn->GetDefinition()->GetParticleName()
00621 <<" : " <<dyn->Get4Momentum()
00622 <<G4endl;
00623 }
00624 G4cout <<"-----------------------" <<G4endl;
00625 G4cout <<"The target prefragment:" <<G4endl;
00626 G4cout <<"-----------------------" <<G4endl;
00627 if (fragmentT != NULL)
00628 G4cout <<*fragmentT <<G4endl;
00629 else
00630 G4cout <<"(No residual prefragment)" <<G4endl;
00631 }
00632
00633
00634
00635
00636
00637 if (fragmentP !=NULL)
00638 {
00639 G4ReactionProductVector *products = NULL;
00640 if (fragmentP->GetZ() != fragmentP->GetA())
00641 products = theExcitationHandler->BreakItUp(*fragmentP);
00642 else
00643 products = theExcitationHandlerx->BreakItUp(*fragmentP);
00644 delete fragmentP;
00645 fragmentP = NULL;
00646
00647 G4ReactionProductVector::iterator iter;
00648 for (iter = products->begin(); iter != products->end(); ++iter)
00649 {
00650 G4DynamicParticle *secondary =
00651 new G4DynamicParticle((*iter)->GetDefinition(),
00652 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
00653 theParticleChange.AddSecondary (secondary);
00654 G4String particleName = (*iter)->GetDefinition()->GetParticleName();
00655 delete (*iter);
00656 if (verboseLevel >= 2 && particleName.find("[",0) < particleName.size())
00657 {
00658 G4cout <<"------------------------" <<G4endl;
00659 G4cout <<"The projectile fragment:" <<G4endl;
00660 G4cout <<"------------------------" <<G4endl;
00661 G4cout <<" fragmentP = " <<particleName
00662 <<" Energy = " <<secondary->GetKineticEnergy()
00663 <<G4endl;
00664 }
00665 }
00666 delete products;
00667 }
00668
00669
00670
00671
00672
00673
00674 if (fragmentT != NULL)
00675 {
00676 G4ReactionProductVector *products = NULL;
00677 if (fragmentT->GetZ() != fragmentT->GetA())
00678 products = theExcitationHandler->BreakItUp(*fragmentT);
00679 else
00680 products = theExcitationHandlerx->BreakItUp(*fragmentT);
00681 delete fragmentT;
00682 fragmentT = NULL;
00683
00684 G4ReactionProductVector::iterator iter;
00685 for (iter = products->begin(); iter != products->end(); ++iter)
00686 {
00687 G4DynamicParticle *secondary =
00688 new G4DynamicParticle((*iter)->GetDefinition(),
00689 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
00690 theParticleChange.AddSecondary (secondary);
00691 G4String particleName = (*iter)->GetDefinition()->GetParticleName();
00692 delete (*iter);
00693 if (verboseLevel >= 2 && particleName.find("[",0) < particleName.size())
00694 {
00695 G4cout <<"--------------------" <<G4endl;
00696 G4cout <<"The target fragment:" <<G4endl;
00697 G4cout <<"--------------------" <<G4endl;
00698 G4cout <<" fragmentT = " <<particleName
00699 <<" Energy = " <<secondary->GetKineticEnergy()
00700 <<G4endl;
00701 }
00702 }
00703 delete products;
00704 }
00705
00706 if (verboseLevel >= 2)
00707 G4cout <<"########################################"
00708 <<"########################################"
00709 <<G4endl;
00710
00711 delete theAbrasionGeometry;
00712
00713 return &theParticleChange;
00714 }
00716
00717 G4Fragment *G4WilsonAbrasionModel::GetAbradedNucleons (G4int Dabr, G4double A,
00718 G4double Z, G4double r)
00719 {
00720
00721
00722
00723
00724
00725
00726
00727 G4double pK = hbarc * std::pow(9.0 * pi / 4.0 * A, third) / (1.29 * r);
00728 if (A <= 24.0) pK *= -0.229*std::pow(A,third) + 1.62;
00729 G4double pKsq = pK * pK;
00730 G4double p1sq = 2.0/5.0 * pKsq;
00731 G4double p2sq = 6.0/5.0 * pKsq;
00732 G4double p3sq = 500.0 * 500.0;
00733 G4double C1 = 1.0;
00734 G4double C2 = 0.03;
00735 G4double C3 = 0.0002;
00736 G4double gamma = 90.0 * MeV;
00737 G4double maxn = C1 + C2 + C3;
00738
00739
00740
00741
00742
00743 G4double Aabr = 0.0;
00744 G4double Zabr = 0.0;
00745 G4ParticleDefinition *typeNucleon = G4Proton::ProtonDefinition();
00746 G4DynamicParticle *dynamicNucleon = NULL;
00747 G4ParticleMomentum pabr(0.0, 0.0, 0.0);
00748
00749
00750
00751
00752 for (G4int i=0; i<Dabr; i++)
00753 {
00754
00755
00756
00757
00758
00759 G4double p = 0.0;
00760 G4bool found = false;
00761 while (!found)
00762 {
00763 while (p <= 0.0) p = npK * pK * G4UniformRand();
00764 G4double psq = p * p;
00765 found = maxn * G4UniformRand() < C1*std::exp(-psq/p1sq/2.0) +
00766 C2*std::exp(-psq/p2sq/2.0) + C3*std::exp(-psq/p3sq/2.0) + p/gamma/std::sinh(p/gamma);
00767 }
00768
00769
00770
00771
00772
00773
00774 G4double prob = (Z-Zabr)/(A-Aabr);
00775 if (G4UniformRand()<prob)
00776 {
00777 Zabr++;
00778 typeNucleon = G4Proton::ProtonDefinition();
00779 }
00780 else
00781 typeNucleon = G4Neutron::NeutronDefinition();
00782 Aabr++;
00783
00784
00785
00786
00787
00788
00789 G4double costheta = 2.*G4UniformRand()-1.0;
00790 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
00791 G4double phi = 2.0*pi*G4UniformRand()*rad;
00792 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
00793 G4double nucleonMass = typeNucleon->GetPDGMass();
00794 G4double E = std::sqrt(p*p + nucleonMass*nucleonMass)-nucleonMass;
00795 dynamicNucleon = new G4DynamicParticle(typeNucleon,direction,E);
00796 theParticleChange.AddSecondary (dynamicNucleon);
00797 pabr += p*direction;
00798 }
00799
00800
00801
00802
00803
00804
00805
00806 G4Fragment *fragment = NULL;
00807 if (Z-Zabr>=1.0)
00808 {
00809 G4double ionMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
00810 GetIonMass(G4lrint(Z-Zabr),G4lrint(A-Aabr));
00811 G4double E = std::sqrt(pabr.mag2() + ionMass*ionMass);
00812 G4LorentzVector lorentzVector = G4LorentzVector(-pabr, E + 1.0*eV);
00813 fragment =
00814 new G4Fragment((G4int) (A-Aabr), (G4int) (Z-Zabr), lorentzVector);
00815 }
00816
00817 return fragment;
00818 }
00820
00821 G4double G4WilsonAbrasionModel::GetNucleonInducedExcitation
00822 (G4double rP, G4double rT, G4double r)
00823 {
00824
00825
00826
00827
00828 G4double Cl = 0.0;
00829 G4double rPsq = rP * rP;
00830 G4double rTsq = rT * rT;
00831 G4double rsq = r * r;
00832
00833
00834
00835
00836
00837 if (r > rT) Cl = 2.0*std::sqrt(rPsq + 2.0*r*rT - rsq - rTsq);
00838 else Cl = 2.0*rP;
00839
00840
00841
00842
00843
00844
00845
00846 G4double Ct = 0.0;
00847 if (rT > rP && rsq < rTsq - rPsq) Ct = 2.0 * rP;
00848 else if (rP > rT && rsq < rPsq - rTsq) Ct = 2.0 * rT;
00849 else {
00850 G4double bP = (rPsq+rsq-rTsq)/2.0/r;
00851 G4double x = rPsq - bP*bP;
00852 if (x < 0.0) {
00853 G4cerr <<"########################################"
00854 <<"########################################"
00855 <<G4endl;
00856 G4cerr <<"ERROR IN G4WilsonAbrasionModel::GetNucleonInducedExcitation"
00857 <<G4endl;
00858 G4cerr <<"rPsq - bP*bP < 0.0 and cannot be square-rooted" <<G4endl;
00859 G4cerr <<"Set to zero instead" <<G4endl;
00860 G4cerr <<"########################################"
00861 <<"########################################"
00862 <<G4endl;
00863 }
00864 Ct = 2.0*std::sqrt(x);
00865 }
00866
00867 G4double Ex = 13.0 * Cl / fermi;
00868 if (Ct > 1.5*fermi)
00869 Ex += 13.0 * Cl / fermi /3.0 * (Ct/fermi - 1.5);
00870
00871 return Ex;
00872 }
00874
00875 void G4WilsonAbrasionModel::SetUseAblation (G4bool useAblation1)
00876 {
00877 if (useAblation != useAblation1)
00878 {
00879 useAblation = useAblation1;
00880 delete theExcitationHandler;
00881 delete theExcitationHandlerx;
00882 theExcitationHandler = new G4ExcitationHandler;
00883 theExcitationHandlerx = new G4ExcitationHandler;
00884 if (useAblation)
00885 {
00886 theAblation = new G4WilsonAblationModel;
00887 theAblation->SetVerboseLevel(verboseLevel);
00888 theExcitationHandler->SetEvaporation(theAblation);
00889 theExcitationHandlerx->SetEvaporation(theAblation);
00890 }
00891 else
00892 {
00893 theAblation = NULL;
00894 G4Evaporation * theEvaporation = new G4Evaporation;
00895 G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp;
00896 G4StatMF * theMF = new G4StatMF;
00897 theExcitationHandler->SetEvaporation(theEvaporation);
00898 theExcitationHandler->SetFermiModel(theFermiBreakUp);
00899 theExcitationHandler->SetMultiFragmentation(theMF);
00900 theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6);
00901 theExcitationHandler->SetMinEForMultiFrag(5.0*MeV);
00902
00903 theEvaporation = new G4Evaporation;
00904 theFermiBreakUp = new G4FermiBreakUp;
00905 theExcitationHandlerx->SetEvaporation(theEvaporation);
00906 theExcitationHandlerx->SetFermiModel(theFermiBreakUp);
00907 theExcitationHandlerx->SetMaxAandZForFermiBreakUp(12, 6);
00908 }
00909 }
00910 return;
00911 }
00913
00914 void G4WilsonAbrasionModel::PrintWelcomeMessage ()
00915 {
00916 G4cout <<G4endl;
00917 G4cout <<" *****************************************************************"
00918 <<G4endl;
00919 G4cout <<" Nuclear abrasion model for nuclear-nuclear interactions activated"
00920 <<G4endl;
00921 G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)"
00922 <<G4endl;
00923 G4cout <<" *****************************************************************"
00924 <<G4endl;
00925 G4cout << G4endl;
00926
00927 return;
00928 }
00930