G4HEAntiProtonInelastic.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 
00029 // G4 Process: Gheisha High Energy Collision model.
00030 // This includes the high energy cascading model, the two-body-resonance model
00031 // and the low energy two-body model. Not included are the low energy stuff
00032 // like nuclear reactions, nuclear fission without any cascading and all
00033 // processes for particles at rest.  
00034 // First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96.  
00035 // H. Fesefeldt, RWTH-Aachen, 23-October-1996
00036  
00037 #include "G4HEAntiProtonInelastic.hh"
00038 #include "globals.hh"
00039 #include "G4ios.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 
00043 void G4HEAntiProtonInelastic::ModelDescription(std::ostream& outFile) const
00044 {
00045   outFile << "G4HEAntiProtonInelastic is one of the High Energy\n"
00046           << "Parameterized (HEP) models used to implement inelastic\n"
00047           << "anti-proton scattering from nuclei.  It is a re-engineered\n"
00048           << "version of the GHEISHA code of H. Fesefeldt.  It divides the\n"
00049           << "initial collision products into backward- and forward-going\n"
00050           << "clusters which are then decayed into final state hadrons.\n"
00051           << "The model does not conserve energy on an event-by-event\n"
00052           << "basis.  It may be applied to anti-protons with initial\n"
00053           << "energies above 20 GeV.\n";
00054 }
00055 
00056 
00057 G4HadFinalState*
00058 G4HEAntiProtonInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00059                                        G4Nucleus& targetNucleus)
00060 {
00061   G4HEVector* pv = new G4HEVector[MAXPART];
00062   const G4HadProjectile* aParticle = &aTrack;
00063   const G4double atomicWeight = targetNucleus.GetA_asInt();
00064   const G4double atomicNumber = targetNucleus.GetZ_asInt();
00065   G4HEVector incidentParticle(aParticle);
00066 
00067   G4int incidentCode = incidentParticle.getCode();
00068   G4double incidentMass = incidentParticle.getMass();
00069   G4double incidentTotalEnergy = incidentParticle.getEnergy();
00070 
00071   // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00072   // DHW 19 may 2011: variable set but not used
00073 
00074   G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
00075 
00076   if (incidentKineticEnergy < 1.)
00077     G4cout << "GHEAntiProtonInelastic: incident energy < 1 GeV" << G4endl;
00078 
00079   if (verboseLevel > 1) {
00080     G4cout << "G4HEAntiProtonInelastic::ApplyYourself" << G4endl;
00081     G4cout << "incident particle " << incidentParticle.getName()
00082            << "mass "              << incidentMass
00083            << "kinetic energy "    << incidentKineticEnergy
00084            << G4endl;
00085     G4cout << "target material with (A,Z) = (" 
00086            << atomicWeight << "," << atomicNumber << ")" << G4endl;
00087   }
00088 
00089   G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 
00090                                               atomicWeight, atomicNumber);
00091   if (verboseLevel > 1)
00092     G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
00093 
00094   incidentKineticEnergy -= inelasticity;
00095     
00096   G4double excitationEnergyGNP = 0.;
00097   G4double excitationEnergyDTA = 0.; 
00098 
00099   G4double excitation = NuclearExcitation(incidentKineticEnergy,
00100                                           atomicWeight, atomicNumber,
00101                                           excitationEnergyGNP,
00102                                           excitationEnergyDTA);
00103   if (verboseLevel > 1)
00104     G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 
00105            << excitationEnergyDTA << G4endl;             
00106 
00107 
00108   incidentKineticEnergy -= excitation;
00109   incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00110   // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
00111   //                                   *(incidentTotalEnergy+incidentMass));
00112   // DHW 19 may 2011: variable set but not used
00113 
00114   G4HEVector targetParticle;
00115   if (G4UniformRand() < atomicNumber/atomicWeight) { 
00116     targetParticle.setDefinition("Proton");
00117   } else { 
00118     targetParticle.setDefinition("Neutron");
00119   }
00120 
00121   G4double targetMass = targetParticle.getMass();
00122   G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
00123                                         + targetMass*targetMass
00124                                         + 2.0*targetMass*incidentTotalEnergy);
00125   G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00126 
00127   G4bool inElastic = true;
00128   vecLength = 0;   
00129         
00130   if (verboseLevel > 1)
00131     G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00132            << incidentCode << G4endl;
00133 
00134   G4bool successful = false; 
00135     
00136   FirstIntInCasAntiProton(inElastic, availableEnergy, pv, vecLength,
00137                           incidentParticle, targetParticle, atomicWeight);
00138 
00139   if (verboseLevel > 1)
00140     G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;  
00141 
00142   if ((vecLength > 0) && (availableEnergy > 1.)) 
00143     StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00144                                   pv, vecLength,
00145                                   incidentParticle, targetParticle);
00146   HighEnergyCascading(successful, pv, vecLength,
00147                       excitationEnergyGNP, excitationEnergyDTA,
00148                       incidentParticle, targetParticle,
00149                       atomicWeight, atomicNumber);
00150   if (!successful)
00151     HighEnergyClusterProduction(successful, pv, vecLength,
00152                                 excitationEnergyGNP, excitationEnergyDTA,
00153                                 incidentParticle, targetParticle,
00154                                 atomicWeight, atomicNumber);
00155   if (!successful) 
00156     MediumEnergyCascading(successful, pv, vecLength, 
00157                           excitationEnergyGNP, excitationEnergyDTA, 
00158                           incidentParticle, targetParticle,
00159                           atomicWeight, atomicNumber);
00160 
00161   if (!successful)
00162     MediumEnergyClusterProduction(successful, pv, vecLength,
00163                                   excitationEnergyGNP, excitationEnergyDTA,       
00164                                   incidentParticle, targetParticle,
00165                                   atomicWeight, atomicNumber);
00166   if (!successful)
00167     QuasiElasticScattering(successful, pv, vecLength,
00168                            excitationEnergyGNP, excitationEnergyDTA,
00169                            incidentParticle, targetParticle, 
00170                            atomicWeight, atomicNumber);
00171   if (!successful)
00172     ElasticScattering(successful, pv, vecLength,
00173                       incidentParticle,    
00174                       atomicWeight, atomicNumber);
00175 
00176   if (!successful)
00177     G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 
00178            << G4endl;
00179 
00180   FillParticleChange(pv,  vecLength);
00181   delete [] pv;
00182   theParticleChange.SetStatusChange(stopAndKill);
00183   return & theParticleChange;
00184 }
00185 
00186 
00187 void
00188 G4HEAntiProtonInelastic::FirstIntInCasAntiProton(G4bool& inElastic,
00189                                                  const G4double availableEnergy,
00190                                                  G4HEVector pv[],
00191                                                  G4int& vecLen,
00192                                                  const G4HEVector& incidentParticle,
00193                                                  const G4HEVector& targetParticle,
00194                                                  const G4double atomicWeight)
00195 
00196 // AntiProton undergoes interaction with nucleon within a nucleus.  Check if it is
00197 // energetically possible to produce pions/kaons.  In not, assume nuclear excitation
00198 // occurs and input particle is degraded in energy. No other particles are produced.
00199 // If reaction is possible, find the correct number of pions/protons/neutrons
00200 // produced using an interpolation to multiplicity data.  Replace some pions or
00201 // protons/neutrons by kaons or strange baryons according to the average
00202 // multiplicity per inelastic reaction.
00203 {
00204   static const G4double expxu = 82;      // upper bound for arg. of exp
00205   static const G4double expxl = -expxu;  // lower bound for arg. of exp
00206 
00207   static const G4double protb = 0.7;
00208   static const G4double neutb = 0.7;
00209   static const G4double     c = 1.25;
00210 
00211   static const G4int numMul   = 1200;
00212   static const G4int numMulAn = 400;
00213   static const G4int numSec   = 60;
00214 
00215   G4int neutronCode = Neutron.getCode();
00216   G4int protonCode  = Proton.getCode();
00217 
00218   G4int targetCode = targetParticle.getCode();
00219   G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00220 
00221   static G4bool first = true;
00222   static G4double protmul[numMul],  protnorm[numSec];   // proton constants
00223   static G4double protmulAn[numMulAn],protnormAn[numSec]; 
00224   static G4double neutmul[numMul],  neutnorm[numSec];   // neutron constants
00225   static G4double neutmulAn[numMulAn],neutnormAn[numSec];
00226 
00227   //  misc. local variables
00228   //  npos = number of pi+,  nneg = number of pi-,  nzero = number of pi0
00229 
00230   G4int i, counter, nt, npos, nneg, nzero;
00231 
00232    if( first ) 
00233      {             // compute normalization constants, this will only be done once
00234        first = false;
00235        for( i=0; i<numMul  ; i++ ) protmul[i]  = 0.0;
00236        for( i=0; i<numSec  ; i++ ) protnorm[i] = 0.0;
00237        counter = -1;
00238        for( npos=0; npos<(numSec/3); npos++ ) 
00239           {
00240             for( nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++ ) 
00241                {
00242                  for( nzero=0; nzero<numSec/3; nzero++ ) 
00243                     {
00244                       if( ++counter < numMul ) 
00245                         {
00246                           nt = npos+nneg+nzero;
00247                           if( (nt>0) && (nt<=numSec) ) 
00248                             {
00249                               protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00250                               protnorm[nt-1] += protmul[counter];
00251                             }
00252                         }
00253                     }
00254                }
00255           }
00256        for( i=0; i<numMul; i++ )neutmul[i]  = 0.0;
00257 
00258        for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00259        counter = -1;
00260        for( npos=0; npos<numSec/3; npos++ ) 
00261           {
00262             for( nneg=npos; nneg<=(npos+2); nneg++ ) 
00263                {
00264                  for( nzero=0; nzero<numSec/3; nzero++ ) 
00265                     {
00266                       if( ++counter < numMul ) 
00267                         {
00268                           nt = npos+nneg+nzero;
00269                           if( (nt>0) && (nt<=numSec) ) 
00270                             {
00271                                neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00272                                neutnorm[nt-1] += neutmul[counter];
00273                             }
00274                         }
00275                     }
00276                }
00277           }
00278        for( i=0; i<numSec; i++ ) 
00279           {
00280             if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00281             if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00282           }
00283                                                                    // annihilation
00284        for( i=0; i<numMulAn  ; i++ ) protmulAn[i]  = 0.0;
00285        for( i=0; i<numSec    ; i++ ) protnormAn[i] = 0.0;
00286        counter = -1;
00287        for( npos=1; npos<(numSec/3); npos++ ) 
00288           {
00289             nneg = npos; 
00290             for( nzero=0; nzero<numSec/3; nzero++ ) 
00291                {
00292                  if( ++counter < numMulAn ) 
00293                    {
00294                      nt = npos+nneg+nzero;
00295                      if( (nt>0) && (nt<=numSec) ) 
00296                        {
00297                          protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00298                          protnormAn[nt-1] += protmulAn[counter];
00299                        }
00300                    }
00301                }
00302           }
00303        for( i=0; i<numMulAn; i++ ) neutmulAn[i]  = 0.0;
00304        for( i=0; i<numSec;   i++ ) neutnormAn[i] = 0.0;
00305        counter = -1;
00306        for( npos=1; npos<numSec/3; npos++ ) 
00307           {
00308             nneg = npos+1; 
00309             for( nzero=0; nzero<numSec/3; nzero++ ) 
00310                {
00311                  if( ++counter < numMulAn ) 
00312                    {
00313                      nt = npos+nneg+nzero;
00314                      if( (nt>0) && (nt<=numSec) ) 
00315                        {
00316                           neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00317                           neutnormAn[nt-1] += neutmulAn[counter];
00318                        }
00319                    }
00320                }
00321           }
00322        for( i=0; i<numSec; i++ ) 
00323           {
00324             if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
00325             if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
00326           }
00327      }                                          // end of initialization
00328 
00329          
00330                                               // initialize the first two places
00331                                               // the same as beam and target                                    
00332    pv[0] = incidentParticle;
00333    pv[1] = targetParticle;
00334    vecLen = 2;
00335 
00336    if( !inElastic ) 
00337      {                                        // pb p --> nb n
00338        if( targetCode == protonCode ) 
00339          {
00340            G4double cech[] = {0.14, 0.170, 0.180, 0.180, 0.180, 0.170, 0.170, 0.160, 0.155, 0.145,
00341                               0.11, 0.082, 0.065, 0.050, 0.041, 0.035, 0.028, 0.024, 0.010, 0.000};
00342 
00343            G4int iplab = G4int( incidentTotalMomentum*10.);
00344            if (iplab > 9) iplab = Imin(19, G4int( incidentTotalMomentum) + 9);
00345            if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) 
00346              {                                            // charge exchange  pi+ n -> pi0 p
00347                pv[0] = AntiNeutron;
00348                pv[1] = Neutron;
00349              }
00350          }
00351        return;
00352      }
00353    else if (availableEnergy <= PionPlus.getMass())
00354        return;
00355 
00356                                                   //   inelastic scattering
00357 
00358    npos = 0; nneg = 0; nzero = 0;
00359    G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.90,
00360                       0.60, 0.52, 0.47, 0.44, 0.41, 0.39, 0.37, 0.35, 0.34, 0.24,
00361                       0.19, 0.15, 0.12, 0.10, 0.09, 0.07, 0.06, 0.05, 0.00};
00362    G4int            iplab =      G4int( incidentTotalMomentum*10.);
00363    if ( iplab >  9) iplab =  9 + G4int( incidentTotalMomentum);          
00364    if ( iplab > 18) iplab = 18 + G4int( incidentTotalMomentum*10.);
00365                     iplab = Imin(28, iplab);
00366 
00367    if ( G4UniformRand() > anhl[iplab] )
00368      {
00369 
00370        G4double eab = availableEnergy;
00371        G4int ieab = G4int( eab*5.0 );
00372    
00373        G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
00374        if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) ) 
00375          {
00376                                               //   suppress high multiplicity events at low momentum
00377                                               //   only one additional pion will be produced
00378            G4double w0, wp, wm, wt, ran;
00379            if( targetCode == neutronCode )                    // target is a neutron 
00380              {
00381                w0 = - sqr(1.+neutb)/(2.*c*c);
00382                w0 = std::exp(w0);
00383                wm = - sqr(-1.+neutb)/(2.*c*c);  
00384                wm = std::exp(wm);
00385                if( G4UniformRand() < w0/(w0+wm) )  
00386                  { npos = 0; nneg = 0; nzero = 1; }
00387                else 
00388                  { npos = 0; nneg = 1; nzero = 0; }       
00389              } 
00390            else 
00391              {                                               // target is a proton
00392                w0 = -sqr(1.+protb)/(2.*c*c);
00393                w0 = std::exp(w0);
00394                wp = w0;
00395                wm = -sqr(-1.+protb)/(2.*c*c);
00396                wm = std::exp(wm);
00397                wt = w0+wp+wm;
00398                wp = w0+wp;
00399                ran = G4UniformRand();
00400                if( ran < w0/wt)
00401                  { npos = 0; nneg = 0; nzero = 1; }       
00402                else if( ran < wp/wt)
00403                  { npos = 1; nneg = 0; nzero = 0; }       
00404                else
00405                  { npos = 0; nneg = 1; nzero = 0; }       
00406              }
00407          }
00408        else
00409          {
00410                          //  number of total particles vs. centre of mass Energy - 2*proton mass
00411    
00412            G4double aleab = std::log(availableEnergy);
00413            G4double n     = 3.62567+aleab*(0.665843+aleab*(0.336514
00414                             + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00415    
00416                          // normalization constant for kno-distribution.
00417                          // calculate first the sum of all constants, check for numerical problems.   
00418            G4double test, dum, anpn = 0.0;
00419 
00420            for (nt=1; nt<=numSec; nt++) {
00421              test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00422              dum = pi*nt/(2.0*n*n);
00423              if (std::fabs(dum) < 1.0) { 
00424                if( test >= 1.0e-10 )anpn += dum*test;
00425              } else { 
00426                anpn += dum*test;
00427              }
00428            }
00429    
00430            G4double ran = G4UniformRand();
00431            G4double excs = 0.0;
00432            if( targetCode == protonCode ) 
00433              {
00434                counter = -1;
00435                for( npos=0; npos<numSec/3; npos++ ) 
00436                   {
00437                     for( nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++ ) 
00438                        {
00439                          for( nzero=0; nzero<numSec/3; nzero++ ) 
00440                             {
00441                               if( ++counter < numMul ) 
00442                                 {
00443                                   nt = npos+nneg+nzero;
00444                                   if ( (nt>0) && (nt<=numSec) ) {
00445                                     test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00446                                     dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00447                                     if (std::fabs(dum) < 1.0) { 
00448                                       if( test >= 1.0e-10 )excs += dum*test;
00449                                     } else { 
00450                                       excs += dum*test;
00451                                     }
00452 
00453                                     if (ran < excs) goto outOfLoop;      //----------------------->
00454                                   }   
00455                                 }    
00456                             }     
00457                        }                                                                                  
00458                   }
00459        
00460                                               // 3 previous loops continued to the end
00461                inElastic = false;                 // quasi-elastic scattering   
00462                return;
00463              }
00464            else   
00465              {                                         // target must be a neutron
00466                counter = -1;
00467                for( npos=0; npos<numSec/3; npos++ ) 
00468                   {
00469                     for( nneg=npos; nneg<=(npos+2); nneg++ ) 
00470                        {
00471                          for( nzero=0; nzero<numSec/3; nzero++ ) 
00472                             {
00473                               if( ++counter < numMul ) 
00474                                 {
00475                                   nt = npos+nneg+nzero;
00476                                   if ( (nt>=1) && (nt<=numSec) ) {
00477                                     test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00478                                     dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00479                                     if (std::fabs(dum) < 1.0) { 
00480                                       if( test >= 1.0e-10 )excs += dum*test;
00481                                     } else { 
00482                                       excs += dum*test;
00483                                     }
00484 
00485                                     if (ran < excs) goto outOfLoop;       // -------------------------->
00486                                   }
00487                                 }
00488                             }
00489                        }
00490                   }
00491                                                       // 3 previous loops continued to the end
00492                inElastic = false;                     // quasi-elastic scattering.
00493                return;
00494              }
00495          } 
00496        outOfLoop:           //  <------------------------------------------------------------------------   
00497     
00498        if( targetCode == neutronCode)
00499          {
00500            if( npos == nneg)
00501              {
00502              }
00503            else if (npos == (nneg-1))
00504              {
00505                if( G4UniformRand() < 0.5)
00506                  {
00507                    pv[1] = Proton;
00508                  }
00509                else
00510                  {
00511                    pv[0] = AntiNeutron;
00512                  }
00513              }
00514            else      
00515              {
00516                pv[0] = AntiNeutron;
00517                pv[1] = Proton;
00518              } 
00519          }  
00520        else
00521          {
00522            if( npos == nneg)
00523              {
00524                if( G4UniformRand() < 0.25)
00525                  {
00526                    pv[0] = AntiNeutron;
00527                    pv[1] = Neutron;
00528                  }
00529                else
00530                  {
00531                  }
00532              } 
00533            else if ( npos == (1+nneg))
00534              {
00535                pv[1] = Neutron;
00536              }
00537            else
00538              {
00539                pv[0] = AntiNeutron;
00540              }
00541          }      
00542      
00543      }
00544    else                                                               // annihilation
00545      {  
00546        if ( availableEnergy > 2. * PionPlus.getMass() )
00547          {
00548 
00549            G4double aleab = std::log(availableEnergy);
00550            G4double n     = 3.62567+aleab*(0.665843+aleab*(0.336514
00551                             + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00552    
00553            // normalization constant for kno-distribution.
00554            // calculate first the sum of all constants, check for numerical problems.   
00555            G4double test, dum, anpn = 0.0;
00556 
00557            for (nt=2; nt<=numSec; nt++) {
00558              test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00559              dum = pi*nt/(2.0*n*n);
00560              if (std::fabs(dum) < 1.0) { 
00561                if( test >= 1.0e-10 )anpn += dum*test;
00562              } else { 
00563                anpn += dum*test;
00564              }
00565            }
00566    
00567            G4double ran = G4UniformRand();
00568            G4double excs = 0.0;
00569            if( targetCode == protonCode ) 
00570              {
00571                counter = -1;
00572                for( npos=1; npos<numSec/3; npos++ ) 
00573                   {
00574                     nneg = npos; 
00575                     for( nzero=0; nzero<numSec/3; nzero++ ) 
00576                       {
00577                         if( ++counter < numMulAn ) 
00578                           {
00579                             nt = npos+nneg+nzero;
00580                             if ( (nt>0) && (nt<=numSec) ) {
00581                               test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00582                               dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
00583                               if (std::fabs(dum) < 1.0) { 
00584                                 if( test >= 1.0e-10 )excs += dum*test;
00585                               } else { 
00586                                 excs += dum*test;
00587                               }
00588 
00589                               if (ran < excs) goto outOfLoopAn;      //----------------------->
00590                             }  
00591                           }    
00592                       }     
00593                  }                                                                                  
00594                                               // 3 previous loops continued to the end
00595                inElastic = false;                 // quasi-elastic scattering   
00596                return;
00597              }
00598            else   
00599              {                                         // target must be a neutron
00600                counter = -1;
00601                for( npos=1; npos<numSec/3; npos++ ) 
00602                  { 
00603                    nneg = npos+1; 
00604                    for( nzero=0; nzero<numSec/3; nzero++ ) 
00605                       {
00606                         if( ++counter < numMulAn ) 
00607                           {
00608                             nt = npos+nneg+nzero;
00609                             if ( (nt>=1) && (nt<=numSec) ) {
00610                               test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00611                               dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
00612                               if (std::fabs(dum) < 1.0) { 
00613                                 if( test >= 1.0e-10 )excs += dum*test;
00614                               } else { 
00615                                 excs += dum*test;
00616                               }
00617 
00618                               if (ran < excs) goto outOfLoopAn;       // -------------------------->
00619                             }
00620                           }
00621                       }
00622                  }
00623                inElastic = false;                     // quasi-elastic scattering.
00624                return;
00625              } 
00626        outOfLoopAn:           //  <------------------------------------------------------------------   
00627        vecLen = 0;
00628          }
00629      }
00630 
00631    nt = npos + nneg + nzero;
00632    while ( nt > 0)
00633        {
00634          G4double ran = G4UniformRand();
00635          if ( ran < (G4double)npos/nt)
00636             { 
00637               if( npos > 0 ) 
00638                 { pv[vecLen++] = PionPlus;
00639                   npos--;
00640                 }
00641             }
00642          else if ( ran < (G4double)(npos+nneg)/nt)
00643             {   
00644               if( nneg > 0 )
00645                 { 
00646                   pv[vecLen++] = PionMinus;
00647                   nneg--;
00648                 }
00649             }
00650          else
00651             {
00652               if( nzero > 0 )
00653                 { 
00654                   pv[vecLen++] = PionZero;
00655                   nzero--;
00656                 }
00657             }
00658          nt = npos + nneg + nzero;
00659        } 
00660    if (verboseLevel > 1)
00661       {
00662         G4cout << "Particles produced: " ;
00663         G4cout << pv[0].getName() << " " ;
00664         G4cout << pv[1].getName() << " " ;
00665         for (i=2; i < vecLen; i++)   
00666             { 
00667               G4cout << pv[i].getName() << " " ;
00668             }
00669          G4cout << G4endl;
00670       }
00671    return;
00672  }
00673 
00674 
00675 
00676 
00677 
00678 
00679 
00680 
00681 

Generated on Mon May 27 17:48:28 2013 for Geant4 by  doxygen 1.4.7