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