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

Last change on this file since 1194 was 1055, checked in by garnier, 17 years ago

maj sur la beta de geant 4.9.3
File size: 17.7 KB
RevLine 
[819]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"
[1055]45#include "G4HadronicProcessStore.hh"
[819]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 }
[962]54 SetProcessSubType(fHadronAtRest);
[819]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;
[1055]83
84 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
[819]85}
86
87
88G4KaonMinusAbsorptionAtRest::~G4KaonMinusAbsorptionAtRest()
[1055]89{
90 G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
91}
[819]92
[1055]93void G4KaonMinusAbsorptionAtRest::PreparePhysicsTable(const G4ParticleDefinition& p)
94{
95 G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
96}
[819]97
[1055]98void G4KaonMinusAbsorptionAtRest::BuildPhysicsTable(const G4ParticleDefinition& p)
99{
100 G4HadronicProcessStore::Instance()->PrintInfo(&p);
101}
102
[819]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 = pProducts + (*absorptionProducts)[i]->GetMomentum();
193 productEnergy += (*absorptionProducts)[i]->GetKineticEnergy();
194 }
195
196 G4double newZ = nucleus->GetZ();
197 G4double newA = nucleus->GetN();
198
[962]199 G4double bDiff = G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(A),static_cast<G4int>(Z)) -
200 G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(newA), static_cast<G4int>(newZ));
[819]201
202 G4double pNucleus = pProducts.mag();
203
204 G4StopDeexcitationAlgorithm* nucleusAlgorithm = new G4StopTheoDeexcitation();
205 G4StopDeexcitation stopDeexcitation(nucleusAlgorithm);
206
207 nucleus->AddExcitationEnergy(bDiff);
208
209 // returns excitation energy for the moment ..
210 G4double energyDeposit = nucleus->GetEnergyDeposit();
211 if (verboseLevel>0)
212 {
213 G4cout << " -- KaonAtRest -- excitation = "
214 << energyDeposit
215 << ", pNucleus = "
216 << pNucleus
217 << ", A: "
218 << A
219 << ", "
220 << newA
221 << ", Z: "
222 << Z
223 << ", "
224 << newZ
225 << G4endl;
226 }
227
228 if (energyDeposit < 0.)
229 {
230 G4Exception("G4KaonMinusAbsorptionAtRest", "007", FatalException,
231 "AtRestDoIt -- excitation energy < 0");
232 }
233 delete nucleus;
234
235 G4ReactionProductVector* fragmentationProducts = stopDeexcitation.DoBreakUp(newA,newZ,energyDeposit,pNucleus);
236
237 unsigned int nFragmentationProducts = 0;
238 if (fragmentationProducts != 0) nFragmentationProducts = fragmentationProducts->size();
239
240 //Initialize ParticleChange
241 aParticleChange.Initialize(track);
242 aParticleChange.SetNumberOfSecondaries(G4int(nAbsorptionProducts+nFragmentationProducts) );
243
244 // update List of alive particles. put energy deposit at the right place ...
245 for (i = 0; i < nAbsorptionProducts; i++)
246 {aParticleChange.AddSecondary((*absorptionProducts)[i]); }
247 if (absorptionProducts != 0) delete absorptionProducts;
248
249// for (i = 0; i < nFragmentationProducts; i++)
250// { aParticleChange.AddSecondary(fragmentationProducts->at(i)); }
251 for(i=0; i<nFragmentationProducts; i++)
252 {
253 G4DynamicParticle * aNew =
254 new G4DynamicParticle((*fragmentationProducts)[i]->GetDefinition(),
255 (*fragmentationProducts)[i]->GetTotalEnergy(),
256 (*fragmentationProducts)[i]->GetMomentum());
257 G4double newTime = aParticleChange.GetGlobalTime((*fragmentationProducts)[i]->GetFormationTime());
258 aParticleChange.AddSecondary(aNew, newTime);
259 delete (*fragmentationProducts)[i];
260 }
261 if (fragmentationProducts != 0) delete fragmentationProducts;
262
263 // finally ...
264 aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the incident Kaon
265 return &aParticleChange;
266}
267
268
269G4DynamicParticle G4KaonMinusAbsorptionAtRest::GetAbsorbingNucleon()
270{
271 G4DynamicParticle aNucleon;
272
273 // Get nucleon definition, based on Z,N of current Nucleus
274 aNucleon.SetDefinition(SelectAbsorbingNucleon());
275
276 // Fermi momentum distribution in three dimensions
277 G4ThreeVector pFermi = nucleus->GetFermiMomentum();
278 aNucleon.SetMomentum(pFermi);
279
280 return aNucleon;
281}
282
283G4ParticleDefinition* G4KaonMinusAbsorptionAtRest::SelectAbsorbingNucleon()
284{
285 // (Ch. Voelcker) extended from ReturnTargetParticle():
286 // Choose a proton or a neutron as the absorbing particle,
287 // taking weight into account!
288 // Update nucleon's atomic numbers.
289
290 G4ParticleDefinition* absorbingParticleDef;
291
292 G4double ranflat = G4UniformRand();
293
294 G4double myZ = nucleus->GetZ(); // number of protons
295 G4double myN = nucleus->GetN(); // number of nucleons (not neutrons!!)
296
297 // See VanderVelde-Wilquet et al, Nuov.Cim.39A(1978)538;
298 G4double carbonRatioNP = 0.18; // (Rn/Rp)c, see page 544
299
300 G4double neutronProtonRatio = NeutronHaloFactor(myZ,myN)*carbonRatioNP*(myN-myZ)/myZ;
301 G4double protonProbability = 1./(1.+neutronProtonRatio);
302
303 if ( ranflat < protonProbability )
304 {
305 absorbingParticleDef = G4Proton::Proton();
306 myZ-= 1.;
307 }
308 else
309 { absorbingParticleDef = G4Neutron::Neutron(); }
310
311 myN -= 1.;
312 nucleus->SetParameters(myN,myZ);
313 return absorbingParticleDef;
314}
315
316
317G4double G4KaonMinusAbsorptionAtRest::NeutronHaloFactor(G4double Z, G4double N)
318{
319 // this function should take care of the probability for absorption
320 // on neutrons, depending on number of protons Z and number of neutrons N-Z
321 // parametrisation from fit to
322 // VanderVelde-Wilquet et al, Nuov.Cim.39A(1978)538;
323 //
324
325 if (Z == 1.) return 1.389; // deuterium
326 else if (Z == 2.) return 1.78; // helium
327 else if (Z == 10.) return 0.66; // neon
328 else
329 return 0.6742+(N-Z)*0.06524;
330}
331
332
333G4DynamicParticleVector* G4KaonMinusAbsorptionAtRest::KaonNucleonReaction()
334{
335 G4DynamicParticleVector* products = new G4DynamicParticleVector();
336
337 G4double ranflat = G4UniformRand();
338 G4double prob = 0;
339
340 G4ParticleDefinition* producedBaryonDef;
341 G4ParticleDefinition* producedMesonDef;
342
343 G4double iniZ = nucleus->GetZ();
344 G4double iniA = nucleus->GetN();
345
346 G4DynamicParticle aNucleon = GetAbsorbingNucleon();
347
348 G4double nucleonMass;
349
350 if (aNucleon.GetDefinition() == G4Proton::Proton())
351 {
352 nucleonMass = proton_mass_c2+electron_mass_c2;
353 if ( (prob += rateLambdaZeroPiZero) > ranflat)
354 { // lambda pi0
355 producedBaryonDef = G4Lambda::Lambda();
356 producedMesonDef = G4PionZero::PionZero();
357 }
358 else if ((prob += rateSigmaPlusPiMinus) > ranflat)
359 { // sigma+ pi-
360 producedBaryonDef = G4SigmaPlus::SigmaPlus();
361 producedMesonDef = G4PionMinus::PionMinus();
362 }
363 else if ((prob += rateSigmaMinusPiPlus) > ranflat)
364 { // sigma- pi+
365 producedBaryonDef = G4SigmaMinus::SigmaMinus();
366 producedMesonDef = G4PionPlus::PionPlus();
367 }
368 else
369 { // sigma0 pi0
370 producedBaryonDef = G4SigmaZero::SigmaZero();
371 producedMesonDef = G4PionZero::PionZero();
372 }
373 }
374 else if (aNucleon.GetDefinition() == G4Neutron::Neutron())
375 {
376 nucleonMass = neutron_mass_c2;
377 if ((prob += rateLambdaZeroPiMinus) > ranflat)
378 { // lambda pi-
379 producedBaryonDef = G4Lambda::Lambda();
380 producedMesonDef = G4PionMinus::PionMinus();
381 }
382 else if ((prob += rateSigmaZeroPiMinus) > ranflat)
383 { // sigma0 pi-
384 producedBaryonDef = G4SigmaZero::SigmaZero();
385 producedMesonDef = G4PionMinus::PionMinus();
386 }
387 else
388 { // sigma- pi0
389 producedBaryonDef = G4SigmaMinus::SigmaMinus();
390 producedMesonDef = G4PionZero::PionZero();
391 }
392 }
393 else
394 {
395 if (verboseLevel>0)
396 {
397 G4cout
398 << "G4KaonMinusAbsorption::KaonNucleonReaction: "
399 << aNucleon.GetDefinition()->GetParticleName()
400 << " is not a good nucleon - check G4Nucleus::ReturnTargetParticle()!"
401 << G4endl;
402 }
403 return 0;
404 }
405
406 G4double newZ = nucleus->GetZ();
407 G4double newA = nucleus->GetN();
408
409 // Modify the Kaon mass to take nuclear binding energy into account
410 // .. using mas formula ..
411 // .. using mass table ..
412 // equivalent to -'initialBindingEnergy+nucleus.GetBindingEnergy' !
413
414 G4double nucleonBindingEnergy =
[962]415 -G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(iniA), static_cast<G4int>(iniZ) )
416 +G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(newA), static_cast<G4int>(newZ) );
[819]417
418 G4DynamicParticle modifiedHadron = (*stoppedHadron);
419 modifiedHadron.SetMass(stoppedHadron->GetMass() + nucleonBindingEnergy);
420
421 // Setup outgoing dynamic particles
422 G4ThreeVector dummy(0.,0.,0.);
423 G4DynamicParticle* producedBaryon = new G4DynamicParticle(producedBaryonDef,dummy);
424 G4DynamicParticle* producedMeson = new G4DynamicParticle(producedMesonDef,dummy);
425
426 // Produce the secondary particles in a twobody process:
427 G4ReactionKinematics theReactionKinematics;
428 theReactionKinematics.TwoBodyScattering( &modifiedHadron, &aNucleon,
429 producedBaryon, producedMeson);
430
431 products->push_back(producedBaryon);
432 products->push_back(producedMeson);
433
434 if (verboseLevel > 1)
435 {
436 G4cout
437 << "G4KaonMinusAbsorption::KaonNucleonReaction: Number of primaries = "
438 << products->size()
439 << ": " <<producedMesonDef->GetParticleName()
440 << ", " <<producedBaryonDef->GetParticleName() << G4endl;
441 }
442
443 return products;
444}
445
446
447G4bool G4KaonMinusAbsorptionAtRest::AbsorbPionByNucleus(G4DynamicParticle* aPion)
448{
449 // Needs some more investigation!
450
451 G4double ranflat = G4UniformRand();
452
453 if (ranflat < pionAbsorptionRate){
454 // Add pion energy to ExcitationEnergy and NucleusMomentum
455 nucleus->AddExcitationEnergy(aPion->GetTotalEnergy());
456 nucleus->AddMomentum(aPion->GetMomentum());
457 }
458
459 return (ranflat < pionAbsorptionRate);
460}
461
462G4DynamicParticle* G4KaonMinusAbsorptionAtRest::SigmaLambdaConversion(G4DynamicParticle* aSigma)
463{
464 G4double ranflat = G4UniformRand();
465 G4double sigmaLambdaConversionRate;
466
467 G4double A = nucleus->GetN();
468 G4double Z = nucleus->GetZ();
469
470 G4double newZ = Z;
471 G4double nucleonMassDifference = 0;
472
473 G4ParticleDefinition* inNucleonDef=NULL;
474 G4ParticleDefinition* outNucleonDef=NULL;
475
476 // Decide which sigma
477 switch((int) aSigma->GetDefinition()->GetPDGCharge()) {
478
479 case 1:
480 sigmaLambdaConversionRate = sigmaPlusLambdaConversionRate;
481 inNucleonDef = G4Neutron::Neutron();
482 outNucleonDef = G4Proton::Proton();
483 newZ = Z+1;
484 nucleonMassDifference = neutron_mass_c2 - proton_mass_c2-electron_mass_c2;
485 break;
486
487 case -1:
488 sigmaLambdaConversionRate = sigmaMinusLambdaConversionRate;
489 inNucleonDef = G4Proton::Proton();
490 outNucleonDef = G4Neutron::Neutron();
491 newZ = Z-1;
492 nucleonMassDifference = proton_mass_c2+electron_mass_c2 - neutron_mass_c2;
493 break;
494
495 case 0:
496 sigmaLambdaConversionRate = sigmaZeroLambdaConversionRate;
497 // The 'outgoing' nucleon is just virtual, to keep the energy-momentum
498 // balance and will not appear in the ParticleChange. Therefore no need
499 // choose between neutron and proton here!
500 inNucleonDef = G4Neutron::Neutron();
501 outNucleonDef = G4Neutron::Neutron();
502 break;
503
504 default:
505 sigmaLambdaConversionRate = 0.;
506 }
507
508 if (ranflat >= sigmaLambdaConversionRate) return 0;
509
510 G4ThreeVector dummy(0.,0.,0.);
511
512 // Fermi momentum distribution in three dimensions
513 G4ThreeVector momentum = nucleus->GetFermiMomentum();
514
515 G4ParticleDefinition* lambdaDef = G4Lambda::Lambda();
516
517 G4DynamicParticle inNucleon(inNucleonDef,momentum);
518 G4DynamicParticle outNucleon(outNucleonDef,dummy);
519 G4DynamicParticle* outLambda = new G4DynamicParticle(lambdaDef,dummy);
520
521 G4ReactionKinematics theReactionKinematics;
522
523 // Now do the twobody scattering
524 theReactionKinematics.TwoBodyScattering(aSigma, &inNucleon,
525 &outNucleon, outLambda);
526
527 // Binding energy of nucleus has changed. This will change the
528 // ExcitationEnergy.
529 // .. using mass formula ..
530 // .. using mass table ..
531 // equivalent to -'initialBindingEnergy+nucleus.GetBindingEnergy' !
[962]532
[819]533 // Add energy and momentum to nucleus, change Z,A
534 nucleus->AddExcitationEnergy(outNucleon.GetKineticEnergy());
535 nucleus->AddMomentum(outNucleon.GetMomentum());
536 nucleus->SetParameters(A,newZ);
537
538 // The calling routine is responsible to delete the sigma!!
539 return outLambda;
540}
541
542
543
Note: See TracBrowser for help on using the repository browser.