Geant4: G4HENeutronInelastic Class Reference

G4HENeutronInelastic::G4HENeutronInelastic ( ) [inline]

Definition at line 55 of file G4HENeutronInelastic.hh.

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

00055                            : G4HEInelastic("G4HENeutronInelastic")
00056     {
00057       vecLength = 0;
00058       theMinEnergy = 45*CLHEP::GeV;
00059       theMaxEnergy = 10*CLHEP::TeV;
00060       MAXPART      = 2048;
00061       verboseLevel = 0; 
00062       G4cout << "WARNING: model G4HENeutronInelastic is being deprecated and will\n"
00063              << "disappear in Geant4 version 10.0"  << G4endl;  
00064     }

G4HENeutronInelastic::~G4HENeutronInelastic ( ) [inline]

Definition at line 66 of file G4HENeutronInelastic.hh.

00066 {};

G4HadFinalState * G4HENeutronInelastic::ApplyYourself	(	const G4HadProjectile &	aTrack,
		G4Nucleus &	targetNucleus
	)			`[virtual]`

Implements G4HadronicInteraction.

Definition at line 57 of file G4HENeutronInelastic.cc.

References G4HEInelastic::ElasticScattering(), G4HEInelastic::FillParticleChange(), FirstIntInCasNeutron(), 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.

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

void G4HENeutronInelastic::FirstIntInCasNeutron	(	G4bool &	inElastic,
		const G4double	availableEnergy,
		G4HEVector	pv[],
		G4int &	vecLen,
		const G4HEVector &	incidentParticle,
		const G4HEVector &	targetParticle,
		const G4double	atomicWeight
	)

Definition at line 197 of file G4HENeutronInelastic.cc.

References G4HEInelastic::Factorial(), 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().

