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