Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LEAntiKaonZeroInelastic.cc
Go to the documentation of this file.
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 // $Id$
27 //
28 // Hadronic Process: Low Energy KaonZeroLong Inelastic Process
29 // J.L. Chuma, TRIUMF, 11-Feb-1997
30 // Modified by J.L.Chuma 30-Apr-97: added originalTarget for CalculateMomenta
31 
32 #include "G4LEAntiKaonZeroInelastic.hh"
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "Randomize.hh"
36 #include "G4HadReentrentException.hh"
37 
38 
39 void G4LEAntiKaonZeroInelastic::ModelDescription(std::ostream& outFile) const
40 {
41  outFile << "G4LEAntiKaonZeroInelastic is one of the Low Energy\n"
42  << "Parameterized (LEP) models used to implement anti-K0 scattering\n"
43  << "from nuclei. It is a re-engineered version of the GHEISHA code\n"
44  << "of H. Fesefeldt. It divides the initial collision products\n"
45  << "into backward- and forward-going clusters which are then\n"
46  << "decayed into final state hadrons. The model does not conserve\n"
47  << "energy on an event-by-event basis. It may be applied to\n"
48  << "anti-K0s with initial energies between 0 and 25 GeV.\n";
49 }
50 
51 
52 G4HadFinalState*
53 G4LEAntiKaonZeroInelastic::ApplyYourself(const G4HadProjectile& aTrack,
54  G4Nucleus& targetNucleus)
55 {
56  const G4HadProjectile* originalIncident = &aTrack;
57 
58  // create the target particle
59  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
60 
61  if (verboseLevel > 1) {
62  const G4Material *targetMaterial = aTrack.GetMaterial();
63  G4cout << "G4LEAntiKaonZeroInelastic::ApplyYourself called" << G4endl;
64  G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
65  G4cout << "target material = " << targetMaterial->GetName() << ", ";
66  G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
67  << G4endl;
68  }
69 
70  // Fermi motion and evaporation
71  // As of Geant3, the Fermi energy calculation had not been Done
72  G4double ek = originalIncident->GetKineticEnergy()/MeV;
73  G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
74  G4ReactionProduct modifiedOriginal;
75  modifiedOriginal = *originalIncident;
76 
77  G4double tkin = targetNucleus.Cinema( ek );
78  ek += tkin;
79  modifiedOriginal.SetKineticEnergy( ek*MeV );
80  G4double et = ek + amas;
81  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
82  G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
83  if (pp > 0.0) {
84  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
85  modifiedOriginal.SetMomentum( momentum * (p/pp) );
86  }
87 
88  // calculate black track energies
89  tkin = targetNucleus.EvaporationEffects( ek );
90  ek -= tkin;
91  modifiedOriginal.SetKineticEnergy( ek*MeV );
92  et = ek + amas;
93  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
94  pp = modifiedOriginal.GetMomentum().mag()/MeV;
95  if (pp > 0.0) {
96  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
97  modifiedOriginal.SetMomentum( momentum * (p/pp) );
98  }
99 
100  G4ReactionProduct currentParticle = modifiedOriginal;
101  G4ReactionProduct targetParticle;
102  targetParticle = *originalTarget;
103  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
104  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
105  G4bool incidentHasChanged = false;
106  G4bool targetHasChanged = false;
107  G4bool quasiElastic = false;
108  G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec; // vec will contain the secondary particles
109  G4int vecLen = 0;
110  vec.Initialize( 0 );
111 
112  const G4double cutOff = 0.1;
113  if (currentParticle.GetKineticEnergy()/MeV > cutOff)
114  Cascade(vec, vecLen,
115  originalIncident, currentParticle, targetParticle,
116  incidentHasChanged, targetHasChanged, quasiElastic);
117 
118  try {
119  CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
120  modifiedOriginal, targetNucleus, currentParticle,
121  targetParticle, incidentHasChanged, targetHasChanged,
122  quasiElastic);
123  }
124  catch(G4HadReentrentException aR)
125  {
126  aR.Report(G4cout);
127  throw G4HadReentrentException(__FILE__, __LINE__, "Bailing out");
128  }
129  SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
130 
131  if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
132  delete originalTarget;
133  return &theParticleChange;
134 }
135 
136 
137 void G4LEAntiKaonZeroInelastic::Cascade(
138  G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
139  G4int& vecLen,
140  const G4HadProjectile *originalIncident,
141  G4ReactionProduct &currentParticle,
142  G4ReactionProduct &targetParticle,
143  G4bool &incidentHasChanged,
144  G4bool &targetHasChanged,
145  G4bool &quasiElastic)
146 {
147  // derived from original FORTRAN code CASK0B by H. Fesefeldt (13-Sep-1987)
148  //
149  // K0Long undergoes interaction with nucleon within a nucleus. Check if it is
150  // energetically possible to produce pions/kaons. In not, assume nuclear excitation
151  // occurs and input particle is degraded in energy. No other particles are produced.
152  // If reaction is possible, find the correct number of pions/protons/neutrons
153  // produced using an interpolation to multiplicity data. Replace some pions or
154  // protons/neutrons by kaons or strange baryons according to the average
155  // multiplicity per Inelastic reaction.
156 
157  const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
158  const G4double etOriginal = originalIncident->Get4Momentum().e()/MeV;
159  const G4double pOriginal = originalIncident->Get4Momentum().vect().mag()/MeV;
160  const G4double targetMass = targetParticle.GetMass()/MeV;
161  G4double centerofmassEnergy = std::sqrt(mOriginal*mOriginal +
162  targetMass*targetMass +
163  2.0*targetMass*etOriginal );
164  G4double availableEnergy = centerofmassEnergy - (targetMass+mOriginal);
165 
166  static G4bool first = true;
167  const G4int numMul = 1200;
168  const G4int numSec = 60;
169  static G4double protmul[numMul], protnorm[numSec]; // proton constants
170  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
171 
172  // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
173  G4int counter;
174  G4int nt = 0;
175  G4int npos = 0, nneg = 0, nzero = 0;
176  const G4double c = 1.25;
177  const G4double b[] = { 0.7, 0.7 };
178  if (first) { // Computation of normalization constants will only be done once
179  first = false;
180  G4int i;
181  for (i = 0; i < numMul; ++i) protmul[i] = 0.0;
182  for (i = 0; i < numSec; ++i) protnorm[i] = 0.0;
183  counter = -1;
184  for (npos = 0; npos < numSec/3; ++npos) {
185  for (nneg = std::max(0,npos - 2); nneg <= npos; ++nneg) {
186  for (nzero = 0; nzero < numSec/3; ++nzero) {
187  if (++counter < numMul) {
188  nt = npos + nneg + nzero;
189  if (nt > 0 && nt <= numSec) {
190  protmul[counter] = Pmltpc(npos, nneg, nzero, nt, b[0], c);
191  protnorm[nt-1] += protmul[counter];
192  }
193  }
194  }
195  }
196  }
197 
198  for (i = 0; i < numMul; ++i) neutmul[i] = 0.0;
199  for (i = 0; i < numSec; ++i) neutnorm[i] = 0.0;
200  counter = -1;
201  for (npos = 0; npos < (numSec/3); ++npos) {
202  for (nneg = std::max(0,npos-1); nneg <= (npos+1); ++nneg) {
203  for (nzero = 0; nzero < numSec/3; ++nzero) {
204  if (++counter < numMul) {
205  nt = npos + nneg + nzero;
206  if (nt > 0 && nt <= numSec) {
207  neutmul[counter] = Pmltpc(npos, nneg, nzero, nt, b[1], c);
208  neutnorm[nt-1] += neutmul[counter];
209  }
210  }
211  }
212  }
213  }
214 
215  for (i = 0; i < numSec; ++i) {
216  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
217  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
218  }
219  } // end of initialization
220 
221  const G4double expxu = 82.; // upper bound for arg. of exp
222  const G4double expxl = -expxu; // lower bound for arg. of exp
223  G4ParticleDefinition* aKaonMinus = G4KaonMinus::KaonMinus();
224  G4ParticleDefinition* aKaonZS = G4KaonZeroShort::KaonZeroShort();
225  G4ParticleDefinition* aKaonZL = G4KaonZeroLong::KaonZeroLong();
226  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
227  G4ParticleDefinition *aProton = G4Proton::Proton();
228  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
229  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
230  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
231  G4ParticleDefinition *aLambda = G4Lambda::Lambda();
232  G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
233  G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
234  G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
235  const G4double cech[] = {1.,1.,1.,0.70,0.60,0.55,0.35,0.25,0.18,0.15};
236  G4int iplab = G4int(std::min( 9.0, 5.0*pOriginal*MeV/GeV ));
237  if ((pOriginal*MeV/GeV <= 2.0) && (G4UniformRand() < cech[iplab]) ) {
238  npos = nneg = nzero = nt = 0;
239  iplab = G4int(std::min( 19.0, pOriginal*MeV/GeV*10.0 ));
240  const G4double cnk0[] = {0.17,0.18,0.17,0.24,0.26,0.20,0.22,0.21,0.34,0.45,
241  0.58,0.55,0.36,0.29,0.29,0.32,0.32,0.33,0.33,0.33};
242  if (G4UniformRand() > cnk0[iplab]) {
243  G4double ran = G4UniformRand();
244  if (ran < 0.25) {
245  // k0Long n --> pi- s+
246  if (targetParticle.GetDefinition() == aNeutron) {
247  currentParticle.SetDefinitionAndUpdateE(aPiMinus);
248  targetParticle.SetDefinitionAndUpdateE(aSigmaPlus);
249  incidentHasChanged = true;
250  targetHasChanged = true;
251  }
252  } else if( ran < 0.50 ) {
253  // k0Long p --> pi+ s0 or k0Long n --> pi0 s0
254  if( targetParticle.GetDefinition() == aNeutron )
255  currentParticle.SetDefinitionAndUpdateE( aPiZero );
256  else
257  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
258  targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
259  incidentHasChanged = true;
260  targetHasChanged = true;
261  } else if (ran < 0.75) {
262  // k0Long n --> pi+ s-
263  if( targetParticle.GetDefinition() == aNeutron )
264  {
265  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
266  targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
267  incidentHasChanged = true;
268  targetHasChanged = true;
269  }
270  } else {
271  // k0Long p --> pi+ L or k0Long n --> pi0 L
272  if( targetParticle.GetDefinition() == aNeutron )
273  currentParticle.SetDefinitionAndUpdateE( aPiZero );
274  else
275  currentParticle.SetDefinitionAndUpdateE( aPiPlus );
276  targetParticle.SetDefinitionAndUpdateE( aLambda );
277  incidentHasChanged = true;
278  targetHasChanged = true;
279  }
280  } else {
281  // ran <= cnk0
282  quasiElastic = true;
283  if (targetParticle.GetDefinition() == aNeutron) {
284  currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
285  targetParticle.SetDefinitionAndUpdateE( aProton );
286  incidentHasChanged = true;
287  targetHasChanged = true;
288  }
289  }
290  } else {
291  // (pOriginal > 2.0*GeV) || (random number >= cech[iplab])
292  if (availableEnergy < aPiPlus->GetPDGMass()/MeV) {
293  quasiElastic = true;
294  return;
295  }
296  G4double n, anpn;
297  GetNormalizationConstant( availableEnergy, n, anpn );
298  G4double ran = G4UniformRand();
299  G4double dum, test, excs = 0.0;
300  if (targetParticle.GetDefinition() == aProton) {
301  counter = -1;
302  for (npos = 0; (npos < numSec/3) && (ran >= excs); ++npos) {
303  for (nneg = std::max(0,npos-2); nneg <= npos && ran >= excs; ++nneg) {
304  for (nzero = 0; nzero < numSec/3 && ran >= excs; ++nzero) {
305  if (++counter < numMul) {
306  nt = npos + nneg +nzero;
307  if (nt > 0 && nt <= numSec) {
308  test = std::exp(std::min(expxu,
309  std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
310  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
311  if (std::fabs(dum) < 1.0) {
312  if( test >= 1.0e-10 )excs += dum*test;
313  } else {
314  excs += dum*test;
315  }
316  }
317  }
318  }
319  }
320  }
321  if (ran >= excs) {
322  // 3 previous loops continued to the end
323  quasiElastic = true;
324  return;
325  }
326  npos--; nneg--; nzero--;
327  switch (npos - nneg)
328  {
329  case 1:
330  if (G4UniformRand() < 0.5) {
331  currentParticle.SetDefinitionAndUpdateE(aKaonMinus);
332  incidentHasChanged = true;
333  } else {
334  targetParticle.SetDefinitionAndUpdateE(aNeutron);
335  targetHasChanged = true;
336  }
337  case 0:
338  break;
339  default:
340  currentParticle.SetDefinitionAndUpdateE(aKaonMinus);
341  targetParticle.SetDefinitionAndUpdateE(aNeutron);
342  incidentHasChanged = true;
343  targetHasChanged = true;
344  break;
345  }
346  } else {
347  // target must be a neutron
348  counter = -1;
349  for (npos = 0; npos < numSec/3 && ran >= excs; ++npos) {
350  for (nneg = std::max(0,npos-1); nneg <= (npos+1) && ran >= excs; ++nneg) {
351  for (nzero = 0; nzero < numSec/3 && ran >= excs; ++nzero) {
352  if (++counter < numMul) {
353  nt = npos + nneg + nzero;
354  if (nt > 0 && nt <= numSec) {
355  test = std::exp(std::min(expxu,
356  std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
357  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
358  if (std::fabs(dum) < 1.0) {
359  if (test >= 1.0e-10) excs += dum*test;
360  } else {
361  excs += dum*test;
362  }
363  }
364  }
365  }
366  }
367  }
368  if (ran >= excs) {
369  // 3 previous loops continued to the end
370  quasiElastic = true;
371  return;
372  }
373  npos--; nneg--; nzero--;
374  switch (npos - nneg)
375  {
376  case 0:
377  currentParticle.SetDefinitionAndUpdateE(aKaonMinus);
378  targetParticle.SetDefinitionAndUpdateE(aProton);
379  incidentHasChanged = true;
380  targetHasChanged = true;
381  break;
382  case 1:
383  currentParticle.SetDefinitionAndUpdateE(aKaonMinus);
384  incidentHasChanged = true;
385  break;
386  default:
387  targetParticle.SetDefinitionAndUpdateE(aProton);
388  targetHasChanged = true;
389  break;
390  }
391  }
392 
393  if (G4UniformRand() >= 0.5) {
394  if (currentParticle.GetDefinition() == aKaonMinus &&
395  targetParticle.GetDefinition() == aNeutron) {
396  ran = G4UniformRand();
397  if (ran < 0.68) {
398  currentParticle.SetDefinitionAndUpdateE(aPiMinus);
399  targetParticle.SetDefinitionAndUpdateE(aLambda);
400  } else if (ran < 0.84) {
401  currentParticle.SetDefinitionAndUpdateE(aPiMinus);
402  targetParticle.SetDefinitionAndUpdateE(aSigmaZero);
403  } else {
404  currentParticle.SetDefinitionAndUpdateE(aPiZero);
405  targetParticle.SetDefinitionAndUpdateE(aSigmaMinus);
406  }
407  } else if ((currentParticle.GetDefinition() == aKaonZS ||
408  currentParticle.GetDefinition() == aKaonZL) &&
409  targetParticle.GetDefinition() == aProton) {
410  ran = G4UniformRand();
411  if (ran < 0.68) {
412  currentParticle.SetDefinitionAndUpdateE(aPiPlus);
413  targetParticle.SetDefinitionAndUpdateE(aLambda);
414  } else if (ran < 0.84) {
415  currentParticle.SetDefinitionAndUpdateE(aPiZero);
416  targetParticle.SetDefinitionAndUpdateE(aSigmaPlus);
417  } else {
418  currentParticle.SetDefinitionAndUpdateE(aPiPlus);
419  targetParticle.SetDefinitionAndUpdateE(aSigmaZero);
420  }
421  } else {
422  ran = G4UniformRand();
423  if (ran < 0.67) {
424  currentParticle.SetDefinitionAndUpdateE(aPiZero);
425  targetParticle.SetDefinitionAndUpdateE(aLambda);
426  } else if (ran < 0.78) {
427  currentParticle.SetDefinitionAndUpdateE(aPiMinus);
428  targetParticle.SetDefinitionAndUpdateE(aSigmaPlus);
429  } else if (ran < 0.89) {
430  currentParticle.SetDefinitionAndUpdateE(aPiZero);
431  targetParticle.SetDefinitionAndUpdateE(aSigmaZero);
432  } else {
433  currentParticle.SetDefinitionAndUpdateE(aPiPlus);
434  targetParticle.SetDefinitionAndUpdateE(aSigmaMinus);
435  }
436  }
437  incidentHasChanged = true;
438  targetHasChanged = true;
439  }
440  }
441 
442  if (currentParticle.GetDefinition() == aKaonZL) {
443  if (G4UniformRand() >= 0.5) {
444  currentParticle.SetDefinitionAndUpdateE(aKaonZS);
445  incidentHasChanged = true;
446  }
447  }
448  if (targetParticle.GetDefinition() == aKaonZL) {
449  if (G4UniformRand() >= 0.5) {
450  targetParticle.SetDefinitionAndUpdateE(aKaonZS);
451  targetHasChanged = true;
452  }
453  }
454  SetUpPions(npos, nneg, nzero, vec, vecLen);
455  return;
456 }
457 
458  /* end of file */
459 
Generated on Sat May 25 2013 14:34:01 for Geant4 by   doxygen 1.8.4