source: trunk/source/processes/hadronic/models/high_energy/src/G4HEAntiProtonInelastic.cc @ 962

Last change on this file since 962 was 962, checked in by garnier, 15 years ago

update processes
File size: 27.4 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//
27// $Id: G4HEAntiProtonInelastic.cc,v 1.14 2008/03/17 20:49:17 dennis Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
29//
30//
31
32#include "globals.hh"
33#include "G4ios.hh"
34
35//
36// G4 Process: Gheisha High Energy Collision model.
37// This includes the high energy cascading model, the two-body-resonance model
38// and the low energy two-body model. Not included are the low energy stuff like
39// nuclear reactions, nuclear fission without any cascading and all processes for
40// particles at rest. 
41// First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96. 
42// H. Fesefeldt, RWTH-Aachen, 23-October-1996
43// Last modified: 29-July-1998
44 
45#include "G4HEAntiProtonInelastic.hh"
46
47G4HadFinalState *  G4HEAntiProtonInelastic::
48ApplyYourself( const G4HadProjectile &aTrack, G4Nucleus &targetNucleus )
49  {
50    G4HEVector * pv = new G4HEVector[MAXPART];
51    const G4HadProjectile *aParticle = &aTrack;
52//    G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
53    const G4double atomicWeight = targetNucleus.GetN();
54    const G4double atomicNumber  = targetNucleus.GetZ();
55    G4HEVector incidentParticle(aParticle);
56
57    G4int    incidentCode          = incidentParticle.getCode();
58    G4double incidentMass          = incidentParticle.getMass();
59    G4double incidentTotalEnergy   = incidentParticle.getEnergy();
60    G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
61    G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
62
63    if(incidentKineticEnergy < 1.)
64      { 
65        G4cout << "GHEAntiProtonInelastic: incident energy < 1 GeV" << G4endl;
66      }
67    if(verboseLevel > 1)
68      {
69        G4cout << "G4HEAntiProtonInelastic::ApplyYourself" << G4endl;
70        G4cout << "incident particle " << incidentParticle.getName()
71             << "mass "              << incidentMass
72             << "kinetic energy "    << incidentKineticEnergy
73             << G4endl;
74        G4cout << "target material with (A,Z) = (" 
75             << atomicWeight << "," << atomicNumber << ")" << G4endl;
76      }
77
78    G4double inelasticity  = NuclearInelasticity(incidentKineticEnergy, 
79                                                 atomicWeight, atomicNumber);
80    if(verboseLevel > 1)
81        G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
82   
83
84    incidentKineticEnergy -= inelasticity;
85   
86    G4double excitationEnergyGNP = 0.;
87    G4double excitationEnergyDTA = 0.; 
88
89    G4double excitation    = NuclearExcitation(incidentKineticEnergy,
90                                               atomicWeight, atomicNumber,
91                                               excitationEnergyGNP,
92                                               excitationEnergyDTA);
93    if(verboseLevel > 1)
94      G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
95           << excitationEnergyDTA << G4endl;             
96
97
98    incidentKineticEnergy -= excitation;
99    incidentTotalEnergy    = incidentKineticEnergy + incidentMass;
100    incidentTotalMomentum  = std::sqrt( (incidentTotalEnergy-incidentMass)                   
101                                  *(incidentTotalEnergy+incidentMass));
102
103    G4HEVector targetParticle;
104    if(G4UniformRand() < atomicNumber/atomicWeight)
105      { 
106        targetParticle.setDefinition("Proton");
107      }
108    else
109      { 
110        targetParticle.setDefinition("Neutron");
111      }
112
113    G4double targetMass         = targetParticle.getMass();
114    G4double centerOfMassEnergy = std::sqrt( incidentMass*incidentMass + targetMass*targetMass
115                                       + 2.0*targetMass*incidentTotalEnergy);
116    G4double availableEnergy    = centerOfMassEnergy - targetMass - incidentMass;
117
118                                                                // this was the meaning of inElastic in the
119                                                                // original Gheisha stand-alone version.
120//    G4bool   inElastic          = InElasticCrossSectionInFirstInt
121//                                    (availableEnergy, incidentCode, incidentTotalMomentum); 
122                                                                // by unknown reasons, it has been replaced
123                                                                // to the following code in Geant???
124    G4bool inElastic = true;
125//    if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false;   
126
127    vecLength  = 0;           
128       
129    if(verboseLevel > 1)
130      G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
131           << incidentCode << G4endl;
132
133    G4bool successful = false; 
134   
135    if(inElastic || (!inElastic && atomicWeight < 1.5))
136      { 
137        FirstIntInCasAntiProton(inElastic, availableEnergy, pv, vecLength,
138                                incidentParticle, targetParticle, atomicWeight);
139
140        if(verboseLevel > 1)
141           G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl; 
142
143
144        if ((vecLength > 0) && (availableEnergy > 1.)) 
145                   StrangeParticlePairProduction( availableEnergy, centerOfMassEnergy,
146                                                  pv, vecLength,
147                                                  incidentParticle, targetParticle);
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      }
174    if (!successful)
175      { 
176            ElasticScattering( successful, pv, vecLength,
177                               incidentParticle,   
178                               atomicWeight, atomicNumber);
179      }
180
181    if (!successful)
182      { 
183        G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" << G4endl;
184      }
185      FillParticleChange(pv,  vecLength);
186      delete [] pv;
187      theParticleChange.SetStatusChange(stopAndKill);
188      return & theParticleChange;
189  }
190
191void
192G4HEAntiProtonInelastic::FirstIntInCasAntiProton( G4bool &inElastic,
193                                                  const G4double availableEnergy,
194                                                  G4HEVector pv[],
195                                                  G4int &vecLen,
196                                                  G4HEVector incidentParticle,
197                                                  G4HEVector targetParticle,
198                                                  const G4double atomicWeight)
199
200// AntiProton undergoes interaction with nucleon within a nucleus.  Check if it is
201// energetically possible to produce pions/kaons.  In not, assume nuclear excitation
202// occurs and input particle is degraded in energy. No other particles are produced.
203// If reaction is possible, find the correct number of pions/protons/neutrons
204// produced using an interpolation to multiplicity data.  Replace some pions or
205// protons/neutrons by kaons or strange baryons according to the average
206// multiplicity per inelastic reaction.
207
208 {
209   static const G4double expxu =  std::log(MAXFLOAT); // upper bound for arg. of exp
210   static const G4double expxl = -expxu;         // lower bound for arg. of exp
211
212   static const G4double protb = 0.7;
213   static const G4double neutb = 0.7;
214   static const G4double     c = 1.25;
215
216   static const G4int   numMul   = 1200;
217   static const G4int   numMulAn = 400;
218   static const G4int   numSec   = 60;
219
220   G4int              neutronCode = Neutron.getCode();
221   G4int              protonCode  = Proton.getCode();
222
223   G4int               targetCode = targetParticle.getCode();
224//   G4double          incidentMass = incidentParticle.getMass();
225//   G4double        incidentEnergy = incidentParticle.getEnergy();
226   G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
227
228   static G4bool first = true;
229   static G4double protmul[numMul],  protnorm[numSec];   // proton constants
230   static G4double protmulAn[numMulAn],protnormAn[numSec]; 
231   static G4double neutmul[numMul],  neutnorm[numSec];   // neutron constants
232   static G4double neutmulAn[numMulAn],neutnormAn[numSec];
233
234                              //  misc. local variables
235                              //  np = number of pi+,  nm = number of pi-,  nz = number of pi0
236
237   G4int i, counter, nt, np, nm, nz;
238
239   if( first ) 
240     {                         // compute normalization constants, this will only be done once
241       first = false;
242       for( i=0; i<numMul  ; i++ ) protmul[i]  = 0.0;
243       for( i=0; i<numSec  ; i++ ) protnorm[i] = 0.0;
244       counter = -1;
245       for( np=0; np<(numSec/3); np++ ) 
246          {
247            for( nm=Imax(0,np-1); nm<=(np+1); nm++ ) 
248               {
249                 for( nz=0; nz<numSec/3; nz++ ) 
250                    {
251                      if( ++counter < numMul ) 
252                        {
253                          nt = np+nm+nz;
254                          if( (nt>0) && (nt<=numSec) ) 
255                            {
256                              protmul[counter] = pmltpc(np,nm,nz,nt,protb,c);
257                              protnorm[nt-1] += protmul[counter];
258                            }
259                        }
260                    }
261               }
262          }
263       for( i=0; i<numMul; i++ )neutmul[i]  = 0.0;
264
265       for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
266       counter = -1;
267       for( np=0; np<numSec/3; np++ ) 
268          {
269            for( nm=np; nm<=(np+2); nm++ ) 
270               {
271                 for( nz=0; nz<numSec/3; nz++ ) 
272                    {
273                      if( ++counter < numMul ) 
274                        {
275                          nt = np+nm+nz;
276                          if( (nt>0) && (nt<=numSec) ) 
277                            {
278                               neutmul[counter] = pmltpc(np,nm,nz,nt,neutb,c);
279                               neutnorm[nt-1] += neutmul[counter];
280                            }
281                        }
282                    }
283               }
284          }
285       for( i=0; i<numSec; i++ ) 
286          {
287            if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
288            if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
289          }
290                                                                   // annihilation
291       for( i=0; i<numMulAn  ; i++ ) protmulAn[i]  = 0.0;
292       for( i=0; i<numSec    ; i++ ) protnormAn[i] = 0.0;
293       counter = -1;
294       for( np=1; np<(numSec/3); np++ ) 
295          {
296            nm = np; 
297            for( nz=0; nz<numSec/3; nz++ ) 
298               {
299                 if( ++counter < numMulAn ) 
300                   {
301                     nt = np+nm+nz;
302                     if( (nt>0) && (nt<=numSec) ) 
303                       {
304                         protmulAn[counter] = pmltpc(np,nm,nz,nt,protb,c);
305                         protnormAn[nt-1] += protmulAn[counter];
306                       }
307                   }
308               }
309          }
310       for( i=0; i<numMulAn; i++ ) neutmulAn[i]  = 0.0;
311       for( i=0; i<numSec;   i++ ) neutnormAn[i] = 0.0;
312       counter = -1;
313       for( np=1; np<numSec/3; np++ ) 
314          {
315            nm = np+1; 
316            for( nz=0; nz<numSec/3; nz++ ) 
317               {
318                 if( ++counter < numMulAn ) 
319                   {
320                     nt = np+nm+nz;
321                     if( (nt>0) && (nt<=numSec) ) 
322                       {
323                          neutmulAn[counter] = pmltpc(np,nm,nz,nt,neutb,c);
324                          neutnormAn[nt-1] += neutmulAn[counter];
325                       }
326                   }
327               }
328          }
329       for( i=0; i<numSec; i++ ) 
330          {
331            if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
332            if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
333          }
334     }                                          // end of initialization
335
336         
337                                              // initialize the first two places
338                                              // the same as beam and target                                   
339   pv[0] = incidentParticle;
340   pv[1] = targetParticle;
341   vecLen = 2;
342
343   if( !inElastic ) 
344     {                                        // pb p --> nb n
345       if( targetCode == protonCode ) 
346         {
347           G4double cech[] = {0.14, 0.170, 0.180, 0.180, 0.180, 0.170, 0.170, 0.160, 0.155, 0.145,
348                              0.11, 0.082, 0.065, 0.050, 0.041, 0.035, 0.028, 0.024, 0.010, 0.000};
349
350           G4int iplab = G4int( incidentTotalMomentum*10.);
351           if (iplab > 9) iplab = Imin(19, G4int( incidentTotalMomentum) + 9);
352           if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) 
353             {                                            // charge exchange  pi+ n -> pi0 p
354               pv[0] = AntiNeutron;
355               pv[1] = Neutron;
356             }
357         }
358       return;
359     }
360   else if (availableEnergy <= PionPlus.getMass())
361       return;
362
363                                                  //   inelastic scattering
364
365   np = 0; nm = 0; nz = 0;
366   G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.90,
367                      0.60, 0.52, 0.47, 0.44, 0.41, 0.39, 0.37, 0.35, 0.34, 0.24,
368                      0.19, 0.15, 0.12, 0.10, 0.09, 0.07, 0.06, 0.05, 0.00};
369   G4int            iplab =      G4int( incidentTotalMomentum*10.);
370   if ( iplab >  9) iplab =  9 + G4int( incidentTotalMomentum);         
371   if ( iplab > 18) iplab = 18 + G4int( incidentTotalMomentum*10.);
372                    iplab = Imin(28, iplab);
373
374   if ( G4UniformRand() > anhl[iplab] )
375     {
376
377       G4double eab = availableEnergy;
378       G4int ieab = G4int( eab*5.0 );
379   
380       G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
381       if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) ) 
382         {
383                                              //   suppress high multiplicity events at low momentum
384                                              //   only one additional pion will be produced
385           G4double w0, wp, wm, wt, ran;
386           if( targetCode == neutronCode )                    // target is a neutron
387             {
388               w0 = - sqr(1.+neutb)/(2.*c*c);
389               w0 = std::exp(w0);
390               wm = - sqr(-1.+neutb)/(2.*c*c); 
391               wm = std::exp(wm);
392               if( G4UniformRand() < w0/(w0+wm) ) 
393                 { np = 0; nm = 0; nz = 1; }
394               else 
395                 { np = 0; nm = 1; nz = 0; }       
396             } 
397           else 
398             {                                               // target is a proton
399               w0 = -sqr(1.+protb)/(2.*c*c);
400               w0 = std::exp(w0);
401               wp = w0;
402               wm = -sqr(-1.+protb)/(2.*c*c);
403               wm = std::exp(wm);
404               wt = w0+wp+wm;
405               wp = w0+wp;
406               ran = G4UniformRand();
407               if( ran < w0/wt)
408                 { np = 0; nm = 0; nz = 1; }       
409               else if( ran < wp/wt)
410                 { np = 1; nm = 0; nz = 0; }       
411               else
412                 { np = 0; nm = 1; nz = 0; }       
413             }
414         }
415       else
416         {
417                         //  number of total particles vs. centre of mass Energy - 2*proton mass
418   
419           G4double aleab = std::log(availableEnergy);
420           G4double n     = 3.62567+aleab*(0.665843+aleab*(0.336514
421                            + aleab*(0.117712+0.0136912*aleab))) - 2.0;
422   
423                         // normalization constant for kno-distribution.
424                         // calculate first the sum of all constants, check for numerical problems.   
425           G4double test, dum, anpn = 0.0;
426
427           for (nt=1; nt<=numSec; nt++) {
428             test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
429             dum = pi*nt/(2.0*n*n);
430             if (std::fabs(dum) < 1.0) { 
431               if( test >= 1.0e-10 )anpn += dum*test;
432             } else { 
433               anpn += dum*test;
434             }
435           }
436   
437           G4double ran = G4UniformRand();
438           G4double excs = 0.0;
439           if( targetCode == protonCode ) 
440             {
441               counter = -1;
442               for( np=0; np<numSec/3; np++ ) 
443                  {
444                    for( nm=Imax(0,np-1); nm<=(np+1); nm++ ) 
445                       {
446                         for( nz=0; nz<numSec/3; nz++ ) 
447                            {
448                              if( ++counter < numMul ) 
449                                {
450                                  nt = np+nm+nz;
451                                  if ( (nt>0) && (nt<=numSec) ) {
452                                    test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
453                                    dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
454                                    if (std::fabs(dum) < 1.0) { 
455                                      if( test >= 1.0e-10 )excs += dum*test;
456                                    } else { 
457                                      excs += dum*test;
458                                    }
459
460                                    if (ran < excs) goto outOfLoop;      //----------------------->
461                                  }   
462                                }   
463                            }     
464                       }                                                                                 
465                  }
466       
467                                              // 3 previous loops continued to the end
468               inElastic = false;                 // quasi-elastic scattering   
469               return;
470             }
471           else   
472             {                                         // target must be a neutron
473               counter = -1;
474               for( np=0; np<numSec/3; np++ ) 
475                  {
476                    for( nm=np; nm<=(np+2); nm++ ) 
477                       {
478                         for( nz=0; nz<numSec/3; nz++ ) 
479                            {
480                              if( ++counter < numMul ) 
481                                {
482                                  nt = np+nm+nz;
483                                  if ( (nt>=1) && (nt<=numSec) ) {
484                                    test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
485                                    dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
486                                    if (std::fabs(dum) < 1.0) { 
487                                      if( test >= 1.0e-10 )excs += dum*test;
488                                    } else { 
489                                      excs += dum*test;
490                                    }
491
492                                    if (ran < excs) goto outOfLoop;       // -------------------------->
493                                  }
494                                }
495                            }
496                       }
497                  }
498                                                      // 3 previous loops continued to the end
499               inElastic = false;                     // quasi-elastic scattering.
500               return;
501             }
502         } 
503       outOfLoop:           //  <------------------------------------------------------------------------   
504   
505       if( targetCode == neutronCode)
506         {
507           if( np == nm)
508             {
509             }
510           else if (np == (nm-1))
511             {
512               if( G4UniformRand() < 0.5)
513                 {
514                   pv[1] = Proton;
515                 }
516               else
517                 {
518                   pv[0] = AntiNeutron;
519                 }
520             }
521           else     
522             {
523               pv[0] = AntiNeutron;
524               pv[1] = Proton;
525             } 
526         } 
527       else
528         {
529           if( np == nm)
530             {
531               if( G4UniformRand() < 0.25)
532                 {
533                   pv[0] = AntiNeutron;
534                   pv[1] = Neutron;
535                 }
536               else
537                 {
538                 }
539             } 
540           else if ( np == (1+nm))
541             {
542               pv[1] = Neutron;
543             }
544           else
545             {
546               pv[0] = AntiNeutron;
547             }
548         }     
549     
550     }
551   else                                                               // annihilation
552     { 
553       if ( availableEnergy > 2. * PionPlus.getMass() )
554         {
555
556           G4double aleab = std::log(availableEnergy);
557           G4double n     = 3.62567+aleab*(0.665843+aleab*(0.336514
558                            + aleab*(0.117712+0.0136912*aleab))) - 2.0;
559   
560           // normalization constant for kno-distribution.
561           // calculate first the sum of all constants, check for numerical problems.   
562           G4double test, dum, anpn = 0.0;
563
564           for (nt=2; nt<=numSec; nt++) {
565             test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
566             dum = pi*nt/(2.0*n*n);
567             if (std::fabs(dum) < 1.0) { 
568               if( test >= 1.0e-10 )anpn += dum*test;
569             } else { 
570               anpn += dum*test;
571             }
572           }
573   
574           G4double ran = G4UniformRand();
575           G4double excs = 0.0;
576           if( targetCode == protonCode ) 
577             {
578               counter = -1;
579               for( np=1; np<numSec/3; np++ ) 
580                  {
581                    nm = np; 
582                    for( nz=0; nz<numSec/3; nz++ ) 
583                      {
584                        if( ++counter < numMulAn ) 
585                          {
586                            nt = np+nm+nz;
587                            if ( (nt>0) && (nt<=numSec) ) {
588                              test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
589                              dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
590                              if (std::fabs(dum) < 1.0) { 
591                                if( test >= 1.0e-10 )excs += dum*test;
592                              } else { 
593                                excs += dum*test;
594                              }
595
596                              if (ran < excs) goto outOfLoopAn;      //----------------------->
597                            } 
598                          }   
599                      }     
600                 }                                                                                 
601                                              // 3 previous loops continued to the end
602               inElastic = false;                 // quasi-elastic scattering   
603               return;
604             }
605           else   
606             {                                         // target must be a neutron
607               counter = -1;
608               for( np=1; np<numSec/3; np++ ) 
609                 { 
610                   nm = np+1; 
611                   for( nz=0; nz<numSec/3; nz++ ) 
612                      {
613                        if( ++counter < numMulAn ) 
614                          {
615                            nt = np+nm+nz;
616                            if ( (nt>=1) && (nt<=numSec) ) {
617                              test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
618                              dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
619                              if (std::fabs(dum) < 1.0) { 
620                                if( test >= 1.0e-10 )excs += dum*test;
621                              } else { 
622                                excs += dum*test;
623                              }
624
625                              if (ran < excs) goto outOfLoopAn;       // -------------------------->
626                            }
627                          }
628                      }
629                 }
630               inElastic = false;                     // quasi-elastic scattering.
631               return;
632             } 
633       outOfLoopAn:           //  <------------------------------------------------------------------   
634       vecLen = 0;
635         }
636     }
637
638   nt = np + nm + nz;
639   while ( nt > 0)
640       {
641         G4double ran = G4UniformRand();
642         if ( ran < (G4double)np/nt)
643            { 
644              if( np > 0 ) 
645                { pv[vecLen++] = PionPlus;
646                  np--;
647                }
648            }
649         else if ( ran < (G4double)(np+nm)/nt)
650            {   
651              if( nm > 0 )
652                { 
653                  pv[vecLen++] = PionMinus;
654                  nm--;
655                }
656            }
657         else
658            {
659              if( nz > 0 )
660                { 
661                  pv[vecLen++] = PionZero;
662                  nz--;
663                }
664            }
665         nt = np + nm + nz;
666       } 
667   if (verboseLevel > 1)
668      {
669        G4cout << "Particles produced: " ;
670        G4cout << pv[0].getName() << " " ;
671        G4cout << pv[1].getName() << " " ;
672        for (i=2; i < vecLen; i++)   
673            { 
674              G4cout << pv[i].getName() << " " ;
675            }
676         G4cout << G4endl;
677      }
678   return;
679 }
680
681
682
683
684
685
686
687
688
Note: See TracBrowser for help on using the repository browser.