00212 {
00213   static const G4double expxu = 82.;     // upper bound for arg. of exp
00214   static const G4double expxl = -expxu;  // lower bound for arg. of exp
00215 
00216   static const G4double protb = 0.35;
00217   static const G4double neutb = 0.35;              
00218   static const G4double c = 1.25;
00219 
00220   static const G4int numMul = 1200;
00221   static const G4int numSec = 60;
00222 
00223   G4int neutronCode = Neutron.getCode();
00224   G4int protonCode  = Proton.getCode();
00225 
00226   G4int targetCode = targetParticle.getCode();
00227   G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00228 
00229   static G4bool first = true;
00230   static G4double protmul[numMul], protnorm[numSec];  // proton constants
00231   static G4double neutmul[numMul], neutnorm[numSec];  // neutron constants
00232 
00233   //  misc. local variables
00234   //  npos = number of pi+,  nneg = number of pi-,  nzero = number of pi0
00235 
00236   G4int i, counter, nt, npos, nneg, nzero;
00237 
00238   if (first) {
00239     // compute normalization constants, this will only be done once
00240     first = false;
00241     for (i = 0; i < numMul; i++) protmul[i] = 0.0;
00242     for (i = 0; i < numSec; i++) protnorm[i] = 0.0;
00243     counter = -1;
00244     for (npos = 0; npos < (numSec/3); npos++) {
00245       for (nneg = std::max(0,npos-1); nneg <= (npos+1); nneg++) {
00246                  for( nzero=0; nzero<numSec/3; nzero++ ) 
00247                     {
00248                       if( ++counter < numMul ) 
00249                         {
00250                           nt = npos+nneg+nzero;
00251                           if( (nt>0) && (nt<=numSec) ) 
00252                             {
00253                               protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c)
00254                                                  /(Factorial(1-npos+nneg)*Factorial(1+npos-nneg)) ;
00255                               protnorm[nt-1] += protmul[counter];
00256                             }
00257                         }
00258                     }
00259       }
00260     }
00261 
00262     for( i=0; i<numMul; i++ )neutmul[i]  = 0.0;
00263     for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00264     counter = -1;
00265        for( npos=0; npos<numSec/3; npos++ ) 
00266           {
00267             for( nneg=npos; nneg<=(npos+2); nneg++ ) 
00268                {
00269                  for( nzero=0; nzero<numSec/3; nzero++ ) 
00270                     {
00271                       if( ++counter < numMul ) 
00272                         {
00273                           nt = npos+nneg+nzero;
00274                           if( (nt>0) && (nt<=numSec) ) 
00275                             {
00276                                neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c)
00277                                                   /(Factorial(-npos+nneg)*Factorial(2+npos-nneg));
00278                                neutnorm[nt-1] += neutmul[counter];
00279                             }
00280                         }
00281                     }
00282                }
00283           }
00284        for( i=0; i<numSec; i++ ) 
00285           {
00286             if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00287             if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00288           }
00289      }                                          // end of initialization
00290 
00291          
00292                                               // initialize the first two places
00293                                               // the same as beam and target                                    
00294    pv[0] = incidentParticle;
00295    pv[1] = targetParticle;
00296    vecLen = 2;
00297 
00298    if( !inElastic ) 
00299      {                                     // quasi-elastic scattering, no pions produced
00300        if( targetCode == protonCode ) 
00301          {
00302            G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
00303            G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*2.5) );
00304            if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) 
00305              {                                            // charge exchange  n p -> p n
00306                pv[0] = Proton;
00307                pv[1] = Neutron;
00308              }
00309          }
00310        return;
00311      }
00312    else if (availableEnergy <= PionPlus.getMass())
00313        return;
00314 
00315                                                   //   inelastic scattering
00316 
00317    npos = 0, nneg = 0, nzero = 0;
00318    G4double eab = availableEnergy;
00319    G4int ieab = G4int( eab*5.0 );
00320    
00321    G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
00322    if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) ) 
00323      {
00324                                           //   suppress high multiplicity events at low momentum
00325                                           //   only one additional pion will be produced
00326        G4double w0, wp, wm, wt, ran;
00327        if( targetCode == neutronCode )                    // target is a neutron 
00328          {
00329            w0 = - sqr(1.+neutb)/(2.*c*c);
00330            wm = w0 = std::exp(w0);
00331            w0 = w0/2.;
00332            if( G4UniformRand() < w0/(w0+wm) )  
00333              { npos = 0; nneg = 0; nzero = 1; }
00334            else 
00335              { npos = 0; nneg = 1; nzero = 0; }       
00336          } 
00337        else 
00338          {                                               // target is a proton
00339            w0 = -sqr(1.+protb)/(2.*c*c);
00340            w0 = std::exp(w0);
00341            wp = w0/2.;
00342            wm = -sqr(-1.+protb)/(2.*c*c);
00343            wm = std::exp(wm)/2.;
00344            wt = w0+wp+wm;
00345            wp = w0+wp;
00346            ran = G4UniformRand();
00347            if( ran < w0/wt)
00348              { npos = 0; nneg = 0; nzero = 1; }       
00349            else if( ran < wp/wt)
00350              { npos = 1; nneg = 0; nzero = 0; }       
00351            else
00352              { npos = 0; nneg = 1; nzero = 0; }       
00353          }
00354      }
00355    else
00356      {
00357        // number of total particles vs. centre of mass Energy - 2*proton mass
00358    
00359        G4double aleab = std::log(availableEnergy);
00360        G4double n     = 3.62567+aleab*(0.665843+aleab*(0.336514
00361                     + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00362    
00363        // normalization constant for kno-distribution.
00364        // calculate first the sum of all constants, check for numerical problems.   
00365        G4double test, dum, anpn = 0.0;
00366 
00367        for (nt=1; nt<=numSec; nt++) {
00368          test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00369          dum = pi*nt/(2.0*n*n);
00370          if (std::fabs(dum) < 1.0) { 
00371            if( test >= 1.0e-10 )anpn += dum*test;
00372          } else { 
00373            anpn += dum*test;
00374          }
00375        }
00376    
00377        G4double ran = G4UniformRand();
00378        G4double excs = 0.0;
00379        if( targetCode == protonCode ) 
00380          {
00381            counter = -1;
00382            for (npos=0; npos<numSec/3; npos++) {
00383              for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
00384                for (nzero=0; nzero<numSec/3; nzero++) {
00385                  if (++counter < numMul) {
00386                    nt = npos+nneg+nzero;
00387                    if ( (nt>0) && (nt<=numSec) ) {
00388                      test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00389                      dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00390                      if (std::fabs(dum) < 1.0) { 
00391                        if( test >= 1.0e-10 )excs += dum*test;
00392                      } else { 
00393                        excs += dum*test;
00394                      }
00395                      if (ran < excs) goto outOfLoop;      //----------------------->
00396                    }
00397                  }
00398                }
00399              }
00400            }
00401        
00402            // 3 previous loops continued to the end
00403            inElastic = false;                 // quasi-elastic scattering   
00404            return;
00405          }
00406        else   
00407          {                                         // target must be a neutron
00408            counter = -1;
00409            for (npos=0; npos<numSec/3; npos++) {
00410              for (nneg=npos; nneg<=(npos+2); nneg++) {
00411                for (nzero=0; nzero<numSec/3; nzero++) {
00412                  if (++counter < numMul) {
00413                    nt = npos+nneg+nzero;
00414                    if ( (nt>=1) && (nt<=numSec) ) {
00415                      test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00416                      dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00417                      if (std::fabs(dum) < 1.0) { 
00418                        if( test >= 1.0e-10 )excs += dum*test;
00419                      } else { 
00420                        excs += dum*test;
00421                      }
00422                      if (ran < excs) goto outOfLoop;       // -------------------------->
00423                    }
00424                  }
00425                }
00426              }
00427            }
00428            // 3 previous loops continued to the end
00429            inElastic = false;                     // quasi-elastic scattering.
00430            return;
00431          }
00432      } 
00433    outOfLoop:           //  <----------------------------------------------   
00434     
00435    if( targetCode == neutronCode)
00436      {
00437        if( npos == nneg)
00438          {
00439          }
00440        else if (npos == (nneg-1))
00441          {
00442            if( G4UniformRand() < 0.5)
00443              {
00444                pv[1] = Proton;
00445              }
00446            else
00447              {
00448                pv[0] = Proton;
00449              }
00450          }
00451        else      
00452          {
00453            pv[0] = Proton;
00454            pv[1] = Proton;
00455          } 
00456      }  
00457    else
00458      {
00459        if( npos == nneg)
00460          {
00461            if( G4UniformRand() < 0.25)
00462              {
00463                pv[0] = Proton;
00464                pv[1] = Neutron;
00465              }
00466            else
00467              {
00468              }
00469          } 
00470        else if ( npos == (1+nneg))
00471          {
00472            pv[1] = Neutron;
00473          }
00474        else
00475          {
00476            pv[0] = Proton;
00477          }
00478      }      
00479 
00480 
00481    nt = npos + nneg + nzero;
00482    while ( nt > 0)
00483        {
00484          G4double ran = G4UniformRand();
00485          if ( ran < (G4double)npos/nt)
00486             { 
00487               if( npos > 0 ) 
00488                 { pv[vecLen++] = PionPlus;
00489                   npos--;
00490                 }
00491             }
00492          else if ( ran < (G4double)(npos+nneg)/nt)
00493             {   
00494               if( nneg > 0 )
00495                 { 
00496                   pv[vecLen++] = PionMinus;
00497                   nneg--;
00498                 }
00499             }
00500          else
00501             {
00502               if( nzero > 0 )
00503                 { 
00504                   pv[vecLen++] = PionZero;
00505                   nzero--;
00506                 }
00507             }
00508          nt = npos + nneg + nzero;
00509        } 
00510    if (verboseLevel > 1)
00511       {
00512         G4cout << "Particles produced: " ;
00513         G4cout << pv[0].getName() << " " ;
00514         G4cout << pv[1].getName() << " " ;
00515         for (i=2; i < vecLen; i++)   
00516             { 
00517               G4cout << pv[i].getName() << " " ;
00518             }
00519          G4cout << G4endl;
00520       }
00521    return;
00522  }

