Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HEAntiLambdaInelastic.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 "G4HEAntiLambdaInelastic.hh"
38 #include "globals.hh"
39 #include "G4ios.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 
43 G4HEAntiLambdaInelastic::G4HEAntiLambdaInelastic(const G4String& name)
44  : G4HEInelastic(name)
45 {
46  vecLength = 0;
47  theMinEnergy = 20*GeV;
48  theMaxEnergy = 10*TeV;
49  MAXPART = 2048;
50  verboseLevel = 0;
51  G4cout << "WARNING: model G4HEAntiLambdaInelastic is being deprecated and will\n"
52  << "disappear in Geant4 version 10.0" << G4endl;
53 }
54 
55 
56 void G4HEAntiLambdaInelastic::ModelDescription(std::ostream& outFile) const
57 {
58  outFile << "G4HEAntiLambdaInelastic is one of the High Energy\n"
59  << "Parameterized (HEP) models used to implement inelastic\n"
60  << "anti-Lambda scattering from nuclei. It is a re-engineered\n"
61  << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
62  << "initial collision products into backward- and forward-going\n"
63  << "clusters which are then decayed into final state hadrons.\n"
64  << "The model does not conserve energy on an event-by-event\n"
65  << "basis. It may be applied to anti-Lambdas with initial energies\n"
66  << "above 20 GeV.\n";
67 }
68 
69 
70 G4HadFinalState*
71 G4HEAntiLambdaInelastic::ApplyYourself(const G4HadProjectile &aTrack,
72  G4Nucleus &targetNucleus)
73 {
74  G4HEVector* pv = new G4HEVector[MAXPART];
75  const G4HadProjectile *aParticle = &aTrack;
76  const G4double atomicWeight = targetNucleus.GetA_asInt();
77  const G4double atomicNumber = targetNucleus.GetZ_asInt();
78  G4HEVector incidentParticle(aParticle);
79 
80  G4int incidentCode = incidentParticle.getCode();
81  G4double incidentMass = incidentParticle.getMass();
82  G4double incidentTotalEnergy = incidentParticle.getEnergy();
83 
84  G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
85 
86  if (incidentKineticEnergy < 1.)
87  G4cout << "GHEAntiLambdaInelastic: incident energy < 1 GeV" << G4endl;
88 
89  if (verboseLevel > 1) {
90  G4cout << "G4HEAntiLambdaInelastic::ApplyYourself" << G4endl;
91  G4cout << "incident particle " << incidentParticle.getName()
92  << "mass " << incidentMass
93  << "kinetic energy " << incidentKineticEnergy
94  << G4endl;
95  G4cout << "target material with (A,Z) = ("
96  << atomicWeight << "," << atomicNumber << ")" << G4endl;
97  }
98 
99  G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
100  atomicWeight, atomicNumber);
101  if (verboseLevel > 1)
102  G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
103 
104  incidentKineticEnergy -= inelasticity;
105 
106  G4double excitationEnergyGNP = 0.;
107  G4double excitationEnergyDTA = 0.;
108 
109  G4double excitation = NuclearExcitation(incidentKineticEnergy,
110  atomicWeight, atomicNumber,
111  excitationEnergyGNP,
112  excitationEnergyDTA);
113  if (verboseLevel > 1)
114  G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
115  << excitationEnergyDTA << G4endl;
116 
117  incidentKineticEnergy -= excitation;
118  incidentTotalEnergy = incidentKineticEnergy + incidentMass;
119  // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
120  // *(incidentTotalEnergy+incidentMass));
121  // DHW 19 May 2011: variable set but not used
122 
123  G4HEVector targetParticle;
124  if (G4UniformRand() < atomicNumber/atomicWeight) {
125  targetParticle.setDefinition("Proton");
126  } else {
127  targetParticle.setDefinition("Neutron");
128  }
129 
130  G4double targetMass = targetParticle.getMass();
131  G4double centerOfMassEnergy =
132  std::sqrt( incidentMass*incidentMass + targetMass*targetMass
133  + 2.0*targetMass*incidentTotalEnergy);
134  G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
135 
136  G4bool inElastic = true;
137  vecLength = 0;
138 
139  if (verboseLevel > 1)
140  G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
141  << incidentCode << G4endl;
142 
143  G4bool successful = false;
144 
145  FirstIntInCasAntiLambda(inElastic, availableEnergy, pv, vecLength,
146  incidentParticle, targetParticle, atomicWeight);
147 
148  if (verboseLevel > 1)
149  G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
150 
151  if ((vecLength > 0) && (availableEnergy > 1.))
152  StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
153  pv, vecLength,
154  incidentParticle, targetParticle);
155  HighEnergyCascading(successful, pv, vecLength,
156  excitationEnergyGNP, excitationEnergyDTA,
157  incidentParticle, targetParticle,
158  atomicWeight, atomicNumber);
159  if (!successful)
160  HighEnergyClusterProduction(successful, pv, vecLength,
161  excitationEnergyGNP, excitationEnergyDTA,
162  incidentParticle, targetParticle,
163  atomicWeight, atomicNumber);
164  if (!successful)
165  MediumEnergyCascading(successful, pv, vecLength,
166  excitationEnergyGNP, excitationEnergyDTA,
167  incidentParticle, targetParticle,
168  atomicWeight, atomicNumber);
169 
170  if (!successful)
171  MediumEnergyClusterProduction(successful, pv, vecLength,
172  excitationEnergyGNP, excitationEnergyDTA,
173  incidentParticle, targetParticle,
174  atomicWeight, atomicNumber);
175  if (!successful)
176  QuasiElasticScattering(successful, pv, vecLength,
177  excitationEnergyGNP, excitationEnergyDTA,
178  incidentParticle, targetParticle,
179  atomicWeight, atomicNumber);
180  if (!successful)
181  ElasticScattering(successful, pv, vecLength,
182  incidentParticle,
183  atomicWeight, atomicNumber);
184 
185  if (!successful)
186  G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
187  << G4endl;
188 
189  FillParticleChange(pv, vecLength);
190  delete [] pv;
191  theParticleChange.SetStatusChange(stopAndKill);
192  return & theParticleChange;
193 }
194 
195 
196 void
197 G4HEAntiLambdaInelastic::FirstIntInCasAntiLambda(G4bool& inElastic,
198  const G4double availableEnergy,
199  G4HEVector pv[],
200  G4int &vecLen,
201  const G4HEVector& incidentParticle,
202  const G4HEVector& targetParticle,
203  const G4double atomicWeight)
204 
205 // AntiLambda undergoes interaction with nucleon within a nucleus.
206 // Check if it is energetically possible to produce pions/kaons. If not,
207 // assume nuclear excitation occurs and input particle is degraded in
208 // energy. No other particles are produced.
209 // If reaction is possible, find the correct number of pions/protons/neutrons
210 // produced using an interpolation to multiplicity data. Replace some pions or
211 // protons/neutrons by kaons or strange baryons according to the average
212 // multiplicity per inelastic reaction.
213 {
214  static const G4double expxu = 82.; // upper bound for arg. of exp
215  static const G4double expxl = -expxu; // lower bound for arg. of exp
216 
217  static const G4double protb = 0.7;
218  static const G4double neutb = 0.7;
219  static const G4double c = 1.25;
220 
221  static const G4int numMul = 1200;
222  static const G4int numMulAn = 400;
223  static const G4int numSec = 60;
224 
225  G4int protonCode = Proton.getCode();
226 
227  G4int targetCode = targetParticle.getCode();
228  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
229 
230  static G4bool first = true;
231  static G4double protmul[numMul], protnorm[numSec]; // proton constants
232  static G4double protmulAn[numMulAn],protnormAn[numSec];
233  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
234  static G4double neutmulAn[numMulAn],neutnormAn[numSec];
235 
236  // misc. local variables
237  // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
238 
239  G4int i, counter, nt, npos, nneg, nzero;
240 
241  if( first )
242  { // compute normalization constants, this will only be done once
243  first = false;
244  for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
245  for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
246  counter = -1;
247  for( npos=0; npos<(numSec/3); npos++ )
248  {
249  for( nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++ )
250  {
251  for( nzero=0; nzero<numSec/3; nzero++ )
252  {
253  if( ++counter < numMul )
254  {
255  nt = npos+nneg+nzero;
256  if( (nt>0) && (nt<=numSec) )
257  {
258  protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
259  protnorm[nt-1] += protmul[counter];
260  }
261  }
262  }
263  }
264  }
265  for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
266  for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
267  counter = -1;
268  for( npos=0; npos<numSec/3; npos++ )
269  {
270  for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
271  {
272  for( nzero=0; nzero<numSec/3; nzero++ )
273  {
274  if( ++counter < numMul )
275  {
276  nt = npos+nneg+nzero;
277  if( (nt>0) && (nt<=numSec) )
278  {
279  neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
280  neutnorm[nt-1] += neutmul[counter];
281  }
282  }
283  }
284  }
285  }
286  for( i=0; i<numSec; i++ )
287  {
288  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
289  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
290  }
291  // annihilation
292  for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0;
293  for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0;
294  counter = -1;
295  for( npos=1; npos<(numSec/3); npos++ )
296  {
297  nneg = std::max(0,npos-1);
298  for( nzero=0; nzero<numSec/3; nzero++ )
299  {
300  if( ++counter < numMulAn )
301  {
302  nt = npos+nneg+nzero;
303  if( (nt>1) && (nt<=numSec) )
304  {
305  protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
306  protnormAn[nt-1] += protmulAn[counter];
307  }
308  }
309  }
310  }
311  for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0;
312  for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0;
313  counter = -1;
314  for( npos=0; npos<numSec/3; npos++ )
315  {
316  nneg = npos;
317  for( nzero=0; nzero<numSec/3; nzero++ )
318  {
319  if( ++counter < numMulAn )
320  {
321  nt = npos+nneg+nzero;
322  if( (nt>1) && (nt<=numSec) )
323  {
324  neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
325  neutnormAn[nt-1] += neutmulAn[counter];
326  }
327  }
328  }
329  }
330  for( i=0; i<numSec; i++ )
331  {
332  if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
333  if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
334  }
335  } // end of initialization
336 
337 
338  // initialize the first two places
339  // the same as beam and target
340  pv[0] = incidentParticle;
341  pv[1] = targetParticle;
342  vecLen = 2;
343 
344  if( !inElastic )
345  { // some two-body reactions
346  G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
347 
348  G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5));
349  if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
350  {
351  G4double ran = G4UniformRand();
352 
353  if ( targetCode == protonCode)
354  {
355  if(ran < 0.2)
356  {
357  pv[0] = AntiSigmaZero;
358  }
359  else if (ran < 0.4)
360  {
361  pv[0] = AntiSigmaMinus;
362  pv[1] = Neutron;
363  }
364  else if (ran < 0.6)
365  {
366  pv[0] = Proton;
367  pv[1] = AntiLambda;
368  }
369  else if (ran < 0.8)
370  {
371  pv[0] = Proton;
372  pv[1] = AntiSigmaZero;
373  }
374  else
375  {
376  pv[0] = Neutron;
377  pv[1] = AntiSigmaMinus;
378  }
379  }
380  else
381  {
382  if (ran < 0.2)
383  {
384  pv[0] = AntiSigmaZero;
385  }
386  else if (ran < 0.4)
387  {
388  pv[0] = AntiSigmaPlus;
389  pv[1] = Proton;
390  }
391  else if (ran < 0.6)
392  {
393  pv[0] = Neutron;
394  pv[1] = AntiLambda;
395  }
396  else if (ran < 0.8)
397  {
398  pv[0] = Neutron;
399  pv[1] = AntiSigmaZero;
400  }
401  else
402  {
403  pv[0] = Proton;
404  pv[1] = AntiSigmaPlus;
405  }
406  }
407  }
408  return;
409  }
410  else if (availableEnergy <= PionPlus.getMass())
411  return;
412 
413  // inelastic scattering
414 
415  npos = 0; nneg = 0; nzero = 0;
416  G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88,
417  0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40,
418  0.39, 0.36, 0.33, 0.10, 0.01};
419  G4int iplab = G4int( incidentTotalMomentum*10.);
420  if ( iplab > 9) iplab = 10 + G4int( (incidentTotalMomentum -1.)*5. );
421  if ( iplab > 14) iplab = 15 + G4int( incidentTotalMomentum -2. );
422  if ( iplab > 22) iplab = 23 + G4int( (incidentTotalMomentum -10.)/10.);
423  iplab = std::min(24, iplab);
424 
425  if (G4UniformRand() > anhl[iplab]) { // non- annihilation channels
426 
427  // number of total particles vs. centre of mass Energy - 2*proton mass
428  G4double aleab = std::log(availableEnergy);
429  G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
430  + aleab*(0.117712+0.0136912*aleab))) - 2.0;
431 
432  // normalization constant for kno-distribution.
433  // calculate first the sum of all constants, check for numerical problems.
434  G4double test, dum, anpn = 0.0;
435 
436  for (nt = 1; nt <= numSec; nt++) {
437  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
438  dum = pi*nt/(2.0*n*n);
439  if (std::fabs(dum) < 1.0) {
440  if (test >= 1.0e-10) anpn += dum*test;
441  } else {
442  anpn += dum*test;
443  }
444  }
445 
446  G4double ran = G4UniformRand();
447  G4double excs = 0.0;
448  if (targetCode == protonCode) {
449  counter = -1;
450  for (npos = 0; npos < numSec/3; npos++) {
451  for (nneg = std::max(0,npos-2); nneg <= (npos+1); nneg++) {
452  for (nzero = 0; nzero < numSec/3; nzero++) {
453  if (++counter < numMul) {
454  nt = npos+nneg+nzero;
455  if ((nt > 0) && (nt <= numSec) ) {
456  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
457  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
458 
459  if (std::fabs(dum) < 1.0) {
460  if (test >= 1.0e-10) excs += dum*test;
461  } else {
462  excs += dum*test;
463  }
464 
465  if (ran < excs) goto outOfLoop; //----------------------->
466  }
467  }
468  }
469  }
470  }
471  // 3 previous loops continued to the end
472  inElastic = false; // quasi-elastic scattering
473  return;
474  } else { // target must be a neutron
475  counter = -1;
476  for( npos=0; npos<numSec/3; npos++ )
477  {
478  for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
479  {
480  for( nzero=0; nzero<numSec/3; nzero++ )
481  {
482  if( ++counter < numMul )
483  {
484  nt = npos+nneg+nzero;
485  if( (nt>0) && (nt<=numSec) )
486  {
487  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
488  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
489  if (std::fabs(dum) < 1.0) {
490  if( test >= 1.0e-10 )excs += dum*test;
491  } else {
492  excs += dum*test;
493  }
494 
495  if (ran < excs) goto outOfLoop; // -------------------------->
496  }
497  }
498  }
499  }
500  }
501  // 3 previous loops continued to the end
502  inElastic = false; // quasi-elastic scattering.
503  return;
504  }
505 
506  outOfLoop: // <------------------------------------------------------------------------
507 
508  ran = G4UniformRand();
509 
510  if( targetCode == protonCode)
511  {
512  if( npos == nneg)
513  {
514  if (ran < 0.40)
515  {
516  }
517  else if (ran < 0.8)
518  {
519  pv[0] = AntiSigmaZero;
520  }
521  else
522  {
523  pv[0] = AntiSigmaMinus;
524  pv[1] = Neutron;
525  }
526  }
527  else if (npos == (nneg+1))
528  {
529  if( ran < 0.25)
530  {
531  pv[1] = Neutron;
532  }
533  else if (ran < 0.5)
534  {
535  pv[0] = AntiSigmaZero;
536  pv[1] = Neutron;
537  }
538  else
539  {
540  pv[0] = AntiSigmaPlus;
541  }
542  }
543  else if (npos == (nneg-1))
544  {
545  pv[0] = AntiSigmaMinus;
546  }
547  else
548  {
549  pv[0] = AntiSigmaPlus;
550  pv[1] = Neutron;
551  }
552  }
553  else
554  {
555  if( npos == nneg)
556  {
557  if (ran < 0.4)
558  {
559  }
560  else if(ran < 0.8)
561  {
562  pv[0] = AntiSigmaZero;
563  }
564  else
565  {
566  pv[0] = AntiSigmaPlus;
567  pv[1] = Proton;
568  }
569  }
570  else if ( npos == (nneg-1))
571  {
572  if (ran < 0.5)
573  {
574  pv[0] = AntiSigmaMinus;
575  }
576  else if (ran < 0.75)
577  {
578  pv[1] = Proton;
579  }
580  else
581  {
582  pv[0] = AntiSigmaZero;
583  pv[1] = Proton;
584  }
585  }
586  else if (npos == (nneg+1))
587  {
588  pv[0] = AntiSigmaPlus;
589  }
590  else
591  {
592  pv[0] = AntiSigmaMinus;
593  pv[1] = Proton;
594  }
595  }
596 
597  }
598  else // annihilation
599  {
600  if ( availableEnergy > 2. * PionPlus.getMass() )
601  {
602 
603  G4double aleab = std::log(availableEnergy);
604  G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
605  + aleab*(0.117712+0.0136912*aleab))) - 2.0;
606 
607  // normalization constant for kno-distribution.
608  // calculate first the sum of all constants, check for numerical problems.
609  G4double test, dum, anpn = 0.0;
610 
611  for (nt=2; nt<=numSec; nt++) {
612  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
613  dum = pi*nt/(2.0*n*n);
614 
615  if (std::fabs(dum) < 1.0) {
616  if( test >= 1.0e-10 )anpn += dum*test;
617  } else {
618  anpn += dum*test;
619  }
620  }
621 
622  G4double ran = G4UniformRand();
623  G4double excs = 0.0;
624  if (targetCode == protonCode) {
625  counter = -1;
626  for (npos=1; npos<numSec/3; npos++) {
627  nneg = npos-1;
628  for( nzero=0; nzero<numSec/3; nzero++ )
629  {
630  if( ++counter < numMulAn )
631  {
632  nt = npos+nneg+nzero;
633  if( (nt>1) && (nt<=numSec) ) {
634  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
635  dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
636 
637  if (std::fabs(dum) < 1.0) {
638  if( test >= 1.0e-10 )excs += dum*test;
639  } else {
640  excs += dum*test;
641  }
642 
643  if (ran < excs) goto outOfLoopAn; //----------------------->
644  }
645  }
646  }
647  }
648  // 3 previous loops continued to the end
649  inElastic = false; // quasi-elastic scattering
650  return;
651 
652  } else { // target must be a neutron
653  counter = -1;
654  for (npos=0; npos<numSec/3; npos++) {
655  nneg = npos;
656  for( nzero=0; nzero<numSec/3; nzero++ ) {
657  if (++counter < numMulAn) {
658  nt = npos+nneg+nzero;
659  if( (nt>1) && (nt<=numSec) ) {
660  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
661  dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
662 
663  if (std::fabs(dum) < 1.0) {
664  if( test >= 1.0e-10 )excs += dum*test;
665  } else {
666  excs += dum*test;
667  }
668 
669  if (ran < excs) goto outOfLoopAn; // -------------------------->
670  }
671  }
672  }
673  }
674 
675  inElastic = false; // quasi-elastic scattering.
676  return;
677  }
678  outOfLoopAn: // <---------------------------------------------------------
679  vecLen = 0;
680  }
681  }
682 
683  nt = npos + nneg + nzero;
684  while ( nt > 0)
685  {
686  G4double ran = G4UniformRand();
687  if ( ran < (G4double)npos/nt)
688  {
689  if( npos > 0 )
690  { pv[vecLen++] = PionPlus;
691  npos--;
692  }
693  }
694  else if ( ran < (G4double)(npos+nneg)/nt)
695  {
696  if( nneg > 0 )
697  {
698  pv[vecLen++] = PionMinus;
699  nneg--;
700  }
701  }
702  else
703  {
704  if( nzero > 0 )
705  {
706  pv[vecLen++] = PionZero;
707  nzero--;
708  }
709  }
710  nt = npos + nneg + nzero;
711  }
712  if (verboseLevel > 1)
713  {
714  G4cout << "Particles produced: " ;
715  G4cout << pv[0].getName() << " " ;
716  G4cout << pv[1].getName() << " " ;
717  for (i=2; i < vecLen; i++)
718  {
719  G4cout << pv[i].getName() << " " ;
720  }
721  G4cout << G4endl;
722  }
723  return;
724  }
725 
726 
727 
728 
729 
730 
731 
732 
733 
734 
Generated on Sat May 25 2013 14:33:56 for Geant4 by   doxygen 1.8.4