G4HEPionPlusInelastic Class Reference

#include <G4HEPionPlusInelastic.hh>

Inheritance diagram for G4HEPionPlusInelastic:

G4HEInelastic G4HadronicInteraction

Public Member Functions

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


Constructor & Destructor Documentation

G4HEPionPlusInelastic::G4HEPionPlusInelastic (  )  [inline]

Definition at line 55 of file G4HEPionPlusInelastic.hh.

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

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

G4HEPionPlusInelastic::~G4HEPionPlusInelastic (  )  [inline]

Definition at line 66 of file G4HEPionPlusInelastic.hh.

00066 {};


Member Function Documentation

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

Implements G4HadronicInteraction.

Definition at line 57 of file G4HEPionPlusInelastic.cc.

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

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

Definition at line 191 of file G4HEPionPlusInelastic.cc.

References G4HEInelastic::Amax(), G4HEInelastic::Amin(), G4cout, G4endl, G4UniformRand, G4HEVector::getCode(), G4HEVector::getMass(), G4HEVector::getName(), G4HEVector::getTotalMomentum(), G4HEInelastic::Imax(), CLHEP::detail::n, G4HEInelastic::Neutron, neutronCode, G4INCL::Math::pi, G4HEInelastic::PionMinus, G4HEInelastic::PionPlus, G4HEInelastic::PionZero, G4HEInelastic::pmltpc(), G4HEInelastic::Proton, sqr(), and G4HEInelastic::verboseLevel.

Referenced by ApplyYourself().

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

G4int G4HEPionPlusInelastic::GetNumberOfSecondaries (  )  [inline]

Definition at line 75 of file G4HEPionPlusInelastic.hh.

References vecLength.

00075 {return vecLength;}

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

Reimplemented from G4HadronicInteraction.

Definition at line 42 of file G4HEPionPlusInelastic.cc.

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


Field Documentation

G4int G4HEPionPlusInelastic::vecLength

Definition at line 70 of file G4HEPionPlusInelastic.hh.

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