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