Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HEKaonZeroShortInelastic.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 //
35 // New version by D.H. Wright (SLAC) to fix seg fault in old version
36 // 21 January 2010
37 
38 #include "G4HEKaonZeroShortInelastic.hh"
39 #include "globals.hh"
40 #include "G4ios.hh"
41 #include "G4PhysicalConstants.hh"
42 
43 void G4HEKaonZeroShortInelastic::ModelDescription(std::ostream& outFile) const
44 {
45  outFile << "G4HEKaonZeroShortInelastic is one of the High Energy\n"
46  << "Parameterized (HEP) models used to implement inelastic\n"
47  << "K0S 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 K0S with initial energies\n"
53  << "above 20 GeV.\n";
54 }
55 
56 
57 G4HadFinalState*
58 G4HEKaonZeroShortInelastic::ApplyYourself(const G4HadProjectile& aTrack,
59  G4Nucleus& targetNucleus)
60 {
61  G4HEVector* pv = new G4HEVector[MAXPART];
62  const G4HadProjectile* aParticle = &aTrack;
63  const G4double atomicWeight = targetNucleus.GetA_asInt();
64  const G4double atomicNumber = targetNucleus.GetZ_asInt();
65  G4HEVector incidentParticle(aParticle);
66 
67  G4int incidentCode = incidentParticle.getCode();
68  G4double incidentMass = incidentParticle.getMass();
69  G4double incidentTotalEnergy = incidentParticle.getEnergy();
70 
71  // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
72  // DHW 19 May 2011: variable set but not used
73 
74  G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
75 
76  if(incidentKineticEnergy < 1)
77  G4cout << "GHEKaonZeroShortInelastic: incident energy < 1 GeV" << G4endl;
78 
79  if(verboseLevel > 1) {
80  G4cout << "G4HEKaonZeroShortInelastic::ApplyYourself" << G4endl;
81  G4cout << "incident particle " << incidentParticle.getName()
82  << "mass " << incidentMass
83  << "kinetic energy " << incidentKineticEnergy
84  << G4endl;
85  G4cout << "target material with (A,Z) = ("
86  << atomicWeight << "," << atomicNumber << ")" << G4endl;
87  }
88 
89  G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
90  atomicWeight, atomicNumber);
91  if(verboseLevel > 1)
92  G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
93 
94  incidentKineticEnergy -= inelasticity;
95 
96  G4double excitationEnergyGNP = 0.;
97  G4double excitationEnergyDTA = 0.;
98 
99  G4double excitation = NuclearExcitation(incidentKineticEnergy,
100  atomicWeight, atomicNumber,
101  excitationEnergyGNP,
102  excitationEnergyDTA);
103  if(verboseLevel > 1)
104  G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
105  << excitationEnergyDTA << G4endl;
106 
107  incidentKineticEnergy -= excitation;
108  incidentTotalEnergy = incidentKineticEnergy + incidentMass;
109  // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
110  // *(incidentTotalEnergy+incidentMass));
111  // DHW 19 May 2011: variable set but not used
112 
113  G4HEVector targetParticle;
114  if(G4UniformRand() < atomicNumber/atomicWeight) {
115  targetParticle.setDefinition("Proton");
116  } else {
117  targetParticle.setDefinition("Neutron");
118  }
119 
120  G4double targetMass = targetParticle.getMass();
121  G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
122  + targetMass*targetMass
123  + 2.0*targetMass*incidentTotalEnergy);
124  G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
125 
126  G4bool inElastic = true;
127  vecLength = 0;
128 
129  if(verboseLevel > 1)
130  G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
131  << incidentCode << G4endl;
132 
133  G4bool successful = false;
134 
135  // Split K0L into K0 and K0bar
136  if (G4UniformRand() < 0.5)
137  FirstIntInCasAntiKaonZero(inElastic, availableEnergy, pv, vecLength,
138  incidentParticle, targetParticle );
139  else
140  FirstIntInCasKaonZero(inElastic, availableEnergy, pv, vecLength,
141  incidentParticle, targetParticle, atomicWeight );
142 
143  // Do nuclear interaction with either K0 or K0bar
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 
175  if (!successful)
176  ElasticScattering(successful, pv, vecLength,
177  incidentParticle,
178  atomicWeight, atomicNumber);
179 
180  if (!successful)
181  G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
182  << G4endl;
183 
184  // Check for K0, K0bar and change particle types to K0L, K0S if necessary
185  G4int kcode;
186  for (G4int i = 0; i < vecLength; i++) {
187  kcode = pv[i].getCode();
188  if (kcode == KaonZero.getCode() || kcode == AntiKaonZero.getCode()) {
189  if (G4UniformRand() < 0.5)
190  pv[i] = KaonZeroShort;
191  else
192  pv[i] = KaonZeroLong;
193  }
194  }
195 
196  // ................
197 
198  FillParticleChange(pv, vecLength);
199  delete [] pv;
200  theParticleChange.SetStatusChange(stopAndKill);
201  return &theParticleChange;
202 }
203 
204 
205 void
206 G4HEKaonZeroShortInelastic::FirstIntInCasKaonZero(G4bool& inElastic,
207  const G4double availableEnergy,
208  G4HEVector pv[],
209  G4int& vecLen,
210  const G4HEVector& incidentParticle,
211  const G4HEVector& targetParticle,
212  const G4double atomicWeight)
213 
214 // Kaon0 undergoes interaction with nucleon within a nucleus. Check if it is
215 // energetically possible to produce pions/kaons. In not, assume nuclear excitation
216 // occurs and input particle is degraded in energy. No other particles are produced.
217 // If reaction is possible, find the correct number of pions/protons/neutrons
218 // produced using an interpolation to multiplicity data. Replace some pions or
219 // protons/neutrons by kaons or strange baryons according to the average
220 // multiplicity per inelastic reaction.
221 {
222  static const G4double expxu = 82.; // upper bound for arg. of exp
223  static const G4double expxl = -expxu; // lower bound for arg. of exp
224 
225  static const G4double protb = 0.7;
226  static const G4double neutb = 0.7;
227  static const G4double c = 1.25;
228 
229  static const G4int numMul = 1200;
230  static const G4int numSec = 60;
231 
232  G4int neutronCode = Neutron.getCode();
233  G4int protonCode = Proton.getCode();
234 
235  G4int targetCode = targetParticle.getCode();
236  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
237 
238  static G4bool first = true;
239  static G4double protmul[numMul], protnorm[numSec]; // proton constants
240  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
241 
242  // misc. local variables
243  // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
244 
245  G4int i, counter, nt, npos, nneg, nzero;
246 
247  if (first) {
248  // compute normalization constants, this will only be done once
249  first = false;
250  for( i=0; i<numMul; i++ )protmul[i] = 0.0;
251  for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
252  counter = -1;
253  for (npos=0; npos<(numSec/3); npos++) {
254  for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
255  for (nzero=0; nzero<numSec/3; nzero++) {
256  if (++counter < numMul) {
257  nt = npos+nneg+nzero;
258  if( (nt>0) && (nt<=numSec) ) {
259  protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c) ;
260  protnorm[nt-1] += protmul[counter];
261  }
262  }
263  }
264  }
265  }
266 
267  for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
268  for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
269  counter = -1;
270  for (npos=0; npos<numSec/3; npos++) {
271  for (nneg=npos; nneg<=(npos+2); nneg++) {
272  for (nzero=0; nzero<numSec/3; nzero++) {
273  if (++counter < numMul) {
274  nt = npos+nneg+nzero;
275  if( (nt>0) && (nt<=numSec) ) {
276  neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
277  neutnorm[nt-1] += neutmul[counter];
278  }
279  }
280  }
281  }
282  }
283 
284  for (i=0; i<numSec; i++) {
285  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
286  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
287  }
288  } // end of initialization
289 
290 
291  // Initialize the first two particles
292  // the same as beam and target
293  pv[0] = incidentParticle;
294  pv[1] = targetParticle;
295  vecLen = 2;
296 
297  if( !inElastic ) {
298  // quasi-elastic scattering, no pions produced
299  if( targetCode == protonCode) {
300  G4double cech[] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.10,0.09,0.07};
301  G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*5. ) );
302  if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42)) {
303  // charge exchange K+ n -> K0 p
304  pv[0] = KaonPlus;
305  pv[1] = Neutron;
306  }
307  }
308  return;
309  } else if (availableEnergy <= PionPlus.getMass()) {
310  return;
311  }
312 
313  // Inelastic scattering
314 
315  npos = 0, nneg = 0, nzero = 0;
316  G4double eab = availableEnergy;
317  G4int ieab = G4int( eab*5.0 );
318 
319  G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
320  if( (ieab <= 9) && (G4UniformRand() >= supp[ieab])) {
321  // Suppress high multiplicity events at low momentum
322  // only one additional pion will be produced
323  G4double w0, wp, wm, wt, ran;
324  if (targetCode == neutronCode) {
325  // target is a neutron
326  w0 = - sqr(1.+protb)/(2.*c*c);
327  w0 = std::exp(w0);
328  wm = - sqr(-1.+protb)/(2.*c*c);
329  wm = std::exp(wm);
330  w0 = w0/2.;
331  wm = wm*1.5;
332  if (G4UniformRand() < w0/(w0+wm) ) {
333  npos = 0;
334  nneg = 0;
335  nzero = 1;
336  } else {
337  npos = 0;
338  nneg = 1;
339  nzero = 0;
340  }
341 
342  } else {
343  // target is a proton
344  w0 = -sqr(1.+neutb)/(2.*c*c);
345  wp = w0 = std::exp(w0);
346  wm = -sqr(-1.+neutb)/(2.*c*c);
347  wm = std::exp(wm);
348  wt = w0+wp+wm;
349  wp = w0+wp;
350  ran = G4UniformRand();
351  if ( ran < w0/wt) {
352  npos = 0;
353  nneg = 0;
354  nzero = 1;
355  } else if (ran < wp/wt) {
356  npos = 1;
357  nneg = 0;
358  nzero = 0;
359  } else {
360  npos = 0;
361  nneg = 1;
362  nzero = 0;
363  }
364  }
365  } else {
366  // number of total particles vs. centre of mass Energy - 2*proton mass
367 
368  G4double aleab = std::log(availableEnergy);
369  G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
370  + aleab*(0.117712+0.0136912*aleab))) - 2.0;
371 
372  // Normalization constant for kno-distribution.
373  // Calculate first the sum of all constants, check for numerical problems.
374  G4double test, dum, anpn = 0.0;
375 
376  for (nt=1; nt<=numSec; nt++) {
377  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
378  dum = pi*nt/(2.0*n*n);
379  if (std::fabs(dum) < 1.0) {
380  if( test >= 1.0e-10 )anpn += dum*test;
381  } else {
382  anpn += dum*test;
383  }
384  }
385 
386  G4double ran = G4UniformRand();
387  G4double excs = 0.0;
388  if( targetCode == protonCode )
389  {
390  counter = -1;
391  for( npos=0; npos<numSec/3; npos++ )
392  {
393  for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
394  {
395  for (nzero=0; nzero<numSec/3; nzero++) {
396  if (++counter < numMul) {
397  nt = npos+nneg+nzero;
398  if ( (nt>0) && (nt<=numSec) ) {
399  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
400  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
401  if (std::fabs(dum) < 1.0) {
402  if( test >= 1.0e-10 )excs += dum*test;
403  } else {
404  excs += dum*test;
405  }
406  if (ran < excs) goto outOfLoop; //----------------------->
407  }
408  }
409  }
410  }
411  }
412 
413  // 3 previous loops continued to the end
414  inElastic = false; // quasi-elastic scattering
415  return;
416  }
417  else
418  { // target must be a neutron
419  counter = -1;
420  for( npos=0; npos<numSec/3; npos++ )
421  {
422  for( nneg=npos; nneg<=(npos+2); nneg++ )
423  {
424  for (nzero=0; nzero<numSec/3; nzero++) {
425  if (++counter < numMul) {
426  nt = npos+nneg+nzero;
427  if ( (nt>=1) && (nt<=numSec) ) {
428  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
429  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
430  if (std::fabs(dum) < 1.0) {
431  if( test >= 1.0e-10 )excs += dum*test;
432  } else {
433  excs += dum*test;
434  }
435  if (ran < excs) goto outOfLoop; // -------------------------->
436  }
437  }
438  }
439  }
440  }
441  // 3 previous loops continued to the end
442  inElastic = false; // quasi-elastic scattering.
443  return;
444  }
445  }
446  outOfLoop: // <-----------------------------------------------
447 
448  if( targetCode == neutronCode)
449  {
450  if( npos == nneg)
451  {
452  }
453  else if (npos == (nneg-1))
454  {
455  if( G4UniformRand() < 0.5)
456  {
457  pv[0] = KaonPlus;
458  }
459  else
460  {
461  pv[1] = Proton;
462  }
463  }
464  else
465  {
466  pv[0] = KaonPlus;
467  pv[1] = Proton;
468  }
469  }
470  else
471  {
472  if( npos == nneg )
473  {
474  if( G4UniformRand() < 0.25)
475  {
476  pv[0] = KaonPlus;
477  pv[1] = Neutron;
478  }
479  else
480  {
481  }
482  }
483  else if ( npos == (nneg+1))
484  {
485  pv[1] = Neutron;
486  }
487  else
488  {
489  pv[0] = KaonPlus;
490  }
491  }
492 
493  nt = npos + nneg + nzero;
494  while (nt > 0) {
495  G4double ran = G4UniformRand();
496  if (ran < (G4double)npos/nt) {
497  if (npos > 0) {
498  pv[vecLen++] = PionPlus;
499  npos--;
500  }
501  } else if ( ran < (G4double)(npos+nneg)/nt) {
502  if (nneg > 0) {
503  pv[vecLen++] = PionMinus;
504  nneg--;
505  }
506  } else {
507  if (nzero > 0) {
508  pv[vecLen++] = PionZero;
509  nzero--;
510  }
511  }
512  nt = npos + nneg + nzero;
513  }
514 
515  if (verboseLevel > 1) {
516  G4cout << "Particles produced: " ;
517  G4cout << pv[0].getName() << " " ;
518  G4cout << pv[1].getName() << " " ;
519  for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ;
520  G4cout << G4endl;
521  }
522 
523  return;
524 }
525 
526 
527 void
528 G4HEKaonZeroShortInelastic::FirstIntInCasAntiKaonZero(G4bool& inElastic,
529  const G4double availableEnergy,
530  G4HEVector pv[],
531  G4int& vecLen,
532  const G4HEVector& incidentParticle,
533  const G4HEVector& targetParticle)
534 
535 // AntiKaon0 undergoes interaction with nucleon within a nucleus. Check if it is
536 // energetically possible to produce pions/kaons. In not, assume nuclear excitation
537 // occurs and input particle is degraded in energy. No other particles are produced.
538 // If reaction is possible, find the correct number of pions/protons/neutrons
539 // produced using an interpolation to multiplicity data. Replace some pions or
540 // protons/neutrons by kaons or strange baryons according to the average
541 // multiplicity per inelastic reaction.
542 {
543  static const G4double expxu = 82.; // upper bound for arg. of exp
544  static const G4double expxl = -expxu; // lower bound for arg. of exp
545 
546  static const G4double protb = 0.7;
547  static const G4double neutb = 0.7;
548  static const G4double c = 1.25;
549 
550  static const G4int numMul = 1200;
551  static const G4int numSec = 60;
552 
553  G4int neutronCode = Neutron.getCode();
554  G4int protonCode = Proton.getCode();
555  G4int kaonMinusCode = KaonMinus.getCode();
556  G4int kaonZeroCode = KaonZero.getCode();
557  G4int antiKaonZeroCode = AntiKaonZero.getCode();
558 
559  G4int targetCode = targetParticle.getCode();
560  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
561 
562  static G4bool first = true;
563  static G4double protmul[numMul], protnorm[numSec]; // proton constants
564  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
565 
566  // misc. local variables
567  // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
568 
569  G4int i, counter, nt, npos, nneg, nzero;
570 
571  if (first) {
572  // compute normalization constants, this will only be done once
573  first = false;
574  for( i=0; i<numMul; i++ )protmul[i] = 0.0;
575  for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
576  counter = -1;
577  for(npos=0; npos<(numSec/3); npos++) {
578  for(nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
579  for(nzero=0; nzero<numSec/3; nzero++) {
580  if(++counter < numMul) {
581  nt = npos+nneg+nzero;
582  if( (nt>0) && (nt<=numSec) ) {
583  protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c) ;
584  protnorm[nt-1] += protmul[counter];
585  }
586  }
587  }
588  }
589  }
590 
591  for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
592  for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
593  counter = -1;
594  for(npos=0; npos<numSec/3; npos++) {
595  for(nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
596  for(nzero=0; nzero<numSec/3; nzero++) {
597  if(++counter < numMul) {
598  nt = npos+nneg+nzero;
599  if( (nt>0) && (nt<=numSec) ) {
600  neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
601  neutnorm[nt-1] += neutmul[counter];
602  }
603  }
604  }
605  }
606  }
607 
608  for(i=0; i<numSec; i++) {
609  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
610  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
611  }
612  } // end of initialization
613 
614  // initialize the first two particles
615  // the same as beam and target
616  pv[0] = incidentParticle;
617  pv[1] = targetParticle;
618  vecLen = 2;
619 
620  if (!inElastic || (availableEnergy <= PionPlus.getMass()))
621  return;
622 
623  // Inelastic scattering
624 
625  npos = 0, nneg = 0, nzero = 0;
626  G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15};
627  G4int iplab = G4int( incidentTotalMomentum*5.);
628  if( (iplab < 10) && (G4UniformRand() < cech[iplab]) ) {
629  G4int ipl = std::min(19, G4int( incidentTotalMomentum*5.));
630  G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45,
631  0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33};
632  if(G4UniformRand() < cnk0[ipl]) {
633  if(targetCode == protonCode) {
634  return;
635  } else {
636  pv[0] = KaonMinus;
637  pv[1] = Proton;
638  return;
639  }
640  }
641 
642  G4double ran = G4UniformRand();
643  if(targetCode == protonCode) {
644 
645  // target is a proton
646  if( ran < 0.25 ) {
647  ;
648  } else if (ran < 0.50) {
649  pv[0] = PionPlus;
650  pv[1] = SigmaZero;
651  } else if (ran < 0.75) {
652  ;
653  } else {
654  pv[0] = PionPlus;
655  pv[1] = Lambda;
656  }
657  } else {
658 
659  // target is a neutron
660  if( ran < 0.25 ) {
661  pv[0] = PionMinus;
662  pv[1] = SigmaPlus;
663  } else if (ran < 0.50) {
664  pv[0] = PionZero;
665  pv[1] = SigmaZero;
666  } else if (ran < 0.75) {
667  pv[0] = PionPlus;
668  pv[1] = SigmaMinus;
669  } else {
670  pv[0] = PionZero;
671  pv[1] = Lambda;
672  }
673  }
674  return;
675 
676  } else {
677  // number of total particles vs. centre of mass Energy - 2*proton mass
678 
679  G4double aleab = std::log(availableEnergy);
680  G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
681  + aleab*(0.117712+0.0136912*aleab))) - 2.0;
682 
683  // Normalization constant for kno-distribution.
684  // Calculate first the sum of all constants, check for numerical problems.
685  G4double test, dum, anpn = 0.0;
686 
687  for (nt=1; nt<=numSec; nt++) {
688  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
689  dum = pi*nt/(2.0*n*n);
690  if (std::fabs(dum) < 1.0) {
691  if( test >= 1.0e-10 )anpn += dum*test;
692  } else {
693  anpn += dum*test;
694  }
695  }
696 
697  G4double ran = G4UniformRand();
698  G4double excs = 0.0;
699  if (targetCode == protonCode) {
700  counter = -1;
701  for (npos=0; npos<numSec/3; npos++) {
702  for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
703  for (nzero=0; nzero<numSec/3; nzero++) {
704  if (++counter < numMul) {
705  nt = npos+nneg+nzero;
706  if( (nt>0) && (nt<=numSec) ) {
707  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
708  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
709 
710  if (std::fabs(dum) < 1.0) {
711  if( test >= 1.0e-10 )excs += dum*test;
712  } else {
713  excs += dum*test;
714  }
715 
716  if (ran < excs) goto outOfLoop; //----------------------->
717  }
718  }
719  }
720  }
721  }
722  // 3 previous loops continued to the end
723  inElastic = false; // quasi-elastic scattering
724  return;
725 
726  } else { // target must be a neutron
727  counter = -1;
728  for (npos=0; npos<numSec/3; npos++) {
729  for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
730  for (nzero=0; nzero<numSec/3; nzero++) {
731  if (++counter < numMul) {
732  nt = npos+nneg+nzero;
733  if( (nt>=1) && (nt<=numSec) ) {
734  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
735  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
736 
737  if (std::fabs(dum) < 1.0) {
738  if( test >= 1.0e-10 )excs += dum*test;
739  } else {
740  excs += dum*test;
741  }
742 
743  if (ran < excs) goto outOfLoop; // -------------------------->
744  }
745  }
746  }
747  }
748  }
749  // 3 previous loops continued to the end
750  inElastic = false; // quasi-elastic scattering.
751  return;
752  }
753  }
754  outOfLoop: // <------------------------------------------------------------------------
755 
756  if( targetCode == protonCode)
757  {
758  if( npos == nneg)
759  {
760  }
761  else if (npos == (1+nneg))
762  {
763  if( G4UniformRand() < 0.5)
764  {
765  pv[0] = KaonMinus;
766  }
767  else
768  {
769  pv[1] = Neutron;
770  }
771  }
772  else
773  {
774  pv[0] = KaonMinus;
775  pv[1] = Neutron;
776  }
777  }
778  else
779  {
780  if( npos == nneg)
781  {
782  if( G4UniformRand() < 0.75)
783  {
784  }
785  else
786  {
787  pv[0] = KaonMinus;
788  pv[1] = Proton;
789  }
790  }
791  else if ( npos == (1+nneg))
792  {
793  pv[0] = KaonMinus;
794  }
795  else
796  {
797  pv[1] = Proton;
798  }
799  }
800 
801 
802  if( G4UniformRand() < 0.5 )
803  {
804  if( ( (pv[0].getCode() == kaonMinusCode)
805  && (pv[1].getCode() == neutronCode) )
806  || ( (pv[0].getCode() == kaonZeroCode)
807  && (pv[1].getCode() == protonCode) )
808  || ( (pv[0].getCode() == antiKaonZeroCode)
809  && (pv[1].getCode() == protonCode) ) )
810  {
811  G4double ran = G4UniformRand();
812  if( pv[1].getCode() == protonCode)
813  {
814  if(ran < 0.68)
815  {
816  pv[0] = PionPlus;
817  pv[1] = Lambda;
818  }
819  else if (ran < 0.84)
820  {
821  pv[0] = PionZero;
822  pv[1] = SigmaPlus;
823  }
824  else
825  {
826  pv[0] = PionPlus;
827  pv[1] = SigmaZero;
828  }
829  }
830  else
831  {
832  if(ran < 0.68)
833  {
834  pv[0] = PionMinus;
835  pv[1] = Lambda;
836  }
837  else if (ran < 0.84)
838  {
839  pv[0] = PionMinus;
840  pv[1] = SigmaZero;
841  }
842  else
843  {
844  pv[0] = PionZero;
845  pv[1] = SigmaMinus;
846  }
847  }
848  }
849  else
850  {
851  G4double ran = G4UniformRand();
852  if (ran < 0.67)
853  {
854  pv[0] = PionZero;
855  pv[1] = Lambda;
856  }
857  else if (ran < 0.78)
858  {
859  pv[0] = PionMinus;
860  pv[1] = SigmaPlus;
861  }
862  else if (ran < 0.89)
863  {
864  pv[0] = PionZero;
865  pv[1] = SigmaZero;
866  }
867  else
868  {
869  pv[0] = PionPlus;
870  pv[1] = SigmaMinus;
871  }
872  }
873  }
874 
875  nt = npos + nneg + nzero;
876  while ( nt > 0) {
877  G4double ran = G4UniformRand();
878  if ( ran < (G4double)npos/nt) {
879  if( npos > 0 ) {
880  pv[vecLen++] = PionPlus;
881  npos--;
882  }
883  } else if (ran < (G4double)(npos+nneg)/nt) {
884  if( nneg > 0 ) {
885  pv[vecLen++] = PionMinus;
886  nneg--;
887  }
888  } else {
889  if( nzero > 0 ) {
890  pv[vecLen++] = PionZero;
891  nzero--;
892  }
893  }
894  nt = npos + nneg + nzero;
895  }
896 
897  if (verboseLevel > 1) {
898  G4cout << "Particles produced: " ;
899  G4cout << pv[0].getName() << " " ;
900  G4cout << pv[1].getName() << " " ;
901  for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ;
902  G4cout << G4endl;
903  }
904 
905  return;
906 }
Generated on Sat May 25 2013 14:33:56 for Geant4 by   doxygen 1.8.4