G4HEAntiNeutronInelastic Class Reference

#include <G4HEAntiNeutronInelastic.hh>

Inheritance diagram for G4HEAntiNeutronInelastic:

G4HEInelastic G4HadronicInteraction

Public Member Functions

 G4HEAntiNeutronInelastic (const G4String &name="G4HEAntiNeutronInelastic")
 ~G4HEAntiNeutronInelastic ()
virtual void ModelDescription (std::ostream &) const
G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
G4int GetNumberOfSecondaries ()
void FirstIntInCasAntiNeutron (G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, const G4double atomicWeight)

Data Fields

G4int vecLength

Detailed Description

Definition at line 50 of file G4HEAntiNeutronInelastic.hh.


Constructor & Destructor Documentation

G4HEAntiNeutronInelastic::G4HEAntiNeutronInelastic ( const G4String name = "G4HEAntiNeutronInelastic"  ) 

Definition at line 43 of file G4HEAntiNeutronInelastic.cc.

References G4cout, G4endl, G4HEInelastic::MAXPART, G4HadronicInteraction::theMaxEnergy, G4HadronicInteraction::theMinEnergy, vecLength, and G4HEInelastic::verboseLevel.

00044  : G4HEInelastic(name)
00045 {
00046   vecLength = 0;
00047   theMinEnergy = 20*GeV;
00048   theMaxEnergy = 10*TeV;
00049   MAXPART      = 2048;
00050   verboseLevel = 0;
00051   G4cout << "WARNING: model G4HEAntiNeutronInelastic is being deprecated and will\n"
00052          << "disappear in Geant4 version 10.0"  << G4endl;  
00053 }

G4HEAntiNeutronInelastic::~G4HEAntiNeutronInelastic (  )  [inline]

Definition at line 54 of file G4HEAntiNeutronInelastic.hh.

00054 {};


Member Function Documentation

G4HadFinalState * G4HEAntiNeutronInelastic::ApplyYourself ( const G4HadProjectile aTrack,
G4Nucleus targetNucleus 
) [virtual]

Implements G4HadronicInteraction.

Definition at line 71 of file G4HEAntiNeutronInelastic.cc.

References G4HEInelastic::ElasticScattering(), G4HEInelastic::FillParticleChange(), FirstIntInCasAntiNeutron(), G4cout, G4endl, G4UniformRand, G4Nucleus::GetA_asInt(), G4HEVector::getCode(), G4HEVector::getEnergy(), G4HEVector::getMass(), G4HEVector::getName(), G4Nucleus::GetZ_asInt(), G4HEInelastic::HighEnergyCascading(), G4HEInelastic::HighEnergyClusterProduction(), G4HEInelastic::MAXPART, G4HEInelastic::MediumEnergyCascading(), G4HEInelastic::MediumEnergyClusterProduction(), G4HEInelastic::NuclearExcitation(), G4HEInelastic::NuclearInelasticity(), G4HEInelastic::QuasiElasticScattering(), G4HEVector::setDefinition(), G4HadFinalState::SetStatusChange(), stopAndKill, G4HEInelastic::StrangeParticlePairProduction(), G4HadronicInteraction::theParticleChange, vecLength, and G4HEInelastic::verboseLevel.

