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

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

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