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

Last change on this file since 889 was 819, checked in by garnier, 18 years ago

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