00073 {
00074   G4HEVector* pv = new G4HEVector[MAXPART];
00075   const G4HadProjectile* aParticle = &aTrack;
00076   const G4double atomicWeight = targetNucleus.GetA_asInt();
00077   const G4double atomicNumber = targetNucleus.GetZ_asInt();
00078   G4HEVector incidentParticle(aParticle);
00079 
00080   G4int incidentCode = incidentParticle.getCode();
00081   G4double incidentMass = incidentParticle.getMass();
00082   G4double incidentTotalEnergy = incidentParticle.getEnergy();
00083 
00084   // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00085   // DHW 19 May 2011: variable set but not used
00086 
00087   G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
00088 
00089   if (incidentKineticEnergy < 1.)
00090     G4cout << "GHEAntiNeutronInelastic: incident energy < 1 GeV" << G4endl;;
00091 
00092   if (verboseLevel > 1) {
00093     G4cout << "G4HEAntiNeutronInelastic::ApplyYourself" << G4endl;
00094     G4cout << "incident particle " << incidentParticle.getName()
00095            << "mass "              << incidentMass
00096            << "kinetic energy "    << incidentKineticEnergy
00097            << G4endl;
00098     G4cout << "target material with (A,Z) = (" 
00099            << atomicWeight << "," << atomicNumber << ")" << G4endl;
00100   }
00101 
00102   G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 
00103                                               atomicWeight, atomicNumber);
00104   if (verboseLevel > 1)
00105     G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
00106     
00107   incidentKineticEnergy -= inelasticity;
00108     
00109   G4double excitationEnergyGNP = 0.;
00110   G4double excitationEnergyDTA = 0.; 
00111 
00112   G4double excitation = NuclearExcitation(incidentKineticEnergy,
00113                                           atomicWeight, atomicNumber,
00114                                           excitationEnergyGNP,
00115                                           excitationEnergyDTA);
00116   if (verboseLevel > 1)
00117     G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 
00118            << excitationEnergyDTA << G4endl;             
00119 
00120   incidentKineticEnergy -= excitation;
00121   incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00122   // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
00123   //                                   *(incidentTotalEnergy+incidentMass));
00124   //  DHW 19 May 2011: variable set but not used
00125 
00126   G4HEVector targetParticle;
00127   if (G4UniformRand() < atomicNumber/atomicWeight) { 
00128     targetParticle.setDefinition("Proton");
00129   } else { 
00130     targetParticle.setDefinition("Neutron");
00131   }
00132 
00133   G4double targetMass = targetParticle.getMass();
00134   G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
00135                                         + targetMass*targetMass
00136                                         + 2.0*targetMass*incidentTotalEnergy);
00137   G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00138 
00139   G4bool inElastic = true;
00140   vecLength = 0;           
00141         
00142   if (verboseLevel > 1)
00143     G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00144            << incidentCode << G4endl;
00145 
00146   G4bool successful = false; 
00147     
00148   FirstIntInCasAntiNeutron(inElastic, availableEnergy, pv, vecLength,
00149                            incidentParticle, targetParticle, atomicWeight);
00150 
00151   if (verboseLevel > 1)
00152     G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00153 
00154   if ((vecLength > 0) && (availableEnergy > 1.)) 
00155     StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00156                                   pv, vecLength,
00157                                   incidentParticle, targetParticle);
00158   HighEnergyCascading(successful, pv, vecLength,
00159                       excitationEnergyGNP, excitationEnergyDTA,
00160                       incidentParticle, targetParticle,
00161                       atomicWeight, atomicNumber);
00162   if (!successful)
00163     HighEnergyClusterProduction(successful, pv, vecLength,
00164                                 excitationEnergyGNP, excitationEnergyDTA,
00165                                 incidentParticle, targetParticle,
00166                                 atomicWeight, atomicNumber);
00167   if (!successful) 
00168     MediumEnergyCascading(successful, pv, vecLength, 
00169                           excitationEnergyGNP, excitationEnergyDTA, 
00170                           incidentParticle, targetParticle,
00171                           atomicWeight, atomicNumber);
00172 
00173   if (!successful)
00174     MediumEnergyClusterProduction(successful, pv, vecLength,
00175                                   excitationEnergyGNP, excitationEnergyDTA,       
00176                                   incidentParticle, targetParticle,
00177                                   atomicWeight, atomicNumber);
00178   if (!successful)
00179     QuasiElasticScattering(successful, pv, vecLength,
00180                            excitationEnergyGNP, excitationEnergyDTA,
00181                            incidentParticle, targetParticle, 
00182                            atomicWeight, atomicNumber);
00183   if (!successful)
00184     ElasticScattering(successful, pv, vecLength,
00185                       incidentParticle,    
00186                       atomicWeight, atomicNumber);
00187 
00188   if (!successful)
00189     G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 
00190            << G4endl;
00191 
00192   FillParticleChange(pv,  vecLength);
00193   delete [] pv;
00194   theParticleChange.SetStatusChange(stopAndKill);
00195   return &theParticleChange;
00196 }

