G4HEOmegaMinusInelastic Class Reference

#include <G4HEOmegaMinusInelastic.hh>

Inheritance diagram for G4HEOmegaMinusInelastic:

G4HEInelastic G4HadronicInteraction

Public Member Functions

 G4HEOmegaMinusInelastic ()
 ~G4HEOmegaMinusInelastic ()
virtual void ModelDescription (std::ostream &) const
G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
G4int GetNumberOfSecondaries ()
void FirstIntInCasOmegaMinus (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 52 of file G4HEOmegaMinusInelastic.hh.


Constructor & Destructor Documentation

G4HEOmegaMinusInelastic::G4HEOmegaMinusInelastic (  )  [inline]

Definition at line 55 of file G4HEOmegaMinusInelastic.hh.

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

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

G4HEOmegaMinusInelastic::~G4HEOmegaMinusInelastic (  )  [inline]

Definition at line 66 of file G4HEOmegaMinusInelastic.hh.

00066 { };


Member Function Documentation

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

Implements G4HadronicInteraction.

Definition at line 57 of file G4HEOmegaMinusInelastic.cc.

References G4HEInelastic::ElasticScattering(), G4HEInelastic::FillParticleChange(), FirstIntInCasOmegaMinus(), 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 << "GHEOmegaMinusInelastic: incident energy < 1 GeV" << G4endl;
00080 
00081   if (verboseLevel > 1) {
00082     G4cout << "G4HEOmegaMinusInelastic::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:  variable set 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   G4bool inElastic = true;
00129   vecLength = 0;
00130         
00131   if (verboseLevel > 1)
00132     G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00133            << incidentCode << G4endl;
00134 
00135   G4bool successful = false; 
00136     
00137   FirstIntInCasOmegaMinus(inElastic, availableEnergy, pv, vecLength,
00138                           incidentParticle, targetParticle, atomicWeight);
00139 
00140   if (verboseLevel > 1)
00141     G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00142 
00143   if ((vecLength > 0) && (availableEnergy > 1.)) 
00144     StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00145                                   pv, vecLength,
00146                                   incidentParticle, targetParticle);
00147 
00148   HighEnergyCascading(successful, pv, vecLength,
00149                       excitationEnergyGNP, excitationEnergyDTA,
00150                       incidentParticle, targetParticle,
00151                       atomicWeight, atomicNumber);
00152   if (!successful)
00153     HighEnergyClusterProduction(successful, pv, vecLength,
00154                                 excitationEnergyGNP, excitationEnergyDTA,
00155                                 incidentParticle, targetParticle,
00156                                 atomicWeight, atomicNumber);
00157   if (!successful) 
00158     MediumEnergyCascading(successful, pv, vecLength, 
00159                           excitationEnergyGNP, excitationEnergyDTA, 
00160                           incidentParticle, targetParticle,
00161                           atomicWeight, atomicNumber);
00162 
00163   if (!successful)
00164     MediumEnergyClusterProduction(successful, pv, vecLength,
00165                                   excitationEnergyGNP, excitationEnergyDTA,       
00166                                   incidentParticle, targetParticle,
00167                                   atomicWeight, atomicNumber);
00168   if (!successful)
00169     QuasiElasticScattering(successful, pv, vecLength,
00170                            excitationEnergyGNP, excitationEnergyDTA,
00171                            incidentParticle, targetParticle, 
00172                            atomicWeight, atomicNumber);
00173   if (!successful)
00174     ElasticScattering(successful, pv, vecLength,
00175                       incidentParticle,    
00176                       atomicWeight, atomicNumber);
00177 
00178   if (!successful)
00179     G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
00180            << G4endl;
00181 
00182   FillParticleChange(pv,  vecLength);
00183   delete [] pv;
00184   theParticleChange.SetStatusChange(stopAndKill);
00185   return &theParticleChange;
00186 }

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

Definition at line 190 of file G4HEOmegaMinusInelastic.cc.

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

Referenced by ApplyYourself().

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

G4int G4HEOmegaMinusInelastic::GetNumberOfSecondaries (  )  [inline]

Definition at line 75 of file G4HEOmegaMinusInelastic.hh.

References vecLength.

00075 {return vecLength;}

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

Reimplemented from G4HadronicInteraction.

Definition at line 42 of file G4HEOmegaMinusInelastic.cc.

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


Field Documentation

G4int G4HEOmegaMinusInelastic::vecLength

Definition at line 70 of file G4HEOmegaMinusInelastic.hh.

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


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