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