Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HEKaonMinusInelastic.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
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 "G4HEKaonMinusInelastic.hh"
38 #include "globals.hh"
39 #include "G4ios.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 
43 void G4HEKaonMinusInelastic::ModelDescription(std::ostream& outFile) const
44 {
45  outFile << "G4HEKaonMinusInelastic is one of the High Energy\n"
46  << "Parameterized (HEP) models used to implement inelastic\n"
47  << "K- scattering from nuclei. It is a re-engineered\n"
48  << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
49  << "initial collision products into backward- and forward-going\n"
50  << "clusters which are then decayed into final state hadrons.\n"
51  << "The model does not conserve energy on an event-by-event\n"
52  << "basis. It may be applied to K- with initial energies\n"
53  << "above 20 GeV.\n";
54 }
55 
56 
57 G4HadFinalState*
58 G4HEKaonMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
59  G4Nucleus& targetNucleus)
60 {
61  G4HEVector* pv = new G4HEVector[MAXPART];
62  const G4HadProjectile* aParticle = &aTrack;
63  const G4double A = targetNucleus.GetA_asInt();
64  const G4double Z = targetNucleus.GetZ_asInt();
65  G4HEVector incidentParticle(aParticle);
66 
67  G4double atomicNumber = Z;
68  G4double atomicWeight = A;
69 
70  G4int incidentCode = incidentParticle.getCode();
71  G4double incidentMass = incidentParticle.getMass();
72  G4double incidentTotalEnergy = incidentParticle.getEnergy();
73 
74  // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
75  // DHW 19 May 2011: variable set but not used
76 
77  G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
78 
79  if (incidentKineticEnergy < 1.)
80  G4cout << "GHEKaonMinusInelastic: incident energy < 1 GeV" << G4endl;
81 
82  if (verboseLevel > 1) {
83  G4cout << "G4HEKaonMinusInelastic::ApplyYourself" << G4endl;
84  G4cout << "incident particle " << incidentParticle.getName()
85  << "mass " << incidentMass
86  << "kinetic energy " << incidentKineticEnergy
87  << G4endl;
88  G4cout << "target material with (A,Z) = ("
89  << atomicWeight << "," << atomicNumber << ")" << G4endl;
90  }
91 
92  G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
93  atomicWeight, atomicNumber);
94  if (verboseLevel > 1)
95  G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
96 
97  incidentKineticEnergy -= inelasticity;
98 
99  G4double excitationEnergyGNP = 0.;
100  G4double excitationEnergyDTA = 0.;
101 
102  G4double excitation = NuclearExcitation(incidentKineticEnergy,
103  atomicWeight, atomicNumber,
104  excitationEnergyGNP,
105  excitationEnergyDTA);
106  if (verboseLevel > 1)
107  G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
108  << excitationEnergyDTA << G4endl;
109 
110  incidentKineticEnergy -= excitation;
111  incidentTotalEnergy = incidentKineticEnergy + incidentMass;
112  // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
113  // *(incidentTotalEnergy+incidentMass));
114  // DHW 19 May 2011: variable set but not used
115 
116  G4HEVector targetParticle;
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 
131  vecLength = 0;
132 
133  if (verboseLevel > 1)
134  G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
135  << incidentCode << G4endl;
136 
137  G4bool successful = false;
138 
139  FirstIntInCasKaonMinus(inElastic, availableEnergy, pv, vecLength,
140  incidentParticle, targetParticle);
141 
142  if (verboseLevel > 1)
143  G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
144 
145  if ((vecLength > 0) && (availableEnergy > 1.))
146  StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
147  pv, vecLength,
148  incidentParticle, targetParticle);
149 
150  HighEnergyCascading(successful, pv, vecLength,
151  excitationEnergyGNP, excitationEnergyDTA,
152  incidentParticle, targetParticle,
153  atomicWeight, atomicNumber);
154  if (!successful)
155  HighEnergyClusterProduction(successful, pv, vecLength,
156  excitationEnergyGNP, excitationEnergyDTA,
157  incidentParticle, targetParticle,
158  atomicWeight, atomicNumber);
159  if (!successful)
160  MediumEnergyCascading(successful, pv, vecLength,
161  excitationEnergyGNP, excitationEnergyDTA,
162  incidentParticle, targetParticle,
163  atomicWeight, atomicNumber);
164 
165  if (!successful)
166  MediumEnergyClusterProduction(successful, pv, vecLength,
167  excitationEnergyGNP, excitationEnergyDTA,
168  incidentParticle, targetParticle,
169  atomicWeight, atomicNumber);
170  if (!successful)
171  QuasiElasticScattering(successful, pv, vecLength,
172  excitationEnergyGNP, excitationEnergyDTA,
173  incidentParticle, targetParticle,
174  atomicWeight, atomicNumber);
175  if (!successful)
176  ElasticScattering(successful, pv, vecLength,
177  incidentParticle,
178  atomicWeight, atomicNumber);
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 G4HEKaonMinusInelastic::FirstIntInCasKaonMinus(G4bool& inElastic,
192  const G4double availableEnergy,
193  G4HEVector pv[],
194  G4int& vecLen,
195  const G4HEVector& incidentParticle,
196  const G4HEVector& targetParticle)
197 
198 // Kaon- undergoes interaction with nucleon within a nucleus. Check if it is
199 // energetically possible to produce pions/kaons. In not, assume nuclear excitation
200 // occurs and input particle is degraded in energy. No other particles are produced.
201 // If reaction is possible, find the correct number of pions/protons/neutrons
202 // produced using an interpolation to multiplicity data. Replace some pions or
203 // protons/neutrons by kaons or strange baryons according to the average
204 // multiplicity per inelastic reaction.
205 {
206  static const G4double expxu = 82.; // upper bound for arg. of exp
207  static const G4double expxl = -expxu; // lower bound for arg. of exp
208 
209  static const G4double protb = 0.7;
210  static const G4double neutb = 0.7;
211  static const G4double c = 1.25;
212 
213  static const G4int numMul = 1200;
214  static const G4int numSec = 60;
215 
216  G4int neutronCode = Neutron.getCode();
217  G4int protonCode = Proton.getCode();
218  G4int kaonMinusCode = KaonMinus.getCode();
219  G4int kaonZeroCode = KaonZero.getCode();
220  G4int antiKaonZeroCode = AntiKaonZero.getCode();
221 
222  G4int targetCode = targetParticle.getCode();
223  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
224 
225  static G4bool first = true;
226  static G4double protmul[numMul], protnorm[numSec]; // proton constants
227  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
228 
229  // misc. local variables
230  // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
231 
232  G4int i, counter, nt, npos, nneg, nzero;
233 
234  if (first) {
235  // compute normalization constants, this will only be done once
236  first = false;
237  for (i = 0; i < numMul; i++) protmul[i] = 0.0;
238  for (i = 0; i < numSec; i++) protnorm[i] = 0.0;
239  counter = -1;
240  for (npos = 0; npos < (numSec/3); npos++) {
241  for( nneg=Imax(0,npos-1); nneg<=npos+1; nneg++ )
242  {
243  for( nzero=0; nzero<numSec/3; nzero++ )
244  {
245  if( ++counter < numMul )
246  {
247  nt = npos+nneg+nzero;
248  if( (nt>0) && (nt<=numSec) )
249  {
250  protmul[counter] =
251  pmltpc(npos,nneg,nzero,nt,protb,c) ;
252  protnorm[nt-1] += protmul[counter];
253  }
254  }
255  }
256  }
257  }
258  for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
259  for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
260  counter = -1;
261  for (npos = 0; npos < numSec/3; npos++) {
262  for (nneg = npos; nneg <= (npos+2); nneg++) {
263  for( nzero=0; nzero<numSec/3; nzero++ )
264  {
265  if( ++counter < numMul )
266  {
267  nt = npos+nneg+nzero;
268  if( (nt>0) && (nt<=numSec) )
269  {
270  neutmul[counter] =
271  pmltpc(npos,nneg,nzero,nt,neutb,c);
272  neutnorm[nt-1] += neutmul[counter];
273  }
274  }
275  }
276  }
277  }
278 
279  for (i = 0; i < numSec; i++) {
280  if (protnorm[i] > 0.0) protnorm[i] = 1.0/protnorm[i];
281  if (neutnorm[i] > 0.0) neutnorm[i] = 1.0/neutnorm[i];
282  }
283  } // end of initialization
284 
285  pv[0] = incidentParticle; // initialize the first two places
286  pv[1] = targetParticle; // the same as beam and target
287  vecLen = 2;
288 
289  if (!inElastic || (availableEnergy <= PionPlus.getMass())) return;
290 
291  // inelastic scattering
292  npos = 0, nneg = 0, nzero = 0;
293  G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15};
294  G4int iplab = G4int( incidentTotalMomentum*5.);
295  if ( (iplab < 10) && (G4UniformRand() < cech[iplab])) {
296  G4int ipl = Imin(19, G4int( incidentTotalMomentum*5.));
297  G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45,
298  0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33};
299  if (G4UniformRand() < cnk0[ipl]) {
300  if (targetCode == protonCode) {
301  pv[0] = AntiKaonZero;
302  pv[1] = Neutron;
303  return;
304  } else {
305  return;
306  }
307  }
308  G4double ran = G4UniformRand();
309  if (targetCode == protonCode) { // target is a proton
310  if( ran < 0.25 )
311  {
312  pv[0] = PionMinus;
313  pv[1] = SigmaPlus;
314  }
315  else if (ran < 0.50)
316  {
317  pv[0] = PionZero;
318  pv[1] = SigmaZero;
319  }
320  else if (ran < 0.75)
321  {
322  pv[0] = PionPlus;
323  pv[1] = SigmaMinus;
324  }
325  else
326  {
327  pv[0] = PionZero;
328  pv[1] = Lambda;
329  }
330  } else { // target is a neutron
331  if( ran < 0.25 )
332  {
333  }
334  else if (ran < 0.50)
335  {
336  pv[0] = PionMinus;
337  pv[1] = SigmaZero;
338  }
339  else if (ran < 0.75)
340  {
341  pv[0] = PionZero;
342  pv[1] = SigmaMinus;
343  }
344  else
345  {
346  pv[0] = PionMinus;
347  pv[1] = Lambda;
348  }
349  }
350  return;
351  } else {
352  // number of total particles vs. centre of mass Energy - 2*proton mass
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-1); nneg<=npos+1; 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 
392  if (ran < excs) goto outOfLoop; //----------------------->
393  }
394  }
395  }
396  }
397  }
398 
399  // 3 previous loops continued to the end
400  inElastic = false; // quasi-elastic scattering
401  return;
402  }
403  else
404  { // target must be a neutron
405  counter = -1;
406  for( npos=0; npos<numSec/3; npos++ )
407  {
408  for( nneg=npos; nneg<=(npos+2); nneg++ )
409  {
410  for (nzero=0; nzero<numSec/3; nzero++) {
411  if (++counter < numMul) {
412  nt = npos+nneg+nzero;
413  if ( (nt>=1) && (nt<=numSec) ) {
414  test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
415  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
416  if (std::fabs(dum) < 1.0) {
417  if( test >= 1.0e-10 )excs += dum*test;
418  } else {
419  excs += dum*test;
420  }
421  if (ran < excs) goto outOfLoop; // -------------------------->
422  }
423  }
424  }
425  }
426  }
427  // 3 previous loops continued to the end
428  inElastic = false; // quasi-elastic scattering.
429  return;
430  }
431  }
432  outOfLoop: // <---------------------------------------------
433 
434  if( targetCode == protonCode)
435  {
436  if( npos == (1+nneg))
437  {
438  pv[1] = Neutron;
439  }
440  else if (npos == nneg)
441  {
442  if( G4UniformRand() < 0.75)
443  {
444  }
445  else
446  {
447  pv[0] = AntiKaonZero;
448  pv[1] = Neutron;
449  }
450  }
451  else
452  {
453  pv[0] = AntiKaonZero;
454  }
455  }
456  else
457  {
458  if( npos == (nneg-1))
459  {
460  if( G4UniformRand() < 0.5)
461  {
462  pv[1] = Proton;
463  }
464  else
465  {
466  pv[0] = AntiKaonZero;
467  }
468  }
469  else if ( npos == nneg)
470  {
471  }
472  else
473  {
474  pv[0] = AntiKaonZero;
475  pv[1] = Proton;
476  }
477  }
478 
479  if( G4UniformRand() < 0.5 )
480  {
481  if( ( (pv[0].getCode() == kaonMinusCode)
482  && (pv[1].getCode() == neutronCode) )
483  || ( (pv[0].getCode() == kaonZeroCode)
484  && (pv[1].getCode() == protonCode) )
485  || ( (pv[0].getCode() == antiKaonZeroCode)
486  && (pv[1].getCode() == protonCode) ) )
487  {
488  G4double ran = G4UniformRand();
489  if( pv[1].getCode() == protonCode)
490  {
491  if(ran < 0.68)
492  {
493  pv[0] = PionPlus;
494  pv[1] = Lambda;
495  }
496  else if (ran < 0.84)
497  {
498  pv[0] = PionZero;
499  pv[1] = SigmaPlus;
500  }
501  else
502  {
503  pv[0] = PionPlus;
504  pv[1] = SigmaZero;
505  }
506  }
507  else
508  {
509  if(ran < 0.68)
510  {
511  pv[0] = PionMinus;
512  pv[1] = Lambda;
513  }
514  else if (ran < 0.84)
515  {
516  pv[0] = PionMinus;
517  pv[1] = SigmaZero;
518  }
519  else
520  {
521  pv[0] = PionZero;
522  pv[1] = SigmaMinus;
523  }
524  }
525  }
526  else
527  {
528  G4double ran = G4UniformRand();
529  if (ran < 0.67)
530  {
531  pv[0] = PionZero;
532  pv[1] = Lambda;
533  }
534  else if (ran < 0.78)
535  {
536  pv[0] = PionMinus;
537  pv[1] = SigmaPlus;
538  }
539  else if (ran < 0.89)
540  {
541  pv[0] = PionZero;
542  pv[1] = SigmaZero;
543  }
544  else
545  {
546  pv[0] = PionPlus;
547  pv[1] = SigmaMinus;
548  }
549  }
550  }
551 
552 
553  nt = npos + nneg + nzero;
554  while ( nt > 0)
555  {
556  G4double ran = G4UniformRand();
557  if ( ran < (G4double)npos/nt)
558  {
559  if( npos > 0 )
560  { pv[vecLen++] = PionPlus;
561  npos--;
562  }
563  }
564  else if ( ran < (G4double)(npos+nneg)/nt)
565  {
566  if( nneg > 0 )
567  {
568  pv[vecLen++] = PionMinus;
569  nneg--;
570  }
571  }
572  else
573  {
574  if( nzero > 0 )
575  {
576  pv[vecLen++] = PionZero;
577  nzero--;
578  }
579  }
580  nt = npos + nneg + nzero;
581  }
582  if (verboseLevel > 1)
583  {
584  G4cout << "Particles produced: " ;
585  G4cout << pv[0].getName() << " " ;
586  G4cout << pv[1].getName() << " " ;
587  for (i=2; i < vecLen; i++)
588  {
589  G4cout << pv[i].getName() << " " ;
590  }
591  G4cout << G4endl;
592  }
593  return;
594  }
595 
596 
597 
598 
599 
600 
601 
602 
603 
Generated on Sat May 25 2013 14:33:56 for Geant4 by   doxygen 1.8.4