G4HEXiMinusInelastic Class Reference

#include <G4HEXiMinusInelastic.hh>

Inheritance diagram for G4HEXiMinusInelastic:

G4HEInelastic G4HadronicInteraction

Public Member Functions

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


Constructor & Destructor Documentation

G4HEXiMinusInelastic::G4HEXiMinusInelastic (  )  [inline]

Definition at line 55 of file G4HEXiMinusInelastic.hh.

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

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

G4HEXiMinusInelastic::~G4HEXiMinusInelastic (  )  [inline]

Definition at line 66 of file G4HEXiMinusInelastic.hh.

00066 {};


Member Function Documentation

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

Implements G4HadronicInteraction.

Definition at line 56 of file G4HEXiMinusInelastic.cc.

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

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

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

Definition at line 189 of file G4HEXiMinusInelastic.cc.

References G4cout, G4endl, G4UniformRand, G4HEVector::getCode(), G4HEVector::getMass(), G4HEVector::getName(), 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, G4HEInelastic::verboseLevel, G4HEInelastic::XiMinus, and G4HEInelastic::XiZero.

Referenced by ApplyYourself().

00204 {
00205   static const G4double expxu = 82.;     // upper bound for arg. of exp
00206   static const G4double expxl = -expxu;  // lower bound for arg. of exp
00207 
00208   static const G4double protb = 0.7;
00209   static const G4double neutb = 0.7;              
00210   static const G4double c = 1.25;
00211 
00212   static const G4int numMul = 1200;
00213   static const G4int numSec = 60;
00214 
00215   G4int neutronCode = Neutron.getCode();
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+1); 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+2); 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 == neutronCode)
00297              {
00298                if (ran < 0.2)
00299                  {
00300                    pv[0] = SigmaMinus;
00301                    pv[1] = SigmaZero;
00302                  } 
00303                else if (ran < 0.4)
00304                  {
00305                    pv[0] = SigmaZero;
00306                    pv[1] = SigmaMinus;
00307                  }
00308                else if (ran < 0.6)
00309                  {
00310                    pv[0] = SigmaMinus;
00311                    pv[1] = Lambda;
00312                  }
00313                else if (ran < 0.8)
00314                  {
00315                    pv[0] = Lambda;
00316                    pv[1] = SigmaMinus;
00317                  }
00318                else
00319                  {
00320                    pv[0] = Neutron;
00321                    pv[1] = XiMinus;
00322                  }         
00323              }              
00324            else
00325              { 
00326                if (ran < 0.2)
00327                  {
00328                    pv[0] = SigmaZero;
00329                    pv[1] = SigmaZero;  
00330                  }
00331                else if (ran < 0.3)
00332                  {
00333                    pv[0] = Lambda;
00334                    pv[1] = Lambda;
00335                  }
00336                else if (ran < 0.4)
00337                  {
00338                    pv[0] = SigmaZero;
00339                    pv[1] = Lambda;
00340                  } 
00341                else if (ran < 0.5)
00342                  {
00343                    pv[0] = Lambda;
00344                    pv[1] = SigmaZero;
00345                  }
00346                else if (ran < 0.7)
00347                  {
00348                    pv[0] = SigmaZero;
00349                    pv[1] = Neutron;
00350                  }  
00351                else if (ran < 0.8)
00352                  {
00353                    pv[0] = Neutron;
00354                    pv[1] = SigmaZero;
00355                  }
00356                else
00357                  {
00358                    pv[0] = Proton;
00359                    pv[1] = XiMinus;
00360                  }  
00361              }
00362          }
00363        return;
00364      }
00365    else if (availableEnergy <= PionPlus.getMass())
00366        return;
00367 
00368                                                   //   inelastic scattering
00369 
00370    npos = 0; nneg = 0; nzero = 0;
00371 
00372 //                    number of total particles vs. centre of mass Energy - 2*proton mass
00373    
00374    G4double aleab = std::log(availableEnergy);
00375    G4double n     = 3.62567+aleab*(0.665843+aleab*(0.336514
00376                     + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00377    
00378 //                    normalization constant for kno-distribution.
00379 //                    calculate first the sum of all constants, check for numerical problems.             
00380    G4double test, dum, anpn = 0.0;
00381 
00382    for (nt=1; nt<=numSec; nt++) {
00383      test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00384      dum = pi*nt/(2.0*n*n);
00385      if (std::fabs(dum) < 1.0) {
00386        if (test >= 1.0e-10) anpn += dum*test;
00387      } else { 
00388        anpn += dum*test;
00389      }
00390    }
00391    
00392    G4double ran = G4UniformRand();
00393    G4double excs = 0.0;
00394    if( targetCode == protonCode ) 
00395      {
00396        counter = -1;
00397        for (npos=0; npos<numSec/3; npos++) {
00398          for (nneg=std::max(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>0) && (nt<=numSec) ) {
00403                  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00404                  dum = (pi/anpn)*nt*protmul[counter]*protnorm[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        
00417        // 3 previous loops continued to the end
00418        inElastic = false;                 // quasi-elastic scattering   
00419        return;
00420      }
00421    else   
00422      {                                         // target must be a neutron
00423        counter = -1;
00424        for (npos=0; npos<numSec/3; npos++) {
00425          for (nneg=npos; nneg<=(npos+2); nneg++) { 
00426            for( nzero=0; nzero<numSec/3; nzero++) {
00427              if (++counter < numMul) {
00428                nt = npos+nneg+nzero;
00429                if ( (nt>=1) && (nt<=numSec) ) {
00430                  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00431                  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00432                  if (std::fabs(dum) < 1.0) { 
00433                    if( test >= 1.0e-10 )excs += dum*test;
00434                  } else { 
00435                    excs += dum*test;
00436                  }
00437                  if (ran < excs) goto outOfLoop;       // --------------------->
00438                }
00439              }
00440            }
00441          }
00442        }
00443        // 3 previous loops continued to the end
00444        inElastic = false;                         // quasi-elastic scattering.
00445        return;
00446      }
00447  
00448   outOfLoop:           //  <---------------------------------------------------   
00449     
00450   // in the following we do not consider
00451   // strangeness transfer in high multiplicity
00452   // events. YK combinations are added in
00453   // StrangeParticlePairProduction
00454   ran = G4UniformRand(); 
00455   if (targetCode == neutronCode) {
00456        if( npos == nneg)
00457          { 
00458          }
00459        else if (npos == (nneg-1))
00460          {
00461            if( ran < 0.50)
00462              {
00463                pv[0] = XiZero;
00464              }
00465            else
00466              {
00467                pv[1] = Proton;
00468              }
00469          }
00470        else
00471          {
00472            pv[0] = XiZero;
00473            pv[1] = Proton;
00474          }   
00475   } else {
00476        if (npos == nneg)
00477          {
00478            if (ran < 0.5) 
00479              {
00480              }
00481            else
00482              {
00483                pv[0] = XiZero;
00484                pv[1] = Neutron;
00485              }  
00486          }  
00487        else if (npos == (nneg+1))
00488          {  
00489            pv[1] = Neutron;
00490          } 
00491        else 
00492          {
00493            pv[0] = XiZero;
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 G4HEXiMinusInelastic::GetNumberOfSecondaries (  )  [inline]

Definition at line 75 of file G4HEXiMinusInelastic.hh.

References vecLength.

00075 {return vecLength;}

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

Reimplemented from G4HadronicInteraction.

Definition at line 41 of file G4HEXiMinusInelastic.cc.

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


Field Documentation

G4int G4HEXiMinusInelastic::vecLength

Definition at line 70 of file G4HEXiMinusInelastic.hh.

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


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