Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HEAntiKaonZeroInelastic.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 
29 // G4 Process: Gheisha High Energy Collision model.
30 // This includes the high energy cascading model, the two-body-resonance model
31 // and the low energy two-body model. Not included is the low energy stuff like
32 // nuclear reactions, nuclear fission without any cascading and all processes
33 // for particles at rest.
34 // First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96.
35 // H. Fesefeldt, RWTH-Aachen, 23-October-1996
36 //
37 
39 #include "globals.hh"
40 #include "G4ios.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 
45  : G4HEInelastic(name)
46 {
47  vecLength = 0;
48  theMinEnergy = 20*GeV;
49  theMaxEnergy = 10*TeV;
50  MAXPART = 2048;
51  verboseLevel = 0;
52  G4cout << "WARNING: model G4HEAntiKaonZeroInelastic is being deprecated and will\n"
53  << "disappear in Geant4 version 10.0" << G4endl;
54 }
55 
56 
58 {
59  outFile << "G4HEAntiKaonZeroInelastic is one of the High Energy\n"
60  << "Parameterized (HEP) models used to implement inelastic\n"
61  << "anti-K0 scattering from nuclei. It is a re-engineered\n"
62  << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
63  << "initial collision products into backward- and forward-going\n"
64  << "clusters which are then decayed into final state hadrons.\n"
65  << "The model does not conserve energy on an event-by-event\n"
66  << "basis. It may be applied to anti-K0 with initial energies\n"
67  << "above 20 GeV.\n";
68 }
69 
70 
73  G4Nucleus& targetNucleus)
74 {
75  G4HEVector* pv = new G4HEVector[MAXPART];
76  const G4HadProjectile* aParticle = &aTrack;
77  const G4double atomicWeight = targetNucleus.GetA_asInt();
78  const G4double atomicNumber = targetNucleus.GetZ_asInt();
79  G4HEVector incidentParticle(aParticle);
80 
81  G4int incidentCode = incidentParticle.getCode();
82  G4double incidentMass = incidentParticle.getMass();
83  G4double incidentTotalEnergy = incidentParticle.getEnergy();
84 
85  // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
86  // DHW 19 May 2011: variable set but not used
87 
88  G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
89 
90  if (incidentKineticEnergy < 1.)
91  G4cout << "GHEAntiKaonZeroInelastic: incident energy < 1 GeV" << G4endl;
92 
93  if (verboseLevel > 1) {
94  G4cout << "G4HEAntiKaonZeroInelastic::ApplyYourself" << G4endl;
95  G4cout << "incident particle " << incidentParticle.getName()
96  << "mass " << incidentMass
97  << "kinetic energy " << incidentKineticEnergy
98  << G4endl;
99  G4cout << "target material with (A,Z) = ("
100  << atomicWeight << "," << atomicNumber << ")" << G4endl;
101  }
102 
103  G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
104  atomicWeight, atomicNumber);
105  if (verboseLevel > 1)
106  G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
107 
108  incidentKineticEnergy -= inelasticity;
109 
110  G4double excitationEnergyGNP = 0.;
111  G4double excitationEnergyDTA = 0.;
112 
113  G4double excitation = NuclearExcitation(incidentKineticEnergy,
114  atomicWeight, atomicNumber,
115  excitationEnergyGNP,
116  excitationEnergyDTA);
117  if (verboseLevel > 1)
118  G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
119  << excitationEnergyDTA << G4endl;
120 
121 
122  incidentKineticEnergy -= excitation;
123  incidentTotalEnergy = incidentKineticEnergy + incidentMass;
124  // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
125  // *(incidentTotalEnergy+incidentMass));
126  // DHW 19 May 2011: variable set but not used
127 
128  G4HEVector targetParticle;
129  if (G4UniformRand() < atomicNumber/atomicWeight) {
130  targetParticle.setDefinition("Proton");
131  } else {
132  targetParticle.setDefinition("Neutron");
133  }
134 
135  G4double targetMass = targetParticle.getMass();
136  G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
137  + targetMass*targetMass
138  + 2.0*targetMass*incidentTotalEnergy);
139  G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
140 
141  vecLength = 0;
142 
143  if (verboseLevel > 1)
144  G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
145  << incidentCode << G4endl;
146 
147  G4bool successful = false;
148 
149  G4bool inElastic = true;
150  FirstIntInCasAntiKaonZero(inElastic, availableEnergy, pv, vecLength,
151  incidentParticle, targetParticle );
152 
153  if (verboseLevel > 1)
154  G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
155 
156  if ((vecLength > 0) && (availableEnergy > 1.))
157  StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
158  pv, vecLength,
159  incidentParticle, targetParticle);
160  HighEnergyCascading(successful, pv, vecLength,
161  excitationEnergyGNP, excitationEnergyDTA,
162  incidentParticle, targetParticle,
163  atomicWeight, atomicNumber);
164  if (!successful)
165  HighEnergyClusterProduction(successful, pv, vecLength,
166  excitationEnergyGNP, excitationEnergyDTA,
167  incidentParticle, targetParticle,
168  atomicWeight, atomicNumber);
169  if (!successful)
170  MediumEnergyCascading(successful, pv, vecLength,
171  excitationEnergyGNP, excitationEnergyDTA,
172  incidentParticle, targetParticle,
173  atomicWeight, atomicNumber);
174  if (!successful)
176  excitationEnergyGNP, excitationEnergyDTA,
177  incidentParticle, targetParticle,
178  atomicWeight, atomicNumber);
179  if (!successful)
180  QuasiElasticScattering(successful, pv, vecLength,
181  excitationEnergyGNP, excitationEnergyDTA,
182  incidentParticle, targetParticle,
183  atomicWeight, atomicNumber);
184  if (!successful)
185  ElasticScattering(successful, pv, vecLength,
186  incidentParticle,
187  atomicWeight, atomicNumber);
188 
189  if (!successful)
190  G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
191  << G4endl;
192 
194  delete [] pv;
196  return &theParticleChange;
197 }
198 
199 
200 void
202  const G4double availableEnergy,
203  G4HEVector pv[],
204  G4int &vecLen,
205  const G4HEVector& incidentParticle,
206  const G4HEVector& targetParticle)
207 
208 // AntiKaon0 undergoes interaction with nucleon within a nucleus. Check if it
209 // is energetically possible to produce pions/kaons. In not, assume nuclear
210 // excitation occurs and input particle is degraded in energy. No other
211 // particles are produced.
212 // If reaction is possible, find the correct number of pions/protons/neutrons
213 // produced using an interpolation to multiplicity data. Replace some pions or
214 // protons/neutrons by kaons or strange baryons according to the average
215 // multiplicity per inelastic reaction.
216 {
217  static const G4double expxu = 82.; // upper bound for arg. of exp
218  static const G4double expxl = -expxu; // lower bound for arg. of exp
219 
220  static const G4double protb = 0.7;
221  static const G4double neutb = 0.7;
222  static const G4double c = 1.25;
223 
224  static const G4int numMul = 1200;
225  static const G4int numSec = 60;
226 
228  G4int protonCode = Proton.getCode();
229  G4int kaonMinusCode = KaonMinus.getCode();
230  G4int kaonZeroCode = KaonZero.getCode();
231  G4int antiKaonZeroCode = AntiKaonZero.getCode();
232 
233  G4int targetCode = targetParticle.getCode();
234  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
235 
236  static G4bool first = true;
237  static G4double protmul[numMul], protnorm[numSec]; // proton constants
238  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
239 
240  // misc. local variables
241  // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
242 
243  G4int i, counter, nt, npos, nneg, nzero;
244 
245  if (first) { // Computation of normalization constants will only be done once
246  first = false;
247  for (i = 0; i < numMul; i++) protmul[i] = 0.0;
248  for (i = 0; i < numSec; i++) protnorm[i] = 0.0;
249  counter = -1;
250  for (npos = 0; npos < (numSec/3); npos++) {
251  for (nneg = std::max(0,npos-2); nneg <= npos; nneg++) {
252  for (nzero = 0; nzero < numSec/3; nzero++) {
253  if (++counter < numMul) {
254  nt = npos+nneg+nzero;
255  if ((nt > 0) && (nt <= numSec) ) {
256  protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
257  protnorm[nt-1] += protmul[counter];
258  }
259  }
260  }
261  }
262  }
263 
264  for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
265  for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
266  counter = -1;
267  for (npos = 0; npos < numSec/3; npos++) {
268  for (nneg = std::max(0,npos-1); nneg <= (npos+1); nneg++) {
269  for (nzero = 0; nzero < numSec/3; nzero++) {
270  if (++counter < numMul) {
271  nt = npos+nneg+nzero;
272  if ((nt > 0) && (nt <= numSec) ) {
273  neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
274  neutnorm[nt-1] += neutmul[counter];
275  }
276  }
277  }
278  }
279  }
280  for (i = 0; i < numSec; i++) {
281  if (protnorm[i] > 0.0) protnorm[i] = 1.0/protnorm[i];
282  if (neutnorm[i] > 0.0) neutnorm[i] = 1.0/neutnorm[i];
283  }
284  } // end of initialization
285 
286  pv[0] = incidentParticle; // initialize the first two places
287  pv[1] = targetParticle; // the same as beam and target
288  vecLen = 2;
289 
290  if (!inElastic || (availableEnergy <= PionPlus.getMass())) return;
291 
292  // inelastic scattering
293  npos = 0, nneg = 0, nzero = 0;
294  G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15};
295  G4int iplab = G4int( incidentTotalMomentum*5.);
296  if ( (iplab < 10) && (G4UniformRand() < cech[iplab]) ) {
297  G4int ipl = std::min(19, G4int( incidentTotalMomentum*5.));
298  G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45,
299  0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33};
300  if (G4UniformRand() < cnk0[ipl] ) {
301  if (targetCode == protonCode) {
302  return;
303  } else {
304  pv[0] = KaonMinus;
305  pv[1] = Proton;
306  return;
307  }
308  }
309 
310  G4double ran = G4UniformRand();
311  if (targetCode == protonCode) { // target is a proton
312  if (ran < 0.25) {
313  } else if (ran < 0.50) {
314  pv[0] = PionPlus;
315  pv[1] = SigmaZero;
316  } else if (ran < 0.75) {
317  } else {
318  pv[0] = PionPlus;
319  pv[1] = Lambda;
320  }
321  } else { // target is a neutron
322  if( ran < 0.25 )
323  {
324  pv[0] = PionMinus;
325  pv[1] = SigmaPlus;
326  }
327  else if (ran < 0.50)
328  {
329  pv[0] = PionZero;
330  pv[1] = SigmaZero;
331  }
332  else if (ran < 0.75)
333  {
334  pv[0] = PionPlus;
335  pv[1] = SigmaMinus;
336  }
337  else
338  {
339  pv[0] = PionZero;
340  pv[1] = Lambda;
341  }
342  }
343  return;
344 
345  } else {
346  // number of total particles vs. centre of mass Energy - 2*proton mass
347  G4double aleab = std::log(availableEnergy);
348  G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
349  + aleab*(0.117712+0.0136912*aleab))) - 2.0;
350 
351  // normalization constant for kno-distribution.
352  // calculate first the sum of all constants, check for numerical problems.
353  G4double test, dum, anpn = 0.0;
354 
355  for (nt=1; nt<=numSec; nt++) {
356  test = std::exp(std::min(expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
357  dum = pi*nt/(2.0*n*n);
358  if (std::fabs(dum) < 1.0) {
359  if( test >= 1.0e-10 )anpn += dum*test;
360  } else {
361  anpn += dum*test;
362  }
363  }
364 
365  G4double ran = G4UniformRand();
366  G4double excs = 0.0;
367  if (targetCode == protonCode) {
368  counter = -1;
369  for (npos=0; npos<numSec/3; npos++) {
370  for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
371  for (nzero=0; nzero<numSec/3; nzero++) {
372  if (++counter < numMul) {
373  nt = npos+nneg+nzero;
374  if ((nt>0) && (nt<=numSec) ) {
375  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
376  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
377 
378  if (std::fabs(dum) < 1.0) {
379  if( test >= 1.0e-10 )excs += dum*test;
380  } else {
381  excs += dum*test;
382  }
383 
384  if (ran < excs) goto outOfLoop; //----------------------->
385  }
386  }
387  }
388  }
389  }
390  // 3 previous loops continued to the end
391  inElastic = false; // quasi-elastic scattering
392  return;
393 
394  } else { // target must be a neutron
395  counter = -1;
396  for (npos=0; npos<numSec/3; npos++) {
397  for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
398  for (nzero=0; nzero<numSec/3; nzero++) {
399  if (++counter < numMul) {
400  nt = npos+nneg+nzero;
401  if( (nt>=1) && (nt<=numSec) ) {
402  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
403  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
404 
405  if (std::fabs(dum) < 1.0) {
406  if( test >= 1.0e-10 )excs += dum*test;
407  } else {
408  excs += dum*test;
409  }
410 
411  if (ran < excs) goto outOfLoop; // -------------------------->
412  }
413  }
414  }
415  }
416  }
417  // 3 previous loops continued to the end
418  inElastic = false; // quasi-elastic scattering.
419  return;
420  }
421  } // if (iplab < 10 .... )
422 
423  outOfLoop: // <------------------------------------------------------------------------
424 
425  if (targetCode == protonCode) {
426  if (npos == nneg) {
427 
428  } else if (npos == (1+nneg)) {
429  if (G4UniformRand() < 0.5) {
430  pv[0] = KaonMinus;
431  } else {
432  pv[1] = Neutron;
433  }
434  } else {
435  pv[0] = KaonMinus;
436  pv[1] = Neutron;
437  }
438 
439  } else {
440  if( npos == nneg)
441  {
442  if( G4UniformRand() < 0.75)
443  {
444  }
445  else
446  {
447  pv[0] = KaonMinus;
448  pv[1] = Proton;
449  }
450  }
451  else if ( npos == (1+nneg))
452  {
453  pv[0] = KaonMinus;
454  }
455  else
456  {
457  pv[1] = Proton;
458  }
459  }
460 
461  if (G4UniformRand() < 0.5) {
462  if (((pv[0].getCode() == kaonMinusCode)
463  && (pv[1].getCode() == neutronCode) )
464  || ((pv[0].getCode() == kaonZeroCode)
465  && (pv[1].getCode() == protonCode) )
466  || ((pv[0].getCode() == antiKaonZeroCode)
467  && (pv[1].getCode() == protonCode) ) ) {
468 
469  G4double ran = G4UniformRand();
470  if (pv[1].getCode() == protonCode) {
471  if (ran < 0.68) {
472  pv[0] = PionPlus;
473  pv[1] = Lambda;
474  } else if (ran < 0.84) {
475  pv[0] = PionZero;
476  pv[1] = SigmaPlus;
477  } else {
478  pv[0] = PionPlus;
479  pv[1] = SigmaZero;
480  }
481  } else {
482  if(ran < 0.68)
483  {
484  pv[0] = PionMinus;
485  pv[1] = Lambda;
486  }
487  else if (ran < 0.84)
488  {
489  pv[0] = PionMinus;
490  pv[1] = SigmaZero;
491  }
492  else
493  {
494  pv[0] = PionZero;
495  pv[1] = SigmaMinus;
496  }
497  }
498  } else {
499  G4double ran = G4UniformRand();
500  if (ran < 0.67)
501  {
502  pv[0] = PionZero;
503  pv[1] = Lambda;
504  }
505  else if (ran < 0.78)
506  {
507  pv[0] = PionMinus;
508  pv[1] = SigmaPlus;
509  }
510  else if (ran < 0.89)
511  {
512  pv[0] = PionZero;
513  pv[1] = SigmaZero;
514  }
515  else
516  {
517  pv[0] = PionPlus;
518  pv[1] = SigmaMinus;
519  }
520  }
521  } // if rand < 0.5
522 
523  nt = npos + nneg + nzero;
524  while (nt > 0) {
525  G4double ran = G4UniformRand();
526  if (ran < (G4double)npos/nt) {
527  if (npos > 0) {
528  pv[vecLen++] = PionPlus;
529  npos--;
530  }
531  } else if (ran < (G4double)(npos+nneg)/nt) {
532  if (nneg > 0) {
533  pv[vecLen++] = PionMinus;
534  nneg--;
535  }
536  } else {
537  if (nzero > 0) {
538  pv[vecLen++] = PionZero;
539  nzero--;
540  }
541  }
542  nt = npos + nneg + nzero;
543  }
544 
545  if (verboseLevel > 1) {
546  G4cout << "Particles produced: " ;
547  G4cout << pv[0].getName() << " " ;
548  G4cout << pv[1].getName() << " " ;
549  for (i=2; i < vecLen; i++) {
550  G4cout << pv[i].getName() << " " ;
551  }
552  G4cout << G4endl;
553  }
554 
555  return;
556  }
557