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 #include "G4KaonMinusAbsorptionAtRest.hh"
00040
00041 #include "G4StopDeexcitation.hh"
00042 #include "G4StopTheoDeexcitation.hh"
00043 #include "G4StopDeexcitationAlgorithm.hh"
00044 #include "G4ReactionKinematics.hh"
00045 #include "G4HadronicProcessStore.hh"
00046 #include "G4HadronicDeprecate.hh"
00047
00048 G4KaonMinusAbsorptionAtRest::G4KaonMinusAbsorptionAtRest(const G4String& processName,
00049 G4ProcessType aType ) :
00050 G4VRestProcess (processName, aType)
00051 {
00052 G4HadronicDeprecate("G4KaonMinusAbsorptionAtRest");
00053 if (verboseLevel>0) {
00054 G4cout << GetProcessName() << " is created "<< G4endl;
00055 }
00056 SetProcessSubType(fHadronAtRest);
00057
00058
00059
00060
00061 pionAbsorptionRate = 0.07;
00062
00063
00064
00065
00066
00067
00068
00069 rateLambdaZeroPiZero = 0.052;
00070 rateSigmaMinusPiPlus = 0.199;
00071 rateSigmaPlusPiMinus = 0.446;
00072 rateSigmaZeroPiZero = 0.303;
00073 rateLambdaZeroPiMinus = 0.568;
00074 rateSigmaZeroPiMinus = 0.216;
00075 rateSigmaMinusPiZero = 0.216;
00076
00077
00078
00079
00080
00081
00082 sigmaPlusLambdaConversionRate = 0.55;
00083 sigmaMinusLambdaConversionRate = 0.55;
00084 sigmaZeroLambdaConversionRate = 0.55;
00085
00086 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
00087 }
00088
00089
00090 G4KaonMinusAbsorptionAtRest::~G4KaonMinusAbsorptionAtRest()
00091 {
00092 G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
00093 }
00094
00095 void G4KaonMinusAbsorptionAtRest::PreparePhysicsTable(const G4ParticleDefinition& p)
00096 {
00097 G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
00098 }
00099
00100 void G4KaonMinusAbsorptionAtRest::BuildPhysicsTable(const G4ParticleDefinition& p)
00101 {
00102 G4HadronicProcessStore::Instance()->PrintInfo(&p);
00103 }
00104
00105 G4VParticleChange* G4KaonMinusAbsorptionAtRest::AtRestDoIt
00106 (const G4Track& track, const G4Step& )
00107 {
00108 stoppedHadron = track.GetDynamicParticle();
00109
00110
00111
00112 if (!IsApplicable(*(stoppedHadron->GetDefinition())))
00113 {
00114 G4cerr <<"G4KaonMinusAbsorptionAtRest:ERROR, particle must be a Kaon!" <<G4endl;
00115 return 0;
00116 }
00117
00118 G4Material* material;
00119 material = track.GetMaterial();
00120 nucleus = 0;
00121 do
00122 {
00123
00124 nucleus = new G4Nucleus(material);
00125 if (nucleus->GetA_asInt() < 1.5)
00126 {
00127 delete nucleus;
00128 nucleus = 0;
00129 }
00130 } while(nucleus == 0);
00131
00132 G4double Z = nucleus->GetZ_asInt();
00133 G4double A = nucleus->GetA_asInt();
00134
00135
00136 G4DynamicParticleVector* absorptionProducts = KaonNucleonReaction();
00137
00138
00139 if ( ! absorptionProducts ) {
00140 G4Exception("G4KaonMinusAbsorptionAtRest::AtRestDoIt()", "HAD_STOP_0001",
00141 FatalException, "NULL absorptionProducts");
00142 return 0;
00143 }
00144
00145
00146
00147 G4DynamicParticle* thePion;
00148 unsigned int i;
00149 for(i = 0; i < absorptionProducts->size(); i++)
00150 {
00151 thePion = (*absorptionProducts)[i];
00152 if (thePion->GetDefinition() == G4PionMinus::PionMinus()
00153 || thePion->GetDefinition() == G4PionPlus::PionPlus()
00154 || thePion->GetDefinition() == G4PionZero::PionZero())
00155 {
00156 if (AbsorbPionByNucleus(thePion))
00157 {
00158 absorptionProducts->erase(absorptionProducts->begin()+i);
00159 i--;
00160 delete thePion;
00161 if (verboseLevel > 1)
00162 G4cout << "G4KaonMinusAbsorption::AtRestDoIt: Pion absorbed in Nucleus"
00163 << G4endl;
00164 }
00165 }
00166 }
00167
00168 G4DynamicParticle* theSigma;
00169 G4DynamicParticle* theLambda;
00170 for (i = 0; i < absorptionProducts->size(); i++)
00171 {
00172 theSigma = (*absorptionProducts)[i];
00173 if (theSigma->GetDefinition() == G4SigmaMinus::SigmaMinus()
00174 || theSigma->GetDefinition() == G4SigmaPlus::SigmaPlus()
00175 || theSigma->GetDefinition() == G4SigmaZero::SigmaZero())
00176 {
00177 theLambda = SigmaLambdaConversion(theSigma);
00178 if (theLambda != 0){
00179 absorptionProducts->erase(absorptionProducts->begin()+i);
00180 i--;
00181 delete theSigma;
00182 absorptionProducts->push_back(theLambda);
00183
00184 if (verboseLevel > 1)
00185 G4cout << "G4KaonMinusAbsorption::AtRestDoIt: SigmaLambdaConversion Done"
00186 << G4endl;
00187 }
00188 }
00189 }
00190
00191
00192
00193 G4double productEnergy = 0.;
00194 G4ThreeVector pProducts(0.,0.,0.);
00195
00196 unsigned int nAbsorptionProducts = 0;
00197 if (absorptionProducts != 0) nAbsorptionProducts = absorptionProducts->size();
00198
00199 for ( i = 0; i<nAbsorptionProducts; i++)
00200 {
00201 pProducts += (*absorptionProducts)[i]->GetMomentum();
00202 productEnergy += (*absorptionProducts)[i]->GetKineticEnergy();
00203 }
00204
00205 G4double newZ = nucleus->GetZ_asInt();
00206 G4double newA = nucleus->GetA_asInt();
00207
00208 G4double bDiff = G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(A),static_cast<G4int>(Z)) -
00209 G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(newA), static_cast<G4int>(newZ));
00210
00211 G4StopDeexcitationAlgorithm* nucleusAlgorithm = new G4StopTheoDeexcitation();
00212 G4StopDeexcitation stopDeexcitation(nucleusAlgorithm);
00213
00214 nucleus->AddExcitationEnergy(bDiff);
00215
00216
00217 G4double energyDeposit = nucleus->GetEnergyDeposit();
00218 if (verboseLevel>0)
00219 {
00220 G4cout << " -- KaonAtRest -- excitation = "
00221 << energyDeposit
00222 << ", pNucleus = "
00223 << pProducts
00224 << ", A: "
00225 << A
00226 << ", "
00227 << newA
00228 << ", Z: "
00229 << Z
00230 << ", "
00231 << newZ
00232 << G4endl;
00233 }
00234
00235 if (energyDeposit < 0.)
00236 G4Exception("G4KaonMinusAbsorptionAtRest::AtRestDoIt()", "HAD_STOP_0002",
00237 FatalException, "Excitation energy < 0");
00238 delete nucleus;
00239
00240 G4ReactionProductVector* fragmentationProducts = stopDeexcitation.DoBreakUp(newA,newZ,energyDeposit,pProducts);
00241
00242 unsigned int nFragmentationProducts = 0;
00243 if (fragmentationProducts != 0) nFragmentationProducts = fragmentationProducts->size();
00244
00245
00246 aParticleChange.Initialize(track);
00247 aParticleChange.SetNumberOfSecondaries(G4int(nAbsorptionProducts+nFragmentationProducts) );
00248
00249
00250 for (i = 0; i < nAbsorptionProducts; i++)
00251 {aParticleChange.AddSecondary((*absorptionProducts)[i]); }
00252 if (absorptionProducts != 0) delete absorptionProducts;
00253
00254
00255
00256 for(i=0; i<nFragmentationProducts; i++)
00257 {
00258 G4DynamicParticle * aNew =
00259 new G4DynamicParticle((*fragmentationProducts)[i]->GetDefinition(),
00260 (*fragmentationProducts)[i]->GetTotalEnergy(),
00261 (*fragmentationProducts)[i]->GetMomentum());
00262 G4double newTime = aParticleChange.GetGlobalTime((*fragmentationProducts)[i]->GetFormationTime());
00263 aParticleChange.AddSecondary(aNew, newTime);
00264 delete (*fragmentationProducts)[i];
00265 }
00266 if (fragmentationProducts != 0) delete fragmentationProducts;
00267
00268
00269 aParticleChange.ProposeTrackStatus(fStopAndKill);
00270 return &aParticleChange;
00271 }
00272
00273
00274 G4DynamicParticle G4KaonMinusAbsorptionAtRest::GetAbsorbingNucleon()
00275 {
00276 G4DynamicParticle aNucleon;
00277
00278
00279 aNucleon.SetDefinition(SelectAbsorbingNucleon());
00280
00281
00282 G4ThreeVector pFermi = nucleus->GetFermiMomentum();
00283 aNucleon.SetMomentum(pFermi);
00284
00285 return aNucleon;
00286 }
00287
00288 G4ParticleDefinition* G4KaonMinusAbsorptionAtRest::SelectAbsorbingNucleon()
00289 {
00290
00291
00292
00293
00294
00295 G4ParticleDefinition* absorbingParticleDef;
00296
00297 G4double ranflat = G4UniformRand();
00298
00299 G4double myZ = nucleus->GetZ_asInt();
00300 G4double myN = nucleus->GetA_asInt();
00301
00302
00303 G4double carbonRatioNP = 0.18;
00304
00305 G4double neutronProtonRatio = NeutronHaloFactor(myZ,myN)*carbonRatioNP*(myN-myZ)/myZ;
00306 G4double protonProbability = 1./(1.+neutronProtonRatio);
00307
00308 if ( ranflat < protonProbability )
00309 {
00310 absorbingParticleDef = G4Proton::Proton();
00311 myZ-= 1.;
00312 }
00313 else
00314 { absorbingParticleDef = G4Neutron::Neutron(); }
00315
00316 myN -= 1.;
00317 nucleus->SetParameters(myN,myZ);
00318 return absorbingParticleDef;
00319 }
00320
00321
00322 G4double G4KaonMinusAbsorptionAtRest::NeutronHaloFactor(G4double Z, G4double N)
00323 {
00324
00325
00326
00327
00328
00329
00330 if (Z == 1.) return 1.389;
00331 else if (Z == 2.) return 1.78;
00332 else if (Z == 10.) return 0.66;
00333 else
00334 return 0.6742+(N-Z)*0.06524;
00335 }
00336
00337
00338 G4DynamicParticleVector* G4KaonMinusAbsorptionAtRest::KaonNucleonReaction()
00339 {
00340 G4DynamicParticleVector* products = new G4DynamicParticleVector();
00341
00342 G4double ranflat = G4UniformRand();
00343 G4double prob = 0;
00344
00345 G4ParticleDefinition* producedBaryonDef;
00346 G4ParticleDefinition* producedMesonDef;
00347
00348 G4double iniZ = nucleus->GetZ_asInt();
00349 G4double iniA = nucleus->GetA_asInt();
00350
00351 G4DynamicParticle aNucleon = GetAbsorbingNucleon();
00352
00353
00354
00355 if (aNucleon.GetDefinition() == G4Proton::Proton())
00356 {
00357
00358 if ( (prob += rateLambdaZeroPiZero) > ranflat)
00359 {
00360 producedBaryonDef = G4Lambda::Lambda();
00361 producedMesonDef = G4PionZero::PionZero();
00362 }
00363 else if ((prob += rateSigmaPlusPiMinus) > ranflat)
00364 {
00365 producedBaryonDef = G4SigmaPlus::SigmaPlus();
00366 producedMesonDef = G4PionMinus::PionMinus();
00367 }
00368 else if ((prob += rateSigmaMinusPiPlus) > ranflat)
00369 {
00370 producedBaryonDef = G4SigmaMinus::SigmaMinus();
00371 producedMesonDef = G4PionPlus::PionPlus();
00372 }
00373 else
00374 {
00375 producedBaryonDef = G4SigmaZero::SigmaZero();
00376 producedMesonDef = G4PionZero::PionZero();
00377 }
00378 }
00379 else if (aNucleon.GetDefinition() == G4Neutron::Neutron())
00380 {
00381
00382 if ((prob += rateLambdaZeroPiMinus) > ranflat)
00383 {
00384 producedBaryonDef = G4Lambda::Lambda();
00385 producedMesonDef = G4PionMinus::PionMinus();
00386 }
00387 else if ((prob += rateSigmaZeroPiMinus) > ranflat)
00388 {
00389 producedBaryonDef = G4SigmaZero::SigmaZero();
00390 producedMesonDef = G4PionMinus::PionMinus();
00391 }
00392 else
00393 {
00394 producedBaryonDef = G4SigmaMinus::SigmaMinus();
00395 producedMesonDef = G4PionZero::PionZero();
00396 }
00397 }
00398 else
00399 {
00400 if (verboseLevel>0)
00401 {
00402 G4cout
00403 << "G4KaonMinusAbsorption::KaonNucleonReaction: "
00404 << aNucleon.GetDefinition()->GetParticleName()
00405 << " is not a good nucleon - check G4Nucleus::ReturnTargetParticle()!"
00406 << G4endl;
00407 }
00408
00409
00410 if ( products ) delete products;
00411
00412 return 0;
00413 }
00414
00415 G4double newZ = nucleus->GetZ_asInt();
00416 G4double newA = nucleus->GetA_asInt();
00417
00418
00419
00420
00421
00422
00423 G4double nucleonBindingEnergy =
00424 -G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(iniA), static_cast<G4int>(iniZ) )
00425 +G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(newA), static_cast<G4int>(newZ) );
00426
00427 G4DynamicParticle modifiedHadron = (*stoppedHadron);
00428 modifiedHadron.SetMass(stoppedHadron->GetMass() + nucleonBindingEnergy);
00429
00430
00431 G4ThreeVector dummy(0.,0.,0.);
00432 G4DynamicParticle* producedBaryon = new G4DynamicParticle(producedBaryonDef,dummy);
00433 G4DynamicParticle* producedMeson = new G4DynamicParticle(producedMesonDef,dummy);
00434
00435
00436 G4ReactionKinematics theReactionKinematics;
00437 theReactionKinematics.TwoBodyScattering( &modifiedHadron, &aNucleon,
00438 producedBaryon, producedMeson);
00439
00440 products->push_back(producedBaryon);
00441 products->push_back(producedMeson);
00442
00443 if (verboseLevel > 1)
00444 {
00445 G4cout
00446 << "G4KaonMinusAbsorption::KaonNucleonReaction: Number of primaries = "
00447 << products->size()
00448 << ": " <<producedMesonDef->GetParticleName()
00449 << ", " <<producedBaryonDef->GetParticleName() << G4endl;
00450 }
00451
00452 return products;
00453 }
00454
00455
00456 G4bool G4KaonMinusAbsorptionAtRest::AbsorbPionByNucleus(G4DynamicParticle* aPion)
00457 {
00458
00459
00460 G4double ranflat = G4UniformRand();
00461
00462 if (ranflat < pionAbsorptionRate){
00463
00464 nucleus->AddExcitationEnergy(aPion->GetTotalEnergy());
00465 nucleus->AddMomentum(aPion->GetMomentum());
00466 }
00467
00468 return (ranflat < pionAbsorptionRate);
00469 }
00470
00471 G4DynamicParticle* G4KaonMinusAbsorptionAtRest::SigmaLambdaConversion(G4DynamicParticle* aSigma)
00472 {
00473 G4double ranflat = G4UniformRand();
00474 G4double sigmaLambdaConversionRate;
00475
00476 G4double A = nucleus->GetA_asInt();
00477 G4double Z = nucleus->GetZ_asInt();
00478
00479 G4double newZ = Z;
00480
00481
00482 G4ParticleDefinition* inNucleonDef=NULL;
00483 G4ParticleDefinition* outNucleonDef=NULL;
00484
00485
00486 switch((int) aSigma->GetDefinition()->GetPDGCharge()) {
00487
00488 case 1:
00489 sigmaLambdaConversionRate = sigmaPlusLambdaConversionRate;
00490 inNucleonDef = G4Neutron::Neutron();
00491 outNucleonDef = G4Proton::Proton();
00492 newZ = Z+1;
00493
00494 break;
00495
00496 case -1:
00497 sigmaLambdaConversionRate = sigmaMinusLambdaConversionRate;
00498 inNucleonDef = G4Proton::Proton();
00499 outNucleonDef = G4Neutron::Neutron();
00500 newZ = Z-1;
00501
00502 break;
00503
00504 case 0:
00505 sigmaLambdaConversionRate = sigmaZeroLambdaConversionRate;
00506
00507
00508
00509 inNucleonDef = G4Neutron::Neutron();
00510 outNucleonDef = G4Neutron::Neutron();
00511 break;
00512
00513 default:
00514 sigmaLambdaConversionRate = 0.;
00515
00516 inNucleonDef = G4Proton::Proton();
00517 outNucleonDef = G4Proton::Proton();
00518 }
00519
00520 if (ranflat >= sigmaLambdaConversionRate) return 0;
00521
00522 G4ThreeVector dummy(0.,0.,0.);
00523
00524
00525 G4ThreeVector momentum = nucleus->GetFermiMomentum();
00526
00527 G4ParticleDefinition* lambdaDef = G4Lambda::Lambda();
00528
00529 G4DynamicParticle inNucleon(inNucleonDef,momentum);
00530 G4DynamicParticle outNucleon(outNucleonDef,dummy);
00531 G4DynamicParticle* outLambda = new G4DynamicParticle(lambdaDef,dummy);
00532
00533 G4ReactionKinematics theReactionKinematics;
00534
00535
00536 theReactionKinematics.TwoBodyScattering(aSigma, &inNucleon,
00537 &outNucleon, outLambda);
00538
00539
00540
00541
00542
00543
00544
00545
00546 nucleus->AddExcitationEnergy(outNucleon.GetKineticEnergy());
00547 nucleus->AddMomentum(outNucleon.GetMomentum());
00548 nucleus->SetParameters(A,newZ);
00549
00550
00551 return outLambda;
00552 }
00553
00554
00555