G4int G4HENeutronInelastic::GetNumberOfSecondaries ( ) [inline]

Definition at line 75 of file G4HENeutronInelastic.hh.

References vecLength.

00075 {return vecLength;}

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

Reimplemented from G4HadronicInteraction.

Definition at line 42 of file G4HENeutronInelastic.cc.

00043 {
00044   outFile << "G4HENeutronInelastic is one of the High Energy\n"
00045           << "Parameterized (HEP) models used to implement inelastic\n"
00046           << "neutron scattering from nuclei.  It is a re-engineered\n"
00047           << "version of the GHEISHA code of H. Fesefeldt.  It divides the\n"
00048           << "initial collision products into backward- and forward-going\n"
00049           << "clusters which are then decayed into final state hadrons.\n"
00050           << "The model does not conserve energy on an event-by-event\n"
00051           << "basis.  It may be applied to neutrons with initial energies\n"
00052           << "above 20 GeV.\n";
00053 }

G4int G4HENeutronInelastic::vecLength

Definition at line 70 of file G4HENeutronInelastic.hh.

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


Public Member Functions
	G4HENeutronInelastic ()
	~G4HENeutronInelastic ()
virtual void	ModelDescription (std::ostream &) const
G4HadFinalState *	ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
G4int	GetNumberOfSecondaries ()
void	FirstIntInCasNeutron (G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, const G4double atomicWeight)
Data Fields
G4int	vecLength

G4HENeutronInelastic Class Reference

Public Member Functions

Data Fields

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Field Documentation