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