Geant4: G4HESigmaPlusInelastic Class Reference

G4HESigmaPlusInelastic::G4HESigmaPlusInelastic ( const G4String & name = "G4HESigmaPlusInelastic" )

Definition at line 43 of file G4HESigmaPlusInelastic.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 G4HESigmaPlusInelastic is being deprecated and will\n"
00052          << "disappear in Geant4 version 10.0"  << G4endl;
00053 }

G4HESigmaPlusInelastic::~G4HESigmaPlusInelastic ( ) [inline]

Definition at line 56 of file G4HESigmaPlusInelastic.hh.

00056 {};

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

Implements G4HadronicInteraction.

Definition at line 70 of file G4HESigmaPlusInelastic.cc.

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

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

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

Definition at line 203 of file G4HESigmaPlusInelastic.cc.

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

Referenced by ApplyYourself().

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

G4int G4HESigmaPlusInelastic::GetNumberOfSecondaries ( ) [inline]

Definition at line 65 of file G4HESigmaPlusInelastic.hh.

References vecLength.

00065 {return vecLength;}

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

Reimplemented from G4HadronicInteraction.

Definition at line 55 of file G4HESigmaPlusInelastic.cc.

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

G4int G4HESigmaPlusInelastic::vecLength

Definition at line 60 of file G4HESigmaPlusInelastic.hh.

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


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

G4HESigmaPlusInelastic Class Reference

Public Member Functions

Data Fields

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Field Documentation