Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HEKaonZeroLongInelastic.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 // 26 January 2010
37 
38 #include "G4HEKaonZeroLongInelastic.hh"
39 #include "globals.hh"
40 #include "G4ios.hh"
41 #include "G4PhysicalConstants.hh"
42 
43 void G4HEKaonZeroLongInelastic::ModelDescription(std::ostream& outFile) const
44 {
45  outFile << "G4HEKaonZeroLongInelastic is one of the High Energy\n"
46  << "Parameterized (HEP) models used to implement inelastic\n"
47  << "K0L 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 K0L with initial energies\n"
53  << "above 20 GeV.\n";
54 }
55 
56 
57 G4HadFinalState*
58 G4HEKaonZeroLongInelastic::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 << "GHEKaonZeroLongInelastic: incident energy < 1 GeV " << G4endl;
78 
79  if(verboseLevel > 1) {
80  G4cout << "G4HEKaonZeroLongInelastic::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 G4HEKaonZeroLongInelastic::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 G4HEKaonZeroLongInelastic::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 {
544  static const G4double expxu = 82.; // upper bound for arg. of exp
545  static const G4double expxl = -expxu; // lower bound for arg. of exp
546 
547  static const G4double protb = 0.7;
548  static const G4double neutb = 0.7;
549  static const G4double c = 1.25;
550 
551  static const G4int numMul = 1200;
552  static const G4int numSec = 60;
553 
554  G4int neutronCode = Neutron.getCode();
555  G4int protonCode = Proton.getCode();
556  G4int kaonMinusCode = KaonMinus.getCode();
557  G4int kaonZeroCode = KaonZero.getCode();
558  G4int antiKaonZeroCode = AntiKaonZero.getCode();
559 
560  G4int targetCode = targetParticle.getCode();
561  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
562 
563  static G4bool first = true;
564  static G4double protmul[numMul], protnorm[numSec]; // proton constants
565  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
566 
567  // misc. local variables
568  // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
569 
570  G4int i, counter, nt, npos, nneg, nzero;
571 
572  if (first) {
573  // compute normalization constants, this will only be done once
574  first = false;
575  for( i=0; i<numMul; i++ )protmul[i] = 0.0;
576  for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
577  counter = -1;
578  for(npos=0; npos<(numSec/3); npos++) {
579  for(nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
580  for(nzero=0; nzero<numSec/3; nzero++) {
581  if(++counter < numMul) {
582  nt = npos+nneg+nzero;
583  if( (nt>0) && (nt<=numSec) ) {
584  protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c) ;
585  protnorm[nt-1] += protmul[counter];
586  }
587  }
588  }
589  }
590  }
591 
592  for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
593  for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
594  counter = -1;
595  for(npos=0; npos<numSec/3; npos++) {
596  for(nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
597  for(nzero=0; nzero<numSec/3; nzero++) {
598  if(++counter < numMul) {
599  nt = npos+nneg+nzero;
600  if( (nt>0) && (nt<=numSec) ) {
601  neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
602  neutnorm[nt-1] += neutmul[counter];
603  }
604  }
605  }
606  }
607  }
608 
609  for(i=0; i<numSec; i++) {
610  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
611  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
612  }
613  } // end of initialization
614 
615  // initialize the first two particles
616  // the same as beam and target
617  pv[0] = incidentParticle;
618  pv[1] = targetParticle;
619  vecLen = 2;
620 
621  if (!inElastic || (availableEnergy <= PionPlus.getMass()))
622  return;
623 
624  // Inelastic scattering
625 
626  npos = 0, nneg = 0, nzero = 0;
627  G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15};
628  G4int iplab = G4int( incidentTotalMomentum*5.);
629  if ((iplab < 10) && (G4UniformRand() < cech[iplab]) ) {
630  G4int ipl = std::min(19, G4int( incidentTotalMomentum*5.));
631  G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45,
632  0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33};
633  if (G4UniformRand() < cnk0[ipl]) {
634  if(targetCode == protonCode) {
635  return;
636  } else {
637  pv[0] = KaonMinus;
638  pv[1] = Proton;
639  return;
640  }
641  }
642 
643  G4double ran = G4UniformRand();
644  if(targetCode == protonCode) {
645 
646  // target is a proton
647  if( ran < 0.25 ) {
648  ;
649  } else if (ran < 0.50) {
650  pv[0] = PionPlus;
651  pv[1] = SigmaZero;
652  } else if (ran < 0.75) {
653  ;
654  } else {
655  pv[0] = PionPlus;
656  pv[1] = Lambda;
657  }
658  } else {
659 
660  // target is a neutron
661  if( ran < 0.25 ) {
662  pv[0] = PionMinus;
663  pv[1] = SigmaPlus;
664  } else if (ran < 0.50) {
665  pv[0] = PionZero;
666  pv[1] = SigmaZero;
667  } else if (ran < 0.75) {
668  pv[0] = PionPlus;
669  pv[1] = SigmaMinus;
670  } else {
671  pv[0] = PionZero;
672  pv[1] = Lambda;
673  }
674  }
675  return;
676 
677  } else {
678  // number of total particles vs. centre of mass Energy - 2*proton mass
679 
680  G4double aleab = std::log(availableEnergy);
681  G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
682  + aleab*(0.117712+0.0136912*aleab))) - 2.0;
683 
684  // Normalization constant for kno-distribution.
685  // Calculate first the sum of all constants, check for numerical problems.
686  G4double test, dum, anpn = 0.0;
687 
688  for (nt=1; nt<=numSec; nt++) {
689  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
690  dum = pi*nt/(2.0*n*n);
691  if (std::fabs(dum) < 1.0) {
692  if( test >= 1.0e-10 )anpn += dum*test;
693  } else {
694  anpn += dum*test;
695  }
696  }
697 
698  G4double ran = G4UniformRand();
699  G4double excs = 0.0;
700  if (targetCode == protonCode) {
701  counter = -1;
702  for (npos=0; npos<numSec/3; npos++) {
703  for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
704  for (nzero=0; nzero<numSec/3; nzero++) {
705  if (++counter < numMul) {
706  nt = npos+nneg+nzero;
707  if( (nt>0) && (nt<=numSec) ) {
708  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
709  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
710 
711  if (std::fabs(dum) < 1.0) {
712  if( test >= 1.0e-10 )excs += dum*test;
713  } else {
714  excs += dum*test;
715  }
716 
717  if (ran < excs) goto outOfLoop; //----------------------->
718  }
719  }
720  }
721  }
722  }
723  // 3 previous loops continued to the end
724  inElastic = false; // quasi-elastic scattering
725  return;
726 
727  } else { // target must be a neutron
728  counter = -1;
729  for (npos=0; npos<numSec/3; npos++) {
730  for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
731  for (nzero=0; nzero<numSec/3; nzero++) {
732  if (++counter < numMul) {
733  nt = npos+nneg+nzero;
734  if( (nt>=1) && (nt<=numSec) ) {
735  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
736  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
737 
738  if (std::fabs(dum) < 1.0) {
739  if( test >= 1.0e-10 )excs += dum*test;
740  } else {
741  excs += dum*test;
742  }
743 
744  if (ran < excs) goto outOfLoop; // -------------------------->
745  }
746  }
747  }
748  }
749  }
750  // 3 previous loops continued to the end
751  inElastic = false; // quasi-elastic scattering.
752  return;
753  }
754  }
755  outOfLoop: // <------------------------------------------------------------------------
756 
757  if( targetCode == protonCode)
758  {
759  if( npos == nneg)
760  {
761  }
762  else if (npos == (1+nneg))
763  {
764  if( G4UniformRand() < 0.5)
765  {
766  pv[0] = KaonMinus;
767  }
768  else
769  {
770  pv[1] = Neutron;
771  }
772  }
773  else
774  {
775  pv[0] = KaonMinus;
776  pv[1] = Neutron;
777  }
778  }
779  else
780  {
781  if( npos == nneg)
782  {
783  if( G4UniformRand() < 0.75)
784  {
785  }
786  else
787  {
788  pv[0] = KaonMinus;
789  pv[1] = Proton;
790  }
791  }
792  else if ( npos == (1+nneg))
793  {
794  pv[0] = KaonMinus;
795  }
796  else
797  {
798  pv[1] = Proton;
799  }
800  }
801 
802 
803  if( G4UniformRand() < 0.5 )
804  {
805  if( ( (pv[0].getCode() == kaonMinusCode)
806  && (pv[1].getCode() == neutronCode) )
807  || ( (pv[0].getCode() == kaonZeroCode)
808  && (pv[1].getCode() == protonCode) )
809  || ( (pv[0].getCode() == antiKaonZeroCode)
810  && (pv[1].getCode() == protonCode) ) )
811  {
812  G4double ran = G4UniformRand();
813  if( pv[1].getCode() == protonCode)
814  {
815  if(ran < 0.68)
816  {
817  pv[0] = PionPlus;
818  pv[1] = Lambda;
819  }
820  else if (ran < 0.84)
821  {
822  pv[0] = PionZero;
823  pv[1] = SigmaPlus;
824  }
825  else
826  {
827  pv[0] = PionPlus;
828  pv[1] = SigmaZero;
829  }
830  }
831  else
832  {
833  if(ran < 0.68)
834  {
835  pv[0] = PionMinus;
836  pv[1] = Lambda;
837  }
838  else if (ran < 0.84)
839  {
840  pv[0] = PionMinus;
841  pv[1] = SigmaZero;
842  }
843  else
844  {
845  pv[0] = PionZero;
846  pv[1] = SigmaMinus;
847  }
848  }
849  }
850  else
851  {
852  G4double ran = G4UniformRand();
853  if (ran < 0.67)
854  {
855  pv[0] = PionZero;
856  pv[1] = Lambda;
857  }
858  else if (ran < 0.78)
859  {
860  pv[0] = PionMinus;
861  pv[1] = SigmaPlus;
862  }
863  else if (ran < 0.89)
864  {
865  pv[0] = PionZero;
866  pv[1] = SigmaZero;
867  }
868  else
869  {
870  pv[0] = PionPlus;
871  pv[1] = SigmaMinus;
872  }
873  }
874  }
875 
876  nt = npos + nneg + nzero;
877  while ( nt > 0) {
878  G4double ran = G4UniformRand();
879  if ( ran < (G4double)npos/nt) {
880  if( npos > 0 ) {
881  pv[vecLen++] = PionPlus;
882  npos--;
883  }
884  } else if (ran < (G4double)(npos+nneg)/nt) {
885  if( nneg > 0 ) {
886  pv[vecLen++] = PionMinus;
887  nneg--;
888  }
889  } else {
890  if( nzero > 0 ) {
891  pv[vecLen++] = PionZero;
892  nzero--;
893  }
894  }
895  nt = npos + nneg + nzero;
896  }
897 
898  if (verboseLevel > 1) {
899  G4cout << "Particles produced: " ;
900  G4cout << pv[0].getName() << " " ;
901  G4cout << pv[1].getName() << " " ;
902  for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ;
903  G4cout << G4endl;
904  }
905 
906  return;
907 }
Generated on Sat May 25 2013 14:33:56 for Geant4 by   doxygen 1.8.4