void G4HEAntiNeutronInelastic::FirstIntInCasAntiNeutron ( G4bool inElastic,
const G4double  availableEnergy,
G4HEVector  pv[],
G4int vecLen,
const G4HEVector incidentParticle,
const G4HEVector targetParticle,
const G4double  atomicWeight 
)

Definition at line 200 of file G4HEAntiNeutronInelastic.cc.

References G4HEInelastic::AntiProton, G4cout, G4endl, G4UniformRand, G4HEVector::getCode(), G4HEVector::getMass(), G4HEVector::getName(), G4HEVector::getTotalMomentum(), CLHEP::detail::n, G4HEInelastic::Neutron, neutronCode, G4INCL::Math::pi, G4HEInelastic::PionMinus, G4HEInelastic::PionPlus, G4HEInelastic::PionZero, G4HEInelastic::pmltpc(), G4HEInelastic::Proton, sqr(), and G4HEInelastic::verboseLevel.

Referenced by ApplyYourself().

00216 {
00217   static const G4double expxu = 82.;     // upper bound for arg. of exp
00218   static const G4double expxl = -expxu;  // lower bound for arg. of exp
00219 
00220   static const G4double protb = 0.7;
00221   static const G4double neutb = 0.7;
00222   static const G4double     c = 1.25;
00223 
00224   static const G4int numMul   = 1200;
00225   static const G4int numMulAn = 400;
00226   static const G4int numSec   = 60;
00227 
00228   G4int neutronCode = Neutron.getCode();
00229   G4int protonCode = Proton.getCode();
00230 
00231   G4int targetCode = targetParticle.getCode();
00232   G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00233 
00234   static G4bool first = true;
00235   static G4double protmul[numMul],  protnorm[numSec];   // proton constants
00236   static G4double protmulAn[numMulAn],protnormAn[numSec]; 
00237   static G4double neutmul[numMul],  neutnorm[numSec];   // neutron constants
00238   static G4double neutmulAn[numMulAn],neutnormAn[numSec];
00239 
00240   //  misc. local variables
00241   //  npos = number of pi+,  nneg = number of pi-,  nzero = number of pi0
00242 
00243   G4int i, counter, nt, npos, nneg, nzero;
00244 
00245    if( first ) 
00246      {                 // compute normalization constants, this will only be done once
00247        first = false;
00248        for( i=0; i<numMul  ; i++ ) protmul[i]  = 0.0;
00249        for( i=0; i<numSec  ; i++ ) protnorm[i] = 0.0;
00250        counter = -1;
00251        for( npos=0; npos<(numSec/3); npos++ ) 
00252           {
00253             for( nneg=std::max(0,npos-2); nneg<=npos; nneg++ ) 
00254                {
00255                  for( nzero=0; nzero<numSec/3; nzero++ ) 
00256                     {
00257                       if( ++counter < numMul ) 
00258                         {
00259                           nt = npos+nneg+nzero;
00260                           if( (nt>0) && (nt<=numSec) ) 
00261                             {
00262                               protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00263                               protnorm[nt-1] += protmul[counter];
00264                             }
00265                         }
00266                     }
00267                }
00268           }
00269        for( i=0; i<numMul; i++ )neutmul[i]  = 0.0;
00270        for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00271        counter = -1;
00272        for( npos=0; npos<numSec/3; npos++ ) 
00273           {
00274             for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ ) 
00275                {
00276                  for( nzero=0; nzero<numSec/3; nzero++ ) 
00277                     {
00278                       if( ++counter < numMul ) 
00279                         {
00280                           nt = npos+nneg+nzero;
00281                           if( (nt>0) && (nt<=numSec) ) 
00282                             {
00283                                neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00284                                neutnorm[nt-1] += neutmul[counter];
00285                             }
00286                         }
00287                     }
00288                }
00289           }
00290        for( i=0; i<numSec; i++ ) 
00291           {
00292             if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00293             if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00294           }
00295                                                                    // annihilation
00296        for( i=0; i<numMulAn  ; i++ ) protmulAn[i]  = 0.0;
00297        for( i=0; i<numSec    ; i++ ) protnormAn[i] = 0.0;
00298        counter = -1;
00299        for( npos=1; npos<(numSec/3); npos++ ) 
00300           {
00301             nneg = std::max(0,npos-1); 
00302             for( nzero=0; nzero<numSec/3; nzero++ ) 
00303                {
00304                  if( ++counter < numMulAn ) 
00305                    {
00306                      nt = npos+nneg+nzero;
00307                      if( (nt>1) && (nt<=numSec) ) 
00308                        {
00309                          protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00310                          protnormAn[nt-1] += protmulAn[counter];
00311                        }
00312                    }
00313                }
00314           }
00315        for( i=0; i<numMulAn; i++ ) neutmulAn[i]  = 0.0;
00316        for( i=0; i<numSec;   i++ ) neutnormAn[i] = 0.0;
00317        counter = -1;
00318        for( npos=0; npos<numSec/3; npos++ ) 
00319           {
00320             nneg = npos; 
00321             for( nzero=0; nzero<numSec/3; nzero++ ) 
00322                {
00323                  if( ++counter < numMulAn ) 
00324                    {
00325                      nt = npos+nneg+nzero;
00326                      if( (nt>1) && (nt<=numSec) ) 
00327                        {
00328                           neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00329                           neutnormAn[nt-1] += neutmulAn[counter];
00330                        }
00331                    }
00332                }
00333           }
00334        for( i=0; i<numSec; i++ ) 
00335           {
00336             if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
00337             if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
00338           }
00339      }                                          // end of initialization
00340 
00341          
00342                                               // initialize the first two places
00343                                               // the same as beam and target                                    
00344    pv[0] = incidentParticle;
00345    pv[1] = targetParticle;
00346    vecLen = 2;
00347 
00348    if( !inElastic ) 
00349      {                                        // nb n --> pb p
00350        if( targetCode == neutronCode ) 
00351          {
00352            G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
00353 
00354            G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5));
00355            if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) 
00356              {
00357                pv[0] = AntiProton;
00358                pv[1] = Proton;
00359              }
00360          }
00361        return;
00362      }
00363    else if (availableEnergy <= PionPlus.getMass())
00364        return;
00365 
00366                                                   //   inelastic scattering
00367 
00368    npos = 0, nneg = 0, nzero = 0;
00369    G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88, 
00370                       0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40, 
00371                       0.39, 0.36, 0.33, 0.10, 0.01};
00372    G4int            iplab =      G4int( incidentTotalMomentum*10.);
00373    if ( iplab >  9) iplab = 10 + G4int( (incidentTotalMomentum  -1.)*5. );          
00374    if ( iplab > 14) iplab = 15 + G4int(  incidentTotalMomentum  -2.     );
00375    if ( iplab > 22) iplab = 23 + G4int( (incidentTotalMomentum -10.)/10.); 
00376                     iplab = std::min(24, iplab);
00377 
00378    if ( G4UniformRand() > anhl[iplab] )
00379      {
00380 
00381        G4double eab = availableEnergy;
00382        G4int ieab = G4int( eab*5.0 );
00383    
00384        G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
00385        if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) ) 
00386          {
00387                                               //   suppress high multiplicity events at low momentum
00388                                               //   only one additional pion will be produced
00389            G4double w0, wp, wm, wt, ran;
00390            if( targetCode == protonCode )                    // target is a proton 
00391              {
00392                w0 = - sqr(1.+protb)/(2.*c*c);
00393                w0 = wp = std::exp(w0);
00394                if( G4UniformRand() < w0/(w0+wp) )  
00395                  { npos = 0; nneg = 0; nzero = 1; }
00396                else 
00397                  { npos = 1; nneg = 0; nzero = 0; }       
00398              } 
00399            else 
00400              {                                               // target is a neutron
00401                w0 = -sqr(1.+neutb)/(2.*c*c);
00402                w0 = wp = std::exp(w0);
00403                wm = -sqr(-1.+neutb)/(2.*c*c);
00404                wm = std::exp(wm);
00405                wt = w0+wp+wm;
00406                wp = w0+wp;
00407                ran = G4UniformRand();
00408                if( ran < w0/wt)
00409                  { npos = 0; nneg = 0; nzero = 1; }       
00410                else if( ran < wp/wt)
00411                  { npos = 1; nneg = 0; nzero = 0; }       
00412                else
00413                  { npos = 0; nneg = 1; nzero = 0; }       
00414              }
00415          }
00416        else
00417          {
00418                          //  number of total particles vs. centre of mass Energy - 2*proton mass
00419    
00420            G4double aleab = std::log(availableEnergy);
00421            G4double n     = 3.62567+aleab*(0.665843+aleab*(0.336514
00422                             + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00423    
00424                          // normalization constant for kno-distribution.
00425                          // calculate first the sum of all constants, check for numerical problems.   
00426            G4double test, dum, anpn = 0.0;
00427 
00428            for (nt=1; nt<=numSec; nt++) {
00429                test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00430                dum = pi*nt/(2.0*n*n);
00431 
00432                if (std::fabs(dum) < 1.0) { 
00433                  if( test >= 1.0e-10 )anpn += dum*test;
00434                } else { 
00435                  anpn += dum*test;
00436                }
00437            }
00438    
00439            G4double ran = G4UniformRand();
00440            G4double excs = 0.0;
00441            if (targetCode == protonCode) {
00442              counter = -1;
00443              for (npos=0; npos<numSec/3; npos++) {
00444                for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
00445                  for (nzero=0; nzero<numSec/3; nzero++) {
00446                    if (++counter < numMul) {
00447                      nt = npos+nneg+nzero;
00448                      if ( (nt>0) && (nt<=numSec) ) {
00449                        test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00450                        dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00451 
00452                        if (std::fabs(dum) < 1.0) { 
00453                          if( test >= 1.0e-10 )excs += dum*test;
00454                        } else { 
00455                          excs += dum*test;
00456                        }
00457 
00458                        if (ran < excs) goto outOfLoop;      //----------------------->
00459                      }  
00460                    }    
00461                  }     
00462                }                                                                                  
00463              }
00464        
00465                                            // 3 previous loops continued to the end
00466              inElastic = false;            // quasi-elastic scattering   
00467              return;
00468 
00469            } else {                        // target must be a neutron
00470              counter = -1;
00471              for (npos=0; npos<numSec/3; npos++) {
00472                for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
00473                  for (nzero=0; nzero<numSec/3; nzero++) {
00474                    if (++counter < numMul) {
00475                      nt = npos+nneg+nzero;
00476                      if ((nt>0) && (nt<=numSec) ) {
00477                        test = std::exp( std::min( expxu, std::max( 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 == protonCode)
00499          {
00500            if( npos == nneg)
00501              {
00502              }
00503            else if (npos == (nneg+1))
00504              {
00505                if( G4UniformRand() < 0.5)
00506                  {
00507                    pv[1] = Neutron;
00508                  }
00509                else
00510                  {
00511                    pv[0] = AntiProton;
00512                  }
00513              }
00514            else      
00515              {
00516                pv[0] = AntiProton;
00517                pv[1] = Neutron;
00518              } 
00519          }  
00520        else
00521          {
00522            if( npos == nneg)
00523              {
00524                if( G4UniformRand() < 0.25)
00525                  {
00526                    pv[0] = AntiProton;
00527                    pv[1] = Proton;
00528                  }
00529                else
00530                  {
00531                  }
00532              } 
00533            else if ( npos == (nneg-1))
00534              {
00535                pv[1] = Proton;
00536              }
00537            else
00538              {
00539                pv[0] = AntiProton;
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( std::min( expxu, std::max( 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              counter = -1;
00571              for (npos=1; npos<numSec/3; npos++) {
00572                nneg = npos-1; 
00573                for (nzero=0; nzero<numSec/3; nzero++) {
00574                  if (++counter < numMulAn) {
00575                    nt = npos+nneg+nzero;
00576                    if ( (nt>1) && (nt<=numSec) ) {
00577                      test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00578                      dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
00579 
00580                      if (std::fabs(dum) < 1.0) { 
00581                        if( test >= 1.0e-10 )excs += dum*test;
00582                      } else { 
00583                        excs += dum*test;
00584                      }
00585 
00586                      if (ran < excs) goto outOfLoopAn;      //----------------------->
00587                    }
00588                  }
00589                }
00590              }
00591                                              // 3 previous loops continued to the end
00592              inElastic = false;            // quasi-elastic scattering   
00593              return;
00594 
00595            } else {                          // target must be a neutron
00596              counter = -1;
00597              for (npos=0; npos<numSec/3; npos++) { 
00598                nneg = npos; 
00599                for (nzero=0; nzero<numSec/3; nzero++) {
00600                  if (++counter < numMulAn) {
00601                    nt = npos+nneg+nzero;
00602                    if ( (nt>1) && (nt<=numSec) ) {
00603                      test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00604                      dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
00605 
00606                      if (std::fabs(dum) < 1.0) { 
00607                        if( test >= 1.0e-10 )excs += dum*test;
00608                      } else { 
00609                        excs += dum*test;
00610                      }
00611 
00612                      if (ran < excs) goto outOfLoopAn;       // -------------------------->
00613                    }
00614                  }
00615                }
00616              }
00617              inElastic = false;                     // quasi-elastic scattering.
00618              return;
00619            }
00620        outOfLoopAn:           //  <------------------------------------------------------------------   
00621        vecLen = 0;
00622          }
00623      }
00624 
00625    nt = npos + nneg + nzero;
00626    while ( nt > 0)
00627        {
00628          G4double ran = G4UniformRand();
00629          if ( ran < (G4double)npos/nt)
00630             { 
00631               if( npos > 0 ) 
00632                 { pv[vecLen++] = PionPlus;
00633                   npos--;
00634                 }
00635             }
00636          else if ( ran < (G4double)(npos+nneg)/nt)
00637             {   
00638               if( nneg > 0 )
00639                 { 
00640                   pv[vecLen++] = PionMinus;
00641                   nneg--;
00642                 }
00643             }
00644          else
00645             {
00646               if( nzero > 0 )
00647                 { 
00648                   pv[vecLen++] = PionZero;
00649                   nzero--;
00650                 }
00651             }
00652          nt = npos + nneg + nzero;
00653        } 
00654    if (verboseLevel > 1)
00655       {
00656         G4cout << "Particles produced: " ;
00657         G4cout << pv[0].getName() << " " ;
00658         G4cout << pv[1].getName() << " " ;
00659         for (i=2; i < vecLen; i++)   
00660             { 
00661               G4cout << pv[i].getName() << " " ;
00662             }
00663          G4cout << G4endl;
00664       }
00665    return;
00666  }

G4int G4HEAntiNeutronInelastic::GetNumberOfSecondaries (  )  [inline]

Definition at line 63 of file G4HEAntiNeutronInelastic.hh.

References vecLength.

00063 {return vecLength;}         

void G4HEAntiNeutronInelastic::ModelDescription ( std::ostream &   )  const [virtual]

Reimplemented from G4HadronicInteraction.

Definition at line 56 of file G4HEAntiNeutronInelastic.cc.

00057 {
00058   outFile << "G4HEAntiNeutronInelastic is one of the High Energy\n"
00059           << "Parameterized (HEP) models used to implement inelastic\n"
00060           << "anti-neutron scattering from nuclei.  It is a re-engineered\n"
00061           << "version of the GHEISHA code of H. Fesefeldt.  It divides the\n"
00062           << "initial collision products into backward- and forward-going\n"
00063           << "clusters which are then decayed into final state hadrons.\n"
00064           << "The model does not conserve energy on an event-by-event\n"
00065           << "basis.  It may be applied to anti-neutrons with initial energies\n"
00066           << "above 20 GeV.\n";
00067 }


Field Documentation

G4int G4HEAntiNeutronInelastic::vecLength

Definition at line 58 of file G4HEAntiNeutronInelastic.hh.

Referenced by ApplyYourself(), G4HEAntiNeutronInelastic(), and GetNumberOfSecondaries().


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:10 2013 for Geant4 by  doxygen 1.4.7