G4HEXiZeroInelastic Class Reference

#include <G4HEXiZeroInelastic.hh>

Inheritance diagram for G4HEXiZeroInelastic:

G4HEInelastic G4HadronicInteraction

Public Member Functions

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


Constructor & Destructor Documentation

G4HEXiZeroInelastic::G4HEXiZeroInelastic (  )  [inline]

Definition at line 55 of file G4HEXiZeroInelastic.hh.

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

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

G4HEXiZeroInelastic::~G4HEXiZeroInelastic (  )  [inline]

Definition at line 66 of file G4HEXiZeroInelastic.hh.

00066 {};


Member Function Documentation

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

Implements G4HadronicInteraction.

Definition at line 56 of file G4HEXiZeroInelastic.cc.

References G4HEInelastic::ElasticScattering(), G4HEInelastic::FillParticleChange(), FirstIntInCasXiZero(), 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 << "GHEXiZeroInelastic: incident energy < 1 GeV" << G4endl;
00079 
00080   if (verboseLevel > 1) {
00081     G4cout << "G4HEXiZeroInelastic::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: variables 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   FirstIntInCasXiZero(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 G4HEXiZeroInelastic::FirstIntInCasXiZero ( G4bool inElastic,
const G4double  availableEnergy,
G4HEVector  pv[],
G4int vecLen,
const G4HEVector incidentParticle,
const G4HEVector targetParticle,
const G4double  atomicWeight 
)

Definition at line 189 of file G4HEXiZeroInelastic.cc.

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

G4int G4HEXiZeroInelastic::GetNumberOfSecondaries (  )  [inline]

Definition at line 75 of file G4HEXiZeroInelastic.hh.

References vecLength.

00075 {return vecLength;}

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

Reimplemented from G4HadronicInteraction.

Definition at line 41 of file G4HEXiZeroInelastic.cc.

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


Field Documentation

G4int G4HEXiZeroInelastic::vecLength

Definition at line 70 of file G4HEXiZeroInelastic.hh.

Referenced by ApplyYourself(), G4HEXiZeroInelastic(), 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