source: trunk/source/processes/hadronic/stopping/src/G4KaonMinusAbsorptionAtRest.cc@ 1330

Last change on this file since 1330 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 17.6 KB
Line 
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// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// Author: Christian V"olcker (Christian.Volcker@cern.ch),
27//
28// Creation date: November 1997
29//
30// Testfile: ../G4KaonMinusAbsorptionAtRestTest.cc
31//
32// Modifications:
33// Maria Grazia Pia September 1998
34// Various bug fixes, eliminated several memory leaks
35//
36// -------------------------------------------------------------------
37
38
39#include "G4KaonMinusAbsorptionAtRest.hh"
40
41#include "G4StopDeexcitation.hh"
42#include "G4StopTheoDeexcitation.hh"
43#include "G4StopDeexcitationAlgorithm.hh"
44#include "G4ReactionKinematics.hh"
45#include "G4HadronicProcessStore.hh"
46
47G4KaonMinusAbsorptionAtRest::G4KaonMinusAbsorptionAtRest(const G4String& processName,
48 G4ProcessType aType ) :
49 G4VRestProcess (processName, aType)
50{
51 if (verboseLevel>0) {
52 G4cout << GetProcessName() << " is created "<< G4endl;
53 }
54 SetProcessSubType(fHadronAtRest);
55
56 // see Cohn et al, PLB27(1968) 527;
57 // Davis et al, PLB1(1967) 434;
58
59 pionAbsorptionRate = 0.07;
60
61 // see VanderVelde-Wilquet et al, Nuov.Cim.39A(1978)538;
62 // see VanderVelde-Wilquet et al, Nuov.Cim.38A(1977)178;
63 // see VanderVelde-Wilquet et al, Nucl.Phys.A241(1975)511;
64 // primary production rates ( for absorption on Carbon)
65 // .. other elements are extrapolated by the halo factor.
66
67 rateLambdaZeroPiZero = 0.052;
68 rateSigmaMinusPiPlus = 0.199;
69 rateSigmaPlusPiMinus = 0.446;
70 rateSigmaZeroPiZero = 0.303;
71 rateLambdaZeroPiMinus = 0.568;
72 rateSigmaZeroPiMinus = 0.216;
73 rateSigmaMinusPiZero = 0.216;
74
75 // for sigma- p -> lambda n
76 // sigma+ n -> lambda p
77 // sigma- n -> lambda
78 // all values compatible with 0.55 same literature as above.
79
80 sigmaPlusLambdaConversionRate = 0.55;
81 sigmaMinusLambdaConversionRate = 0.55;
82 sigmaZeroLambdaConversionRate = 0.55;
83
84 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
85}
86
87
88G4KaonMinusAbsorptionAtRest::~G4KaonMinusAbsorptionAtRest()
89{
90 G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
91}
92
93void G4KaonMinusAbsorptionAtRest::PreparePhysicsTable(const G4ParticleDefinition& p)
94{
95 G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
96}
97
98void G4KaonMinusAbsorptionAtRest::BuildPhysicsTable(const G4ParticleDefinition& p)
99{
100 G4HadronicProcessStore::Instance()->PrintInfo(&p);
101}
102
103G4VParticleChange* G4KaonMinusAbsorptionAtRest::AtRestDoIt
104(const G4Track& track, const G4Step& )
105{
106 stoppedHadron = track.GetDynamicParticle();
107
108 // Check applicability
109
110 if (!IsApplicable(*(stoppedHadron->GetDefinition())))
111 {
112 G4cerr <<"G4KaonMinusAbsorptionAtRest:ERROR, particle must be a Kaon!" <<G4endl;
113 return 0;
114 }
115
116 G4Material* material;
117 material = track.GetMaterial();
118 nucleus = 0;
119 do
120 {
121 // Select the nucleus, get nucleon
122 nucleus = new G4Nucleus(material);
123 if (nucleus->GetN() < 1.5)
124 {
125 delete nucleus;
126 nucleus = 0;
127 }
128 } while(nucleus == 0);
129
130 G4double Z = nucleus->GetZ();
131 G4double A = nucleus->GetN();
132
133 // Do the interaction with the nucleon
134 G4DynamicParticleVector* absorptionProducts = KaonNucleonReaction();
135
136 // Secondary interactions
137
138 G4DynamicParticle* thePion;
139 unsigned int i;
140 for(i = 0; i < absorptionProducts->size(); i++)
141 {
142 thePion = (*absorptionProducts)[i];
143 if (thePion->GetDefinition() == G4PionMinus::PionMinus()
144 || thePion->GetDefinition() == G4PionPlus::PionPlus()
145 || thePion->GetDefinition() == G4PionZero::PionZero())
146 {
147 if (AbsorbPionByNucleus(thePion))
148 {
149 absorptionProducts->erase(absorptionProducts->begin()+i);
150 i--;
151 delete thePion;
152 if (verboseLevel > 1)
153 G4cout << "G4KaonMinusAbsorption::AtRestDoIt: Pion absorbed in Nucleus"
154 << G4endl;
155 }
156 }
157 }
158
159 G4DynamicParticle* theSigma;
160 G4DynamicParticle* theLambda;
161 for (i = 0; i < absorptionProducts->size(); i++)
162 {
163 theSigma = (*absorptionProducts)[i];
164 if (theSigma->GetDefinition() == G4SigmaMinus::SigmaMinus()
165 || theSigma->GetDefinition() == G4SigmaPlus::SigmaPlus()
166 || theSigma->GetDefinition() == G4SigmaZero::SigmaZero())
167 {
168 theLambda = SigmaLambdaConversion(theSigma);
169 if (theLambda != 0){
170 absorptionProducts->erase(absorptionProducts->begin()+i);
171 i--;
172 delete theSigma;
173 absorptionProducts->push_back(theLambda);
174
175 if (verboseLevel > 1)
176 G4cout << "G4KaonMinusAbsorption::AtRestDoIt: SigmaLambdaConversion Done"
177 << G4endl;
178 }
179 }
180 }
181
182 // Nucleus deexcitation
183
184 G4double productEnergy = 0.;
185 G4ThreeVector pProducts(0.,0.,0.);
186
187 unsigned int nAbsorptionProducts = 0;
188 if (absorptionProducts != 0) nAbsorptionProducts = absorptionProducts->size();
189
190 for ( i = 0; i<nAbsorptionProducts; i++)
191 {
192 pProducts += (*absorptionProducts)[i]->GetMomentum();
193 productEnergy += (*absorptionProducts)[i]->GetKineticEnergy();
194 }
195
196 G4double newZ = nucleus->GetZ();
197 G4double newA = nucleus->GetN();
198
199 G4double bDiff = G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(A),static_cast<G4int>(Z)) -
200 G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(newA), static_cast<G4int>(newZ));
201
202 G4StopDeexcitationAlgorithm* nucleusAlgorithm = new G4StopTheoDeexcitation();
203 G4StopDeexcitation stopDeexcitation(nucleusAlgorithm);
204
205 nucleus->AddExcitationEnergy(bDiff);
206
207 // returns excitation energy for the moment ..
208 G4double energyDeposit = nucleus->GetEnergyDeposit();
209 if (verboseLevel>0)
210 {
211 G4cout << " -- KaonAtRest -- excitation = "
212 << energyDeposit
213 << ", pNucleus = "
214 << pProducts
215 << ", A: "
216 << A
217 << ", "
218 << newA
219 << ", Z: "
220 << Z
221 << ", "
222 << newZ
223 << G4endl;
224 }
225
226 if (energyDeposit < 0.)
227 {
228 G4Exception("G4KaonMinusAbsorptionAtRest", "007", FatalException,
229 "AtRestDoIt -- excitation energy < 0");
230 }
231 delete nucleus;
232
233 G4ReactionProductVector* fragmentationProducts = stopDeexcitation.DoBreakUp(newA,newZ,energyDeposit,pProducts);
234
235 unsigned int nFragmentationProducts = 0;
236 if (fragmentationProducts != 0) nFragmentationProducts = fragmentationProducts->size();
237
238 //Initialize ParticleChange
239 aParticleChange.Initialize(track);
240 aParticleChange.SetNumberOfSecondaries(G4int(nAbsorptionProducts+nFragmentationProducts) );
241
242 // update List of alive particles. put energy deposit at the right place ...
243 for (i = 0; i < nAbsorptionProducts; i++)
244 {aParticleChange.AddSecondary((*absorptionProducts)[i]); }
245 if (absorptionProducts != 0) delete absorptionProducts;
246
247// for (i = 0; i < nFragmentationProducts; i++)
248// { aParticleChange.AddSecondary(fragmentationProducts->at(i)); }
249 for(i=0; i<nFragmentationProducts; i++)
250 {
251 G4DynamicParticle * aNew =
252 new G4DynamicParticle((*fragmentationProducts)[i]->GetDefinition(),
253 (*fragmentationProducts)[i]->GetTotalEnergy(),
254 (*fragmentationProducts)[i]->GetMomentum());
255 G4double newTime = aParticleChange.GetGlobalTime((*fragmentationProducts)[i]->GetFormationTime());
256 aParticleChange.AddSecondary(aNew, newTime);
257 delete (*fragmentationProducts)[i];
258 }
259 if (fragmentationProducts != 0) delete fragmentationProducts;
260
261 // finally ...
262 aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the incident Kaon
263 return &aParticleChange;
264}
265
266
267G4DynamicParticle G4KaonMinusAbsorptionAtRest::GetAbsorbingNucleon()
268{
269 G4DynamicParticle aNucleon;
270
271 // Get nucleon definition, based on Z,N of current Nucleus
272 aNucleon.SetDefinition(SelectAbsorbingNucleon());
273
274 // Fermi momentum distribution in three dimensions
275 G4ThreeVector pFermi = nucleus->GetFermiMomentum();
276 aNucleon.SetMomentum(pFermi);
277
278 return aNucleon;
279}
280
281G4ParticleDefinition* G4KaonMinusAbsorptionAtRest::SelectAbsorbingNucleon()
282{
283 // (Ch. Voelcker) extended from ReturnTargetParticle():
284 // Choose a proton or a neutron as the absorbing particle,
285 // taking weight into account!
286 // Update nucleon's atomic numbers.
287
288 G4ParticleDefinition* absorbingParticleDef;
289
290 G4double ranflat = G4UniformRand();
291
292 G4double myZ = nucleus->GetZ(); // number of protons
293 G4double myN = nucleus->GetN(); // number of nucleons (not neutrons!!)
294
295 // See VanderVelde-Wilquet et al, Nuov.Cim.39A(1978)538;
296 G4double carbonRatioNP = 0.18; // (Rn/Rp)c, see page 544
297
298 G4double neutronProtonRatio = NeutronHaloFactor(myZ,myN)*carbonRatioNP*(myN-myZ)/myZ;
299 G4double protonProbability = 1./(1.+neutronProtonRatio);
300
301 if ( ranflat < protonProbability )
302 {
303 absorbingParticleDef = G4Proton::Proton();
304 myZ-= 1.;
305 }
306 else
307 { absorbingParticleDef = G4Neutron::Neutron(); }
308
309 myN -= 1.;
310 nucleus->SetParameters(myN,myZ);
311 return absorbingParticleDef;
312}
313
314
315G4double G4KaonMinusAbsorptionAtRest::NeutronHaloFactor(G4double Z, G4double N)
316{
317 // this function should take care of the probability for absorption
318 // on neutrons, depending on number of protons Z and number of neutrons N-Z
319 // parametrisation from fit to
320 // VanderVelde-Wilquet et al, Nuov.Cim.39A(1978)538;
321 //
322
323 if (Z == 1.) return 1.389; // deuterium
324 else if (Z == 2.) return 1.78; // helium
325 else if (Z == 10.) return 0.66; // neon
326 else
327 return 0.6742+(N-Z)*0.06524;
328}
329
330
331G4DynamicParticleVector* G4KaonMinusAbsorptionAtRest::KaonNucleonReaction()
332{
333 G4DynamicParticleVector* products = new G4DynamicParticleVector();
334
335 G4double ranflat = G4UniformRand();
336 G4double prob = 0;
337
338 G4ParticleDefinition* producedBaryonDef;
339 G4ParticleDefinition* producedMesonDef;
340
341 G4double iniZ = nucleus->GetZ();
342 G4double iniA = nucleus->GetN();
343
344 G4DynamicParticle aNucleon = GetAbsorbingNucleon();
345
346 G4double nucleonMass;
347
348 if (aNucleon.GetDefinition() == G4Proton::Proton())
349 {
350 nucleonMass = proton_mass_c2+electron_mass_c2;
351 if ( (prob += rateLambdaZeroPiZero) > ranflat)
352 { // lambda pi0
353 producedBaryonDef = G4Lambda::Lambda();
354 producedMesonDef = G4PionZero::PionZero();
355 }
356 else if ((prob += rateSigmaPlusPiMinus) > ranflat)
357 { // sigma+ pi-
358 producedBaryonDef = G4SigmaPlus::SigmaPlus();
359 producedMesonDef = G4PionMinus::PionMinus();
360 }
361 else if ((prob += rateSigmaMinusPiPlus) > ranflat)
362 { // sigma- pi+
363 producedBaryonDef = G4SigmaMinus::SigmaMinus();
364 producedMesonDef = G4PionPlus::PionPlus();
365 }
366 else
367 { // sigma0 pi0
368 producedBaryonDef = G4SigmaZero::SigmaZero();
369 producedMesonDef = G4PionZero::PionZero();
370 }
371 }
372 else if (aNucleon.GetDefinition() == G4Neutron::Neutron())
373 {
374 nucleonMass = neutron_mass_c2;
375 if ((prob += rateLambdaZeroPiMinus) > ranflat)
376 { // lambda pi-
377 producedBaryonDef = G4Lambda::Lambda();
378 producedMesonDef = G4PionMinus::PionMinus();
379 }
380 else if ((prob += rateSigmaZeroPiMinus) > ranflat)
381 { // sigma0 pi-
382 producedBaryonDef = G4SigmaZero::SigmaZero();
383 producedMesonDef = G4PionMinus::PionMinus();
384 }
385 else
386 { // sigma- pi0
387 producedBaryonDef = G4SigmaMinus::SigmaMinus();
388 producedMesonDef = G4PionZero::PionZero();
389 }
390 }
391 else
392 {
393 if (verboseLevel>0)
394 {
395 G4cout
396 << "G4KaonMinusAbsorption::KaonNucleonReaction: "
397 << aNucleon.GetDefinition()->GetParticleName()
398 << " is not a good nucleon - check G4Nucleus::ReturnTargetParticle()!"
399 << G4endl;
400 }
401 return 0;
402 }
403
404 G4double newZ = nucleus->GetZ();
405 G4double newA = nucleus->GetN();
406
407 // Modify the Kaon mass to take nuclear binding energy into account
408 // .. using mas formula ..
409 // .. using mass table ..
410 // equivalent to -'initialBindingEnergy+nucleus.GetBindingEnergy' !
411
412 G4double nucleonBindingEnergy =
413 -G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(iniA), static_cast<G4int>(iniZ) )
414 +G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(newA), static_cast<G4int>(newZ) );
415
416 G4DynamicParticle modifiedHadron = (*stoppedHadron);
417 modifiedHadron.SetMass(stoppedHadron->GetMass() + nucleonBindingEnergy);
418
419 // Setup outgoing dynamic particles
420 G4ThreeVector dummy(0.,0.,0.);
421 G4DynamicParticle* producedBaryon = new G4DynamicParticle(producedBaryonDef,dummy);
422 G4DynamicParticle* producedMeson = new G4DynamicParticle(producedMesonDef,dummy);
423
424 // Produce the secondary particles in a twobody process:
425 G4ReactionKinematics theReactionKinematics;
426 theReactionKinematics.TwoBodyScattering( &modifiedHadron, &aNucleon,
427 producedBaryon, producedMeson);
428
429 products->push_back(producedBaryon);
430 products->push_back(producedMeson);
431
432 if (verboseLevel > 1)
433 {
434 G4cout
435 << "G4KaonMinusAbsorption::KaonNucleonReaction: Number of primaries = "
436 << products->size()
437 << ": " <<producedMesonDef->GetParticleName()
438 << ", " <<producedBaryonDef->GetParticleName() << G4endl;
439 }
440
441 return products;
442}
443
444
445G4bool G4KaonMinusAbsorptionAtRest::AbsorbPionByNucleus(G4DynamicParticle* aPion)
446{
447 // Needs some more investigation!
448
449 G4double ranflat = G4UniformRand();
450
451 if (ranflat < pionAbsorptionRate){
452 // Add pion energy to ExcitationEnergy and NucleusMomentum
453 nucleus->AddExcitationEnergy(aPion->GetTotalEnergy());
454 nucleus->AddMomentum(aPion->GetMomentum());
455 }
456
457 return (ranflat < pionAbsorptionRate);
458}
459
460G4DynamicParticle* G4KaonMinusAbsorptionAtRest::SigmaLambdaConversion(G4DynamicParticle* aSigma)
461{
462 G4double ranflat = G4UniformRand();
463 G4double sigmaLambdaConversionRate;
464
465 G4double A = nucleus->GetN();
466 G4double Z = nucleus->GetZ();
467
468 G4double newZ = Z;
469 G4double nucleonMassDifference = 0;
470
471 G4ParticleDefinition* inNucleonDef=NULL;
472 G4ParticleDefinition* outNucleonDef=NULL;
473
474 // Decide which sigma
475 switch((int) aSigma->GetDefinition()->GetPDGCharge()) {
476
477 case 1:
478 sigmaLambdaConversionRate = sigmaPlusLambdaConversionRate;
479 inNucleonDef = G4Neutron::Neutron();
480 outNucleonDef = G4Proton::Proton();
481 newZ = Z+1;
482 nucleonMassDifference = neutron_mass_c2 - proton_mass_c2-electron_mass_c2;
483 break;
484
485 case -1:
486 sigmaLambdaConversionRate = sigmaMinusLambdaConversionRate;
487 inNucleonDef = G4Proton::Proton();
488 outNucleonDef = G4Neutron::Neutron();
489 newZ = Z-1;
490 nucleonMassDifference = proton_mass_c2+electron_mass_c2 - neutron_mass_c2;
491 break;
492
493 case 0:
494 sigmaLambdaConversionRate = sigmaZeroLambdaConversionRate;
495 // The 'outgoing' nucleon is just virtual, to keep the energy-momentum
496 // balance and will not appear in the ParticleChange. Therefore no need
497 // choose between neutron and proton here!
498 inNucleonDef = G4Neutron::Neutron();
499 outNucleonDef = G4Neutron::Neutron();
500 break;
501
502 default:
503 sigmaLambdaConversionRate = 0.;
504 }
505
506 if (ranflat >= sigmaLambdaConversionRate) return 0;
507
508 G4ThreeVector dummy(0.,0.,0.);
509
510 // Fermi momentum distribution in three dimensions
511 G4ThreeVector momentum = nucleus->GetFermiMomentum();
512
513 G4ParticleDefinition* lambdaDef = G4Lambda::Lambda();
514
515 G4DynamicParticle inNucleon(inNucleonDef,momentum);
516 G4DynamicParticle outNucleon(outNucleonDef,dummy);
517 G4DynamicParticle* outLambda = new G4DynamicParticle(lambdaDef,dummy);
518
519 G4ReactionKinematics theReactionKinematics;
520
521 // Now do the twobody scattering
522 theReactionKinematics.TwoBodyScattering(aSigma, &inNucleon,
523 &outNucleon, outLambda);
524
525 // Binding energy of nucleus has changed. This will change the
526 // ExcitationEnergy.
527 // .. using mass formula ..
528 // .. using mass table ..
529 // equivalent to -'initialBindingEnergy+nucleus.GetBindingEnergy' !
530
531 // Add energy and momentum to nucleus, change Z,A
532 nucleus->AddExcitationEnergy(outNucleon.GetKineticEnergy());
533 nucleus->AddMomentum(outNucleon.GetMomentum());
534 nucleus->SetParameters(A,newZ);
535
536 // The calling routine is responsible to delete the sigma!!
537 return outLambda;
538}
539
540
541
Note: See TracBrowser for help on using the repository browser.