Geant4: G4HEProtonInelastic Class Reference

G4HEProtonInelastic::G4HEProtonInelastic ( ) [inline]

Definition at line 55 of file G4HEProtonInelastic.hh.

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

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

G4HEProtonInelastic::~G4HEProtonInelastic ( ) [inline]

Definition at line 66 of file G4HEProtonInelastic.hh.

00066 {};

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

Implements G4HadronicInteraction.

Definition at line 57 of file G4HEProtonInelastic.cc.

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

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

Definition at line 201 of file G4HEProtonInelastic.cc.

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

G4int G4HEProtonInelastic::GetNumberOfSecondaries ( ) [inline]

Definition at line 75 of file G4HEProtonInelastic.hh.

References vecLength.

00075 {return vecLength;}

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

Reimplemented from G4HadronicInteraction.

Definition at line 42 of file G4HEProtonInelastic.cc.

00043 {
00044   outFile << "G4HEProtonInelastic is one of the High Energy\n"
00045           << "Parameterized (HEP) models used to implement inelastic\n"
00046           << "proton 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 protons with initial energies\n"
00052           << "above 20 GeV.\n";
00053 }

G4int G4HEProtonInelastic::vecLength

Definition at line 70 of file G4HEProtonInelastic.hh.

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


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

G4HEProtonInelastic Class Reference

Public Member Functions

Data Fields

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Field Documentation