G4HEAntiSigmaMinusInelastic.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 "G4HEAntiSigmaMinusInelastic.hh"
00038 #include "globals.hh"
00039 #include "G4ios.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 
00043 void G4HEAntiSigmaMinusInelastic::ModelDescription(std::ostream& outFile) const
00044 {
00045   outFile << "G4HEAntiSigmaMinusInelastic is one of the High Energy\n"
00046           << "Parameterized (HEP) models used to implement inelastic\n"
00047           << "anti-Sigma- 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-Sigma- with initial\n"
00053           << "energies above 20 GeV.\n";
00054 }
00055 
00056 
00057 G4HadFinalState*
00058 G4HEAntiSigmaMinusInelastic::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 << "GHEAntiSigmaMinusInelastic: incident energy < 1 GeV" << G4endl;
00078 
00079   if (verboseLevel > 1) {
00080     G4cout << "G4HEAntiSigmaMinusInelastic::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   incidentKineticEnergy -= excitation;
00108   incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00109   // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
00110   //                                   *(incidentTotalEnergy+incidentMass));
00111   // DHW 19 May 2011: variable set but not used
00112 
00113   G4HEVector targetParticle;
00114   if (G4UniformRand() < atomicNumber/atomicWeight) { 
00115     targetParticle.setDefinition("Proton");
00116   } else { 
00117     targetParticle.setDefinition("Neutron");
00118   }
00119 
00120   G4double targetMass = targetParticle.getMass();
00121   G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
00122                                         + targetMass*targetMass
00123                                         + 2.0*targetMass*incidentTotalEnergy);
00124   G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00125 
00126   G4bool inElastic = true;
00127   vecLength  = 0;           
00128         
00129   if (verboseLevel > 1)
00130     G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00131            << incidentCode << G4endl;
00132 
00133   G4bool successful = false; 
00134     
00135   FirstIntInCasAntiSigmaMinus(inElastic, availableEnergy, pv, vecLength,
00136                               incidentParticle, targetParticle, atomicWeight);
00137 
00138   if (verboseLevel > 1)
00139     G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00140 
00141   if ((vecLength > 0) && (availableEnergy > 1.)) 
00142     StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00143                                   pv, vecLength,
00144                                   incidentParticle, targetParticle);
00145   HighEnergyCascading(successful, pv, vecLength,
00146                       excitationEnergyGNP, excitationEnergyDTA,
00147                       incidentParticle, targetParticle,
00148                       atomicWeight, atomicNumber);
00149   if (!successful)
00150     HighEnergyClusterProduction(successful, pv, vecLength,
00151                                 excitationEnergyGNP, excitationEnergyDTA,
00152                                 incidentParticle, targetParticle,
00153                                 atomicWeight, atomicNumber);
00154   if (!successful) 
00155     MediumEnergyCascading(successful, pv, vecLength, 
00156                           excitationEnergyGNP, excitationEnergyDTA, 
00157                           incidentParticle, targetParticle,
00158                           atomicWeight, atomicNumber);
00159 
00160   if (!successful)
00161     MediumEnergyClusterProduction(successful, pv, vecLength,
00162                                   excitationEnergyGNP, excitationEnergyDTA,       
00163                                   incidentParticle, targetParticle,
00164                                   atomicWeight, atomicNumber);
00165   if (!successful)
00166     QuasiElasticScattering(successful, pv, vecLength,
00167                            excitationEnergyGNP, excitationEnergyDTA,
00168                            incidentParticle, targetParticle, 
00169                            atomicWeight, atomicNumber);
00170   if (!successful)
00171     ElasticScattering(successful, pv, vecLength,
00172                       incidentParticle,    
00173                       atomicWeight, atomicNumber);
00174 
00175   if (!successful)
00176     G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 
00177            << G4endl;
00178   FillParticleChange(pv,  vecLength);
00179   delete [] pv;
00180   theParticleChange.SetStatusChange(stopAndKill);
00181   return &theParticleChange;
00182 }
00183 
00184 
00185 void
00186 G4HEAntiSigmaMinusInelastic::FirstIntInCasAntiSigmaMinus(G4bool& inElastic,
00187                                                          const G4double availableEnergy,
00188                                                          G4HEVector pv[],
00189                                                          G4int& vecLen,
00190                                                          const G4HEVector& incidentParticle,
00191                                                          const G4HEVector& targetParticle,
00192                                                          const G4double atomicWeight)
00193 
00194 // AntiSigma- undergoes interaction with nucleon within a nucleus.  Check if it is
00195 // energetically possible to produce pions/kaons.  In not, assume nuclear excitation
00196 // occurs and input particle is degraded in energy. No other particles are produced.
00197 // If reaction is possible, find the correct number of pions/protons/neutrons
00198 // produced using an interpolation to multiplicity data.  Replace some pions or
00199 // protons/neutrons by kaons or strange baryons according to the average
00200 // multiplicity per inelastic reaction.
00201 {
00202   static const G4double expxu = 82;      // upper bound for arg. of exp
00203   static const G4double expxl = -expxu;  // lower bound for arg. of exp
00204 
00205   static const G4double protb = 0.7;
00206   static const G4double neutb = 0.7;
00207   static const G4double     c = 1.25;
00208 
00209   static const G4int numMul   = 1200;
00210   static const G4int numMulAn = 400;
00211   static const G4int numSec   = 60;
00212 
00213   G4int neutronCode = Neutron.getCode();
00214   G4int protonCode  = Proton.getCode();
00215 
00216   G4int targetCode = targetParticle.getCode();
00217   G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00218 
00219   static G4bool first = true;
00220   static G4double protmul[numMul],  protnorm[numSec];   // proton constants
00221   static G4double protmulAn[numMulAn],protnormAn[numSec]; 
00222   static G4double neutmul[numMul],  neutnorm[numSec];   // neutron constants
00223   static G4double neutmulAn[numMulAn],neutnormAn[numSec];
00224 
00225   //  misc. local variables
00226   //  npos = number of pi+,  nneg = number of pi-,  nzero = number of pi0
00227 
00228   G4int i, counter, nt, npos, nneg, nzero;
00229 
00230    if( first ) 
00231      {               // compute normalization constants, this will only be done once
00232        first = false;
00233        for( i=0; i<numMul  ; i++ ) protmul[i]  = 0.0;
00234        for( i=0; i<numSec  ; i++ ) protnorm[i] = 0.0;
00235        counter = -1;
00236        for( npos=0; npos<(numSec/3); npos++ ) 
00237           {
00238             for( nneg=std::max(0,npos-2); nneg<=npos; nneg++ ) 
00239                {
00240                  for( nzero=0; nzero<numSec/3; nzero++ ) 
00241                     {
00242                       if( ++counter < numMul ) 
00243                         {
00244                           nt = npos+nneg+nzero;
00245                           if( (nt>0) && (nt<=numSec) ) 
00246                             {
00247                               protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00248                               protnorm[nt-1] += protmul[counter];
00249                             }
00250                         }
00251                     }
00252                }
00253           }
00254        for( i=0; i<numMul; i++ )neutmul[i]  = 0.0;
00255        for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00256        counter = -1;
00257        for( npos=0; npos<numSec/3; npos++ ) 
00258           {
00259             for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ ) 
00260                {
00261                  for( nzero=0; nzero<numSec/3; nzero++ ) 
00262                     {
00263                       if( ++counter < numMul ) 
00264                         {
00265                           nt = npos+nneg+nzero;
00266                           if( (nt>0) && (nt<=numSec) ) 
00267                             {
00268                                neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00269                                neutnorm[nt-1] += neutmul[counter];
00270                             }
00271                         }
00272                     }
00273                }
00274           }
00275        for( i=0; i<numSec; i++ ) 
00276           {
00277             if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00278             if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00279           }
00280                                                                    // annihilation
00281        for( i=0; i<numMulAn  ; i++ ) protmulAn[i]  = 0.0;
00282        for( i=0; i<numSec    ; i++ ) protnormAn[i] = 0.0;
00283        counter = -1;
00284        for( npos=1; npos<(numSec/3); npos++ ) 
00285           {
00286             nneg = std::max(0,npos-2); 
00287             for( nzero=0; nzero<numSec/3; nzero++ ) 
00288                {
00289                  if( ++counter < numMulAn ) 
00290                    {
00291                      nt = npos+nneg+nzero;
00292                      if( (nt>1) && (nt<=numSec) ) 
00293                        {
00294                          protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00295                          protnormAn[nt-1] += protmulAn[counter];
00296                        }
00297                    }
00298                }
00299           }
00300        for( i=0; i<numMulAn; i++ ) neutmulAn[i]  = 0.0;
00301        for( i=0; i<numSec;   i++ ) neutnormAn[i] = 0.0;
00302        counter = -1;
00303        for( npos=0; npos<numSec/3; npos++ ) 
00304           {
00305             nneg = npos-1; 
00306             for( nzero=0; nzero<numSec/3; nzero++ ) 
00307                {
00308                  if( ++counter < numMulAn ) 
00309                    {
00310                      nt = npos+nneg+nzero;
00311                      if( (nt>1) && (nt<=numSec) ) 
00312                        {
00313                           neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00314                           neutnormAn[nt-1] += neutmulAn[counter];
00315                        }
00316                    }
00317                }
00318           }
00319        for( i=0; i<numSec; i++ ) 
00320           {
00321             if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
00322             if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
00323           }
00324      }                                          // end of initialization
00325 
00326          
00327                                               // initialize the first two places
00328                                               // the same as beam and target                                    
00329    pv[0] = incidentParticle;
00330    pv[1] = targetParticle;
00331    vecLen = 2;
00332 
00333    if( !inElastic ) 
00334      {                                        // some two-body reactions 
00335        G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
00336 
00337        G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5 ));
00338        if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) 
00339          {           
00340            G4double ran = G4UniformRand();
00341 
00342            if ( targetCode == neutronCode)
00343              {
00344                if(ran < 0.2)
00345                  {
00346                    pv[0] = AntiSigmaZero;
00347                    pv[1] = Proton;
00348                  }
00349                else if (ran < 0.4)
00350                  {
00351                    pv[0] = AntiLambda;
00352                    pv[1] = Proton;
00353                  }
00354                else if (ran < 0.6)
00355                  {
00356                    pv[0] = Proton;
00357                    pv[1] = AntiLambda;
00358                  }
00359                else if (ran < 0.8)
00360                  {
00361                    pv[0] = Proton;
00362                    pv[1] = AntiSigmaZero;
00363                  }
00364                else
00365                  {
00366                    pv[0] = Neutron;
00367                    pv[1] = AntiSigmaMinus;
00368                  }     
00369              }
00370            else
00371              {
00372                pv[0] = Proton;
00373                pv[1] = AntiSigmaMinus;
00374              }  
00375          }   
00376        return;
00377      }
00378    else if (availableEnergy <= PionPlus.getMass())
00379        return;
00380 
00381                                                   //   inelastic scattering
00382 
00383   npos = 0; nneg = 0; nzero = 0;
00384   G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88, 
00385                      0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40, 
00386                      0.39, 0.36, 0.33, 0.10, 0.01};
00387   G4int            iplab =      G4int( incidentTotalMomentum*10.);
00388   if ( iplab >  9) iplab = 10 + G4int( (incidentTotalMomentum  -1.)*5. );          
00389   if ( iplab > 14) iplab = 15 + G4int(  incidentTotalMomentum  -2.     );
00390   if ( iplab > 22) iplab = 23 + G4int( (incidentTotalMomentum -10.)/10.); 
00391                    iplab = std::min(24, iplab);
00392 
00393   if (G4UniformRand() > anhl[iplab]) { // non- annihilation channels
00394 
00395     //  number of total particles vs. centre of mass Energy - 2*proton mass
00396     G4double aleab = std::log(availableEnergy);
00397     G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00398                         + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00399    
00400     // normalization constant for kno-distribution.
00401     // calculate first the sum of all constants, check for numerical problems.   
00402     G4double test, dum, anpn = 0.0;
00403 
00404     for (nt=1; nt<=numSec; nt++) {
00405       test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00406       dum = pi*nt/(2.0*n*n);
00407       if (std::fabs(dum) < 1.0) { 
00408         if( test >= 1.0e-10 )anpn += dum*test;
00409       } else { 
00410         anpn += dum*test;
00411       }
00412     }
00413    
00414     G4double ran = G4UniformRand();
00415     G4double excs = 0.0;
00416     if (targetCode == protonCode) {
00417       counter = -1;
00418       for (npos = 0; npos < numSec/3; npos++) {
00419         for (nneg = std::max(0,npos-2); nneg <= npos; nneg++) {
00420           for (nzero = 0; nzero < numSec/3; nzero++) {
00421             if (++counter < numMul) {
00422               nt = npos+nneg+nzero;
00423               if ((nt > 0) && (nt <= numSec) ) {
00424                 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00425                 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00426                 if (std::fabs(dum) < 1.0) { 
00427                   if( test >= 1.0e-10 )excs += dum*test;
00428                 } else { 
00429                   excs += dum*test;
00430                 }
00431 
00432                 if (ran < excs) goto outOfLoop;      //----------------------->
00433               }   
00434             }    
00435           }     
00436         }                                                                                  
00437       } 
00438                              // 3 previous loops continued to the end
00439       inElastic = false;     // quasi-elastic scattering   
00440       return;
00441     } else {   // target must be a neutron
00442       counter = -1;
00443                for( npos=0; npos<numSec/3; npos++ ) 
00444                   {
00445                     for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ ) 
00446                        {
00447                          for( nzero=0; nzero<numSec/3; nzero++ ) 
00448                             {
00449                               if( ++counter < numMul ) 
00450                                 {
00451                                   nt = npos+nneg+nzero;
00452                                   if ( (nt>0) && (nt<=numSec) ) {
00453                                     test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00454                                     dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00455                                     if (std::fabs(dum) < 1.0) { 
00456                                       if( test >= 1.0e-10 )excs += dum*test;
00457                                     } else { 
00458                                       excs += dum*test;
00459                                     }
00460 
00461                                     if (ran < excs) goto outOfLoop;       // -------------------------->
00462                                   }
00463                                 }
00464                             }
00465                        }
00466                   }
00467                                                       // 3 previous loops continued to the end
00468                inElastic = false;                     // quasi-elastic scattering.
00469                return;
00470              }
00471           
00472        outOfLoop:           //  <------------------------------------------------------------------------   
00473     
00474        ran = G4UniformRand();
00475 
00476        if( targetCode == protonCode)
00477          {
00478            if( npos == nneg)
00479              {
00480              }
00481            else if (npos == (nneg+1))
00482              {
00483                if( ran < 0.50)
00484                  {
00485                    pv[1] = Neutron;
00486                  }
00487                else if (ran < 0.75)
00488                  {
00489                    pv[0] = AntiSigmaZero;
00490                  }
00491                else
00492                  {
00493                    pv[0] = AntiLambda;
00494                  }
00495              }
00496            else      
00497              {
00498                if (ran < 0.5)
00499                  {
00500                    pv[0] = AntiSigmaZero;
00501                    pv[1] = Neutron;
00502                  }
00503                else
00504                  {
00505                    pv[0] = AntiLambda;
00506                    pv[1] = Neutron;
00507                  } 
00508              } 
00509          }  
00510        else
00511          {
00512            if( npos == nneg)
00513              {
00514                if (ran < 0.5)
00515                  {
00516                  }
00517                else if(ran < 0.75)
00518                  {
00519                    pv[0] = AntiSigmaZero;
00520                    pv[1] = Proton;
00521                  }
00522                else
00523                  {
00524                    pv[0] = AntiLambda;
00525                    pv[1] = Proton;
00526                  }
00527              } 
00528            else if ( npos == (nneg+1))
00529              {
00530                if (ran < 0.5)
00531                  {
00532                    pv[0] = AntiSigmaZero;
00533                  }
00534                else 
00535                  {
00536                    pv[0] = AntiLambda;
00537                  }
00538              }
00539            else 
00540              {
00541                pv[1] = Proton;
00542              }
00543          }      
00544      
00545      }
00546    else                                                               // annihilation
00547      {  
00548        if ( availableEnergy > 2. * PionPlus.getMass() )
00549          {
00550 
00551            G4double aleab = std::log(availableEnergy);
00552            G4double n     = 3.62567+aleab*(0.665843+aleab*(0.336514
00553                             + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00554    
00555                       //   normalization constant for kno-distribution.
00556                       //   calculate first the sum of all constants, check for numerical problems.   
00557            G4double test, dum, anpn = 0.0;
00558 
00559            for (nt=2; nt<=numSec; nt++) {
00560              test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00561              dum = pi*nt/(2.0*n*n);
00562              if (std::fabs(dum) < 1.0) { 
00563                if( test >= 1.0e-10 )anpn += dum*test;
00564              } else { 
00565                anpn += dum*test;
00566              }
00567            }
00568    
00569            G4double ran = G4UniformRand();
00570            G4double excs = 0.0;
00571            if( targetCode == protonCode ) 
00572              {
00573                counter = -1;
00574                for( npos=2; npos<numSec/3; npos++ ) 
00575                   {
00576                     nneg = npos-2; 
00577                     for( nzero=0; nzero<numSec/3; nzero++ ) 
00578                       {
00579                         if( ++counter < numMulAn ) 
00580                           {
00581                             nt = npos+nneg+nzero;
00582                             if ( (nt>1) && (nt<=numSec) ) {
00583                               test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00584                               dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
00585                               if (std::fabs(dum) < 1.0) { 
00586                                 if( test >= 1.0e-10 )excs += dum*test;
00587                               } else { 
00588                                 excs += dum*test;
00589                               }
00590 
00591                               if (ran < excs) goto outOfLoopAn;      //----------------------->
00592                             }   
00593                           }    
00594                       }     
00595                  }                                                                                  
00596                                               // 3 previous loops continued to the end
00597                inElastic = false;                 // quasi-elastic scattering   
00598                return;
00599              }
00600            else   
00601              {                                         // target must be a neutron
00602                counter = -1;
00603                for( npos=1; npos<numSec/3; npos++ ) 
00604                  { 
00605                    nneg = npos-1; 
00606                    for( nzero=0; nzero<numSec/3; nzero++ ) 
00607                       {
00608                         if (++counter < numMulAn) {
00609                           nt = npos+nneg+nzero;
00610                           if ( (nt>1) && (nt<=numSec) ) {
00611                             test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00612                             dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
00613                             if (std::fabs(dum) < 1.0) { 
00614                               if( test >= 1.0e-10 )excs += dum*test;
00615                             } else { 
00616                               excs += dum*test;
00617                             }
00618 
00619                             if (ran < excs) goto outOfLoopAn;       // -------------------------->
00620                           }
00621                         }
00622                       }
00623                  }
00624                inElastic = false;                     // quasi-elastic scattering.
00625                return;
00626              }
00627        outOfLoopAn:           //  <------------------------------------------------------------------   
00628        vecLen = 0;
00629          }
00630      }
00631 
00632    nt = npos + nneg + nzero;
00633    while ( nt > 0)
00634        {
00635          G4double ran = G4UniformRand();
00636          if ( ran < (G4double)npos/nt)
00637             { 
00638               if( npos > 0 ) 
00639                 { pv[vecLen++] = PionPlus;
00640                   npos--;
00641                 }
00642             }
00643          else if ( ran < (G4double)(npos+nneg)/nt)
00644             {   
00645               if( nneg > 0 )
00646                 { 
00647                   pv[vecLen++] = PionMinus;
00648                   nneg--;
00649                 }
00650             }
00651          else
00652             {
00653               if( nzero > 0 )
00654                 { 
00655                   pv[vecLen++] = PionZero;
00656                   nzero--;
00657                 }
00658             }
00659          nt = npos + nneg + nzero;
00660        } 
00661    if (verboseLevel > 1)
00662       {
00663         G4cout << "Particles produced: " ;
00664         G4cout << pv[0].getName() << " " ;
00665         G4cout << pv[1].getName() << " " ;
00666         for (i=2; i < vecLen; i++)   
00667             { 
00668               G4cout << pv[i].getName() << " " ;
00669             }
00670          G4cout << G4endl;
00671       }
00672    return;
00673  }
00674 
00675 
00676 
00677 
00678 
00679 
00680 
00681 
00682 
00683 

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