source: trunk/source/processes/hadronic/models/high_energy/src/G4HEAntiOmegaMinusInelastic.cc @ 1348

Last change on this file since 1348 was 1347, checked in by garnier, 14 years ago

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