Geant4: G4HESigmaMinusInelastic Class Reference

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

Definition at line 42 of file G4HESigmaMinusInelastic.cc.

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

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

G4HESigmaMinusInelastic::~G4HESigmaMinusInelastic ( ) [inline]

Definition at line 56 of file G4HESigmaMinusInelastic.hh.

00056 {};

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

Implements G4HadronicInteraction.

Definition at line 69 of file G4HESigmaMinusInelastic.cc.

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

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

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

Definition at line 202 of file G4HESigmaMinusInelastic.cc.

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

Referenced by ApplyYourself().

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

G4int G4HESigmaMinusInelastic::GetNumberOfSecondaries ( ) [inline]

Definition at line 65 of file G4HESigmaMinusInelastic.hh.

References vecLength.

00065 {return vecLength;}

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

Reimplemented from G4HadronicInteraction.

Definition at line 55 of file G4HESigmaMinusInelastic.cc.

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

G4int G4HESigmaMinusInelastic::vecLength

Definition at line 60 of file G4HESigmaMinusInelastic.hh.

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


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

G4HESigmaMinusInelastic Class Reference

Public Member Functions

Data Fields

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Field Documentation