Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HEAntiXiZeroInelastic.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 "G4HEAntiXiZeroInelastic.hh"
38 #include "globals.hh"
39 #include "G4ios.hh"
40 #include "G4PhysicalConstants.hh"
41 
42 void G4HEAntiXiZeroInelastic::ModelDescription(std::ostream& outFile) const
43 {
44  outFile << "G4HEAntiXiZeroInelastic is one of the High Energy\n"
45  << "Parameterized (HEP) models used to implement inelastic\n"
46  << "anti-Xi0 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-Xi0 with initial energies\n"
52  << "above 20 GeV.\n";
53 }
54 
55 
56 G4HadFinalState*
57 G4HEAntiXiZeroInelastic::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 << "GHEAntiXiZeroInelastic: incident energy < 1 GeV" << G4endl;
80 
81  if (verboseLevel > 1) {
82  G4cout << "G4HEAntiXiZeroInelastic::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  FirstIntInCasAntiXiZero(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 
148  HighEnergyCascading(successful, pv, vecLength,
149  excitationEnergyGNP, excitationEnergyDTA,
150  incidentParticle, targetParticle,
151  atomicWeight, atomicNumber);
152  if (!successful)
153  HighEnergyClusterProduction(successful, pv, vecLength,
154  excitationEnergyGNP, excitationEnergyDTA,
155  incidentParticle, targetParticle,
156  atomicWeight, atomicNumber);
157  if (!successful)
158  MediumEnergyCascading(successful, pv, vecLength,
159  excitationEnergyGNP, excitationEnergyDTA,
160  incidentParticle, targetParticle,
161  atomicWeight, atomicNumber);
162 
163  if (!successful)
164  MediumEnergyClusterProduction(successful, pv, vecLength,
165  excitationEnergyGNP, excitationEnergyDTA,
166  incidentParticle, targetParticle,
167  atomicWeight, atomicNumber);
168  if (!successful)
169  QuasiElasticScattering(successful, pv, vecLength,
170  excitationEnergyGNP, excitationEnergyDTA,
171  incidentParticle, targetParticle,
172  atomicWeight, atomicNumber);
173  if (!successful)
174  ElasticScattering(successful, pv, vecLength,
175  incidentParticle,
176  atomicWeight, atomicNumber);
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 G4HEAntiXiZeroInelastic::FirstIntInCasAntiXiZero(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 // AntiXi0 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 // ( decay Xi0 -> L Pi > 99 % )
201 {
202  static const G4double expxu = 82.; // upper bound for arg. of exp
203  static const G4double expxl = -expxu; // lower bound for arg. of exp
204 
205  static const G4double protb = 0.7;
206  static const G4double neutb = 0.7;
207  static const G4double c = 1.25;
208 
209  static const G4int numMul = 1200;
210  static const G4int numMulAn = 400;
211  static const G4int numSec = 60;
212 
213  G4int protonCode = Proton.getCode();
214 
215  G4int targetCode = targetParticle.getCode();
216  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
217 
218  static G4bool first = true;
219  static G4double protmul[numMul], protnorm[numSec]; // proton constants
220  static G4double protmulAn[numMulAn],protnormAn[numSec];
221  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
222  static G4double neutmulAn[numMulAn],neutnormAn[numSec];
223 
224  // misc. local variables
225  // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
226 
227  G4int i, counter, nt, npos, nneg, nzero;
228 
229  if( first )
230  { // compute normalization constants, this will only be done once
231  first = false;
232  for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
233  for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
234  counter = -1;
235  for( npos=0; npos<(numSec/3); npos++ )
236  {
237  for( nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++ )
238  {
239  for( nzero=0; nzero<numSec/3; nzero++ )
240  {
241  if( ++counter < numMul )
242  {
243  nt = npos+nneg+nzero;
244  if( (nt>0) && (nt<=numSec) )
245  {
246  protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
247  protnorm[nt-1] += protmul[counter];
248  }
249  }
250  }
251  }
252  }
253  for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
254  for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
255  counter = -1;
256  for( npos=0; npos<numSec/3; npos++ )
257  {
258  for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
259  {
260  for( nzero=0; nzero<numSec/3; nzero++ )
261  {
262  if( ++counter < numMul )
263  {
264  nt = npos+nneg+nzero;
265  if( (nt>0) && (nt<=numSec) )
266  {
267  neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
268  neutnorm[nt-1] += neutmul[counter];
269  }
270  }
271  }
272  }
273  }
274  for( i=0; i<numSec; i++ )
275  {
276  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
277  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
278  }
279  // annihilation
280  for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0;
281  for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0;
282  counter = -1;
283  for( npos=1; npos<(numSec/3); npos++ )
284  {
285  nneg = std::max(0,npos-1);
286  for( nzero=0; nzero<numSec/3; nzero++ )
287  {
288  if( ++counter < numMulAn )
289  {
290  nt = npos+nneg+nzero;
291  if( (nt>1) && (nt<=numSec) )
292  {
293  protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
294  protnormAn[nt-1] += protmulAn[counter];
295  }
296  }
297  }
298  }
299  for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0;
300  for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0;
301  counter = -1;
302  for( npos=0; npos<numSec/3; npos++ )
303  {
304  nneg = npos;
305  for( nzero=0; nzero<numSec/3; nzero++ )
306  {
307  if( ++counter < numMulAn )
308  {
309  nt = npos+nneg+nzero;
310  if( (nt>1) && (nt<=numSec) )
311  {
312  neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
313  neutnormAn[nt-1] += neutmulAn[counter];
314  }
315  }
316  }
317  }
318  for( i=0; i<numSec; i++ )
319  {
320  if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
321  if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
322  }
323  } // end of initialization
324 
325 
326  // initialize the first two places
327  // the same as beam and target
328  pv[0] = incidentParticle;
329  pv[1] = targetParticle;
330  vecLen = 2;
331 
332  if( !inElastic )
333  { // some two-body reactions
334  G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
335 
336  G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5 ));
337  if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
338  {
339  G4double ran = G4UniformRand();
340 
341  if ( targetCode == protonCode)
342  {
343  if(ran < 0.2)
344  {
345  pv[0] = AntiSigmaZero;
346  }
347  else if (ran < 0.4)
348  {
349  pv[0] = AntiSigmaMinus;
350  pv[1] = Neutron;
351  }
352  else if (ran < 0.6)
353  {
354  pv[0] = Proton;
355  pv[1] = AntiLambda;
356  }
357  else if (ran < 0.8)
358  {
359  pv[0] = Proton;
360  pv[1] = AntiSigmaZero;
361  }
362  else
363  {
364  pv[0] = Neutron;
365  pv[1] = AntiSigmaMinus;
366  }
367  }
368  else
369  {
370  if (ran < 0.2)
371  {
372  pv[0] = AntiSigmaZero;
373  }
374  else if (ran < 0.4)
375  {
376  pv[0] = AntiSigmaPlus;
377  pv[1] = Proton;
378  }
379  else if (ran < 0.6)
380  {
381  pv[0] = Neutron;
382  pv[1] = AntiLambda;
383  }
384  else if (ran < 0.8)
385  {
386  pv[0] = Neutron;
387  pv[1] = AntiSigmaZero;
388  }
389  else
390  {
391  pv[0] = Proton;
392  pv[1] = AntiSigmaPlus;
393  }
394  }
395  }
396  return;
397  }
398  else if (availableEnergy <= PionPlus.getMass())
399  return;
400 
401  // inelastic scattering
402 
403  npos = 0; nneg = 0; nzero = 0;
404  G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88,
405  0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40,
406  0.39, 0.36, 0.33, 0.10, 0.01};
407  G4int iplab = G4int( incidentTotalMomentum*10.);
408  if ( iplab > 9) iplab = 10 + G4int( (incidentTotalMomentum -1.)*5. );
409  if ( iplab > 14) iplab = 15 + G4int( incidentTotalMomentum -2. );
410  if ( iplab > 22) iplab = 23 + G4int( (incidentTotalMomentum -10.)/10.);
411  iplab = std::min(24, iplab);
412 
413  if ( G4UniformRand() > anhl[iplab] )
414  { // non- annihilation channels
415 
416  // number of total particles vs. centre of mass Energy - 2*proton mass
417 
418  G4double aleab = std::log(availableEnergy);
419  G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
420  + aleab*(0.117712+0.0136912*aleab))) - 2.0;
421 
422  // normalization constant for kno-distribution.
423  // calculate first the sum of all constants, check for numerical problems.
424  G4double test, dum, anpn = 0.0;
425 
426  for (nt=1; nt<=numSec; nt++) {
427  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
428  dum = pi*nt/(2.0*n*n);
429  if (std::fabs(dum) < 1.0) {
430  if( test >= 1.0e-10 )anpn += dum*test;
431  } else {
432  anpn += dum*test;
433  }
434  }
435 
436  G4double ran = G4UniformRand();
437  G4double excs = 0.0;
438  if( targetCode == protonCode )
439  {
440  counter = -1;
441  for( npos=0; npos<numSec/3; npos++ )
442  {
443  for( nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++ )
444  {
445  for( nzero=0; nzero<numSec/3; nzero++ )
446  {
447  if( ++counter < numMul )
448  {
449  nt = npos+nneg+nzero;
450  if ( (nt>0) && (nt<=numSec) ) {
451  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
452  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
453  if (std::fabs(dum) < 1.0) {
454  if( test >= 1.0e-10 )excs += dum*test;
455  } else {
456  excs += dum*test;
457  }
458 
459  if (ran < excs) goto outOfLoop; //----------------------->
460  }
461  }
462  }
463  }
464  }
465 
466  // 3 previous loops continued to the end
467  inElastic = false; // quasi-elastic scattering
468  return;
469  }
470  else
471  { // target must be a neutron
472  counter = -1;
473  for( npos=0; npos<numSec/3; npos++ )
474  {
475  for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
476  {
477  for( nzero=0; nzero<numSec/3; nzero++ )
478  {
479  if( ++counter < numMul )
480  {
481  nt = npos+nneg+nzero;
482  if ( (nt>0) && (nt<=numSec) ) {
483  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
484  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
485  if (std::fabs(dum) < 1.0) {
486  if( test >= 1.0e-10 )excs += dum*test;
487  } else {
488  excs += dum*test;
489  }
490 
491  if (ran < excs) goto outOfLoop; // -------------------------->
492  }
493  }
494  }
495  }
496  }
497  // 3 previous loops continued to the end
498  inElastic = false; // quasi-elastic scattering.
499  return;
500  }
501 
502  outOfLoop: // <------------------------------------------------------------------------
503 
504  ran = G4UniformRand();
505 
506  if( targetCode == protonCode)
507  {
508  if( npos == nneg)
509  {
510  if (ran < 0.40)
511  {
512  }
513  else if (ran < 0.8)
514  {
515  pv[0] = AntiSigmaZero;
516  }
517  else
518  {
519  pv[0] = AntiSigmaMinus;
520  pv[1] = Neutron;
521  }
522  }
523  else if (npos == (nneg+1))
524  {
525  if( ran < 0.25)
526  {
527  pv[1] = Neutron;
528  }
529  else if (ran < 0.5)
530  {
531  pv[0] = AntiSigmaZero;
532  pv[1] = Neutron;
533  }
534  else
535  {
536  pv[0] = AntiSigmaPlus;
537  }
538  }
539  else if (npos == (nneg-1))
540  {
541  pv[0] = AntiSigmaMinus;
542  }
543  else
544  {
545  pv[0] = AntiSigmaPlus;
546  pv[1] = Neutron;
547  }
548  }
549  else
550  {
551  if( npos == nneg)
552  {
553  if (ran < 0.4)
554  {
555  }
556  else if(ran < 0.8)
557  {
558  pv[0] = AntiSigmaZero;
559  }
560  else
561  {
562  pv[0] = AntiSigmaPlus;
563  pv[1] = Proton;
564  }
565  }
566  else if ( npos == (nneg-1))
567  {
568  if (ran < 0.5)
569  {
570  pv[0] = AntiSigmaMinus;
571  }
572  else if (ran < 0.75)
573  {
574  pv[1] = Proton;
575  }
576  else
577  {
578  pv[0] = AntiSigmaZero;
579  pv[1] = Proton;
580  }
581  }
582  else if (npos == (nneg+1))
583  {
584  pv[0] = AntiSigmaPlus;
585  }
586  else
587  {
588  pv[0] = AntiSigmaMinus;
589  pv[1] = Proton;
590  }
591  }
592 
593  }
594  else // annihilation
595  {
596  if ( availableEnergy > 2. * PionPlus.getMass() )
597  {
598 
599  G4double aleab = std::log(availableEnergy);
600  G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
601  + aleab*(0.117712+0.0136912*aleab))) - 2.0;
602 
603  // normalization constant for kno-distribution.
604  // calculate first the sum of all constants, check for numerical problems.
605  G4double test, dum, anpn = 0.0;
606 
607  for (nt=2; nt<=numSec; nt++) {
608  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
609  dum = pi*nt/(2.0*n*n);
610  if (std::fabs(dum) < 1.0) {
611  if( test >= 1.0e-10 )anpn += dum*test;
612  } else {
613  anpn += dum*test;
614  }
615  }
616 
617  G4double ran = G4UniformRand();
618  G4double excs = 0.0;
619  if( targetCode == protonCode )
620  {
621  counter = -1;
622  for( npos=1; npos<numSec/3; npos++ )
623  {
624  nneg = npos-1;
625  for( nzero=0; nzero<numSec/3; nzero++ )
626  {
627  if( ++counter < numMulAn )
628  {
629  nt = npos+nneg+nzero;
630  if ( (nt>1) && (nt<=numSec) ) {
631  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
632  dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
633  if (std::fabs(dum) < 1.0) {
634  if( test >= 1.0e-10 )excs += dum*test;
635  } else {
636  excs += dum*test;
637  }
638 
639  if (ran < excs) goto outOfLoopAn; //----------------------->
640  }
641  }
642  }
643  }
644  // 3 previous loops continued to the end
645  inElastic = false; // quasi-elastic scattering
646  return;
647  }
648  else
649  { // target must be a neutron
650  counter = -1;
651  for( npos=0; npos<numSec/3; npos++ )
652  {
653  nneg = npos;
654  for( nzero=0; nzero<numSec/3; nzero++ )
655  {
656  if( ++counter < numMulAn )
657  {
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  if (std::fabs(dum) < 1.0) {
663  if( test >= 1.0e-10 )excs += dum*test;
664  } else {
665  excs += dum*test;
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].getCode() << " " ;
713  G4cout << pv[1].getCode() << " " ;
714  for (i=2; i < vecLen; i++)
715  {
716  G4cout << pv[i].getCode() << " " ;
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