G4HEPionMinusInelastic.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // 11-OCT-2007 F.W. Jones: fixed incorrect Imax (should be Imin) in
00029 //             sampling of charge exchange.
00030 
00031 // G4 Process: Gheisha High Energy Collision model.
00032 // This includes the high energy cascading model, the two-body-resonance model
00033 // and the low energy two-body model. Not included are the low energy stuff
00034 // like nuclear reactions, nuclear fission without any cascading and all
00035 // processes for particles at rest.  
00036 // First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96.  
00037 // H. Fesefeldt, RWTH-Aachen, 23-October-1996
00038  
00039 #include "G4HEPionMinusInelastic.hh"
00040 #include "globals.hh"
00041 #include "G4ios.hh"
00042 #include "G4PhysicalConstants.hh"
00043 
00044 void G4HEPionMinusInelastic::ModelDescription(std::ostream& outFile) const
00045 {
00046   outFile << "G4HEPionMinusInelastic is one of the High Energy\n"
00047           << "Parameterized (HEP) models used to implement inelastic\n"
00048           << "pi- scattering from nuclei.  It is a re-engineered\n"
00049           << "version of the GHEISHA code of H. Fesefeldt.  It divides the\n"
00050           << "initial collision products into backward- and forward-going\n"
00051           << "clusters which are then decayed into final state hadrons.\n"
00052           << "The model does not conserve energy on an event-by-event\n"
00053           << "basis.  It may be applied to pi- with initial energies\n"
00054           << "above 20 GeV.\n";
00055 }
00056 
00057 
00058 G4HadFinalState*
00059 G4HEPionMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00060                                       G4Nucleus& targetNucleus)
00061 {
00062   G4HEVector* pv = new G4HEVector[MAXPART];
00063   const G4HadProjectile* aParticle = &aTrack;
00064   const G4double A = targetNucleus.GetA_asInt();
00065   const G4double Z = targetNucleus.GetZ_asInt();
00066   G4HEVector incidentParticle(aParticle);
00067      
00068   G4double atomicNumber = Z;
00069   G4double atomicWeight = A;
00070 
00071   G4int incidentCode = incidentParticle.getCode();
00072   G4double incidentMass = incidentParticle.getMass();
00073   G4double incidentTotalEnergy = incidentParticle.getEnergy();
00074 
00075   // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00076   // DHW 19 May 2011: variable set but not used
00077 
00078   G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
00079 
00080   if (incidentKineticEnergy < 1.)
00081     G4cout << "GHEPionMinusInelastic: incident energy < 1 GeV" << G4endl;
00082 
00083   if (verboseLevel > 1) {
00084     G4cout << "G4HEPionMinusInelastic::ApplyYourself" << G4endl;
00085     G4cout << "incident particle " << incidentParticle.getName()
00086            << "mass "              << incidentMass
00087            << "kinetic energy "    << incidentKineticEnergy
00088            << G4endl;
00089     G4cout << "target material with (A,Z) = (" 
00090            << atomicWeight << "," << atomicNumber << ")" << G4endl;
00091   }
00092 
00093   G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 
00094                                               atomicWeight, atomicNumber);
00095   if (verboseLevel > 1)
00096     G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
00097 
00098   incidentKineticEnergy -= inelasticity;
00099     
00100   G4double excitationEnergyGNP = 0.;
00101   G4double excitationEnergyDTA = 0.; 
00102 
00103   G4double excitation = NuclearExcitation(incidentKineticEnergy,
00104                                           atomicWeight, atomicNumber,
00105                                           excitationEnergyGNP,
00106                                           excitationEnergyDTA);
00107   if (verboseLevel > 1)
00108     G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 
00109            << excitationEnergyDTA << G4endl;
00110 
00111   incidentKineticEnergy -= excitation;
00112   incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00113   // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)                    
00114   //                                   *(incidentTotalEnergy+incidentMass));
00115   // DHW 19 May 2011: variable set but not used
00116 
00117   G4HEVector targetParticle;
00118  
00119   if (G4UniformRand() < atomicNumber/atomicWeight) { 
00120     targetParticle.setDefinition("Proton");
00121   } else { 
00122     targetParticle.setDefinition("Neutron");
00123   }
00124 
00125   G4double targetMass = targetParticle.getMass();
00126   G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
00127                                         + targetMass*targetMass
00128                                         + 2.0*targetMass*incidentTotalEnergy);
00129   G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00130 
00131   // The value of the inElastic flag was originally defined in the Gheisha
00132   // stand-alone code as follows: 
00133   //   G4bool inElastic = InElasticCrossSectionInFirstInt
00134   //                      (availableEnergy, incidentCode, incidentTotalMomentum);  
00135   // For unknown reasons, it was replaced by the following code in Geant
00136 
00137   G4bool inElastic = true;
00138   //    if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false;    
00139 
00140   vecLength = 0;
00141         
00142   if (verboseLevel > 1)
00143     G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00144            << incidentCode << G4endl;
00145 
00146   G4bool successful = false; 
00147     
00148   FirstIntInCasPionMinus(inElastic, availableEnergy, pv, vecLength,
00149                          incidentParticle, targetParticle);
00150 
00151   if (verboseLevel > 1)
00152     G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00153 
00154   if ((vecLength > 0) && (availableEnergy > 1.)) 
00155     StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00156                                   pv, vecLength,
00157                                   incidentParticle, targetParticle);
00158 
00159   HighEnergyCascading(successful, pv, vecLength,
00160                       excitationEnergyGNP, excitationEnergyDTA,
00161                       incidentParticle, targetParticle,
00162                       atomicWeight, atomicNumber);
00163   if (!successful)
00164     HighEnergyClusterProduction(successful, pv, vecLength,
00165                                 excitationEnergyGNP, excitationEnergyDTA,
00166                                 incidentParticle, targetParticle,
00167                                 atomicWeight, atomicNumber);
00168   if (!successful) 
00169     MediumEnergyCascading(successful, pv, vecLength, 
00170                           excitationEnergyGNP, excitationEnergyDTA, 
00171                           incidentParticle, targetParticle,
00172                           atomicWeight, atomicNumber);
00173 
00174   if (!successful)
00175     MediumEnergyClusterProduction(successful, pv, vecLength,
00176                                   excitationEnergyGNP, excitationEnergyDTA,       
00177                                   incidentParticle, targetParticle,
00178                                   atomicWeight, atomicNumber);
00179   if (!successful)
00180     QuasiElasticScattering(successful, pv, vecLength,
00181                            excitationEnergyGNP, excitationEnergyDTA,
00182                            incidentParticle, targetParticle, 
00183                            atomicWeight, atomicNumber);
00184   if (!successful)
00185     ElasticScattering(successful, pv, vecLength,
00186                       incidentParticle,    
00187                       atomicWeight, atomicNumber);
00188 
00189   if (!successful)
00190     G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 
00191            << G4endl;
00192 
00193   FillParticleChange(pv,  vecLength);
00194   delete [] pv;
00195   theParticleChange.SetStatusChange(stopAndKill);
00196   return &theParticleChange;
00197 }
00198 
00199 
00200 void
00201 G4HEPionMinusInelastic::FirstIntInCasPionMinus(G4bool& inElastic,
00202                                                const G4double availableEnergy,
00203                                                G4HEVector pv[],
00204                                                G4int& vecLen,
00205                                                const G4HEVector& incidentParticle,
00206                                                const G4HEVector& targetParticle)
00207 
00208 // Pion- undergoes interaction with nucleon within a nucleus.  Check if it is
00209 // energetically possible to produce pions/kaons.  In not, assume nuclear excitation
00210 // occurs and input particle is degraded in energy. No other particles are produced.
00211 // If reaction is possible, find the correct number of pions/protons/neutrons
00212 // produced using an interpolation to multiplicity data.  Replace some pions or
00213 // protons/neutrons by kaons or strange baryons according to the average
00214 // multiplicity per inelastic reaction.
00215 {
00216   static const G4double expxu = 82.;     // upper bound for arg. of exp
00217   static const G4double expxl = -expxu;  // lower bound for arg. of exp
00218 
00219   static const G4double protb = 0.7;
00220   static const G4double neutb = 0.7;
00221   static const G4double c = 1.25;
00222 
00223   static const G4int numMul = 1200;
00224   static const G4int numSec = 60;
00225 
00226   G4int protonCode = Proton.getCode();
00227   G4int targetCode = targetParticle.getCode();
00228   G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00229 
00230   static G4bool first = true;
00231   static G4double protmul[numMul], protnorm[numSec];  // proton constants
00232   static G4double neutmul[numMul], neutnorm[numSec];  // neutron constants
00233 
00234   // misc. local variables
00235   // npos = number of pi+,  nneg = number of pi-,  nero = number of pi0
00236 
00237   G4int i, counter, nt, npos, nneg, nero;
00238 
00239   if (first) {
00240     // compute normalization constants, this will only be done once
00241     first = false;
00242     for( i=0; i<numMul; i++ )protmul[i]  = 0.0;
00243     for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
00244     counter = -1;
00245        for( npos=0; npos<(numSec/3); npos++ ) 
00246           {
00247             for( nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++ ) 
00248                {
00249                  for( nero=0; nero<numSec/3; nero++ ) 
00250                     {
00251                       if( ++counter < numMul ) 
00252                         {
00253                           nt = npos+nneg+nero;
00254                           if( (nt>0) && (nt<=numSec) ) 
00255                             {
00256                               protmul[counter] =
00257                                     pmltpc(npos,nneg,nero,nt,protb,c) ;
00258                               protnorm[nt-1] += protmul[counter];
00259                             }
00260                         }
00261                     }
00262                }
00263           }
00264        for( i=0; i<numMul; i++ )neutmul[i]  = 0.0;
00265        for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00266        counter = -1;
00267        for( npos=0; npos<numSec/3; npos++ ) 
00268           {
00269             for( nneg=npos; nneg<=(npos+2); nneg++ ) 
00270                {
00271                  for( nero=0; nero<numSec/3; nero++ ) 
00272                     {
00273                       if( ++counter < numMul ) 
00274                         {
00275                           nt = npos+nneg+nero;
00276                           if( (nt>0) && (nt<=numSec) ) 
00277                             {
00278                                neutmul[counter] =
00279                                       pmltpc(npos,nneg,nero,nt,neutb,c);
00280                                neutnorm[nt-1] += neutmul[counter];
00281                             }
00282                         }
00283                     }
00284                }
00285           }
00286        for( i=0; i<numSec; i++ ) 
00287           {
00288             if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00289             if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00290           }
00291      }                                          // end of initialization
00292 
00293          
00294    // initialize the first two places
00295    // the same as beam and target                                    
00296    pv[0] = incidentParticle;
00297    pv[1] = targetParticle;
00298    vecLen = 2;
00299 
00300    if (!inElastic || (availableEnergy <= PionPlus.getMass())) 
00301       return;
00302 
00303                                         
00304    // inelastic scattering
00305 
00306    npos = 0, nneg = 0, nero = 0;
00307    G4double eab = availableEnergy;
00308    G4int ieab = G4int( eab*5.0 );
00309    
00310    G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
00311    if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) ) 
00312      {
00313        // suppress high multiplicity events at low momentum
00314        // only one additional pion will be produced
00315        G4double cech[] = {1., 0.95, 0.79, 0.32, 0.19, 0.16, 0.14, 0.12, 0.10, 0.08};
00316        G4int iplab = Imin(9, G4int( incidentTotalMomentum*5.));
00317        if( G4UniformRand() < cech[iplab] )
00318          {
00319            if( targetCode == protonCode )
00320              {
00321                pv[0] = PionZero;
00322                pv[1] = Neutron;
00323                return;
00324              }
00325            else
00326              {
00327                return;
00328              }
00329          }
00330 
00331        G4double w0, wp, wm, wt, ran;
00332        if( targetCode == protonCode )                    // target is a proton 
00333          {
00334            w0 = - sqr(1.+protb)/(2.*c*c);
00335            wp = w0 = std::exp(w0);
00336            wm = - sqr(-1.+protb)/(2.*c*c);
00337            wm = std::exp(wm);
00338            wp *= 10.;
00339            wt = w0+wp+wm;
00340            wp = w0+wp;
00341            ran = G4UniformRand();
00342            if( ran < w0/wt )  
00343              { npos = 0; nneg = 0; nero = 1; }
00344            else if ( ran < wp/wt )
00345              { npos = 1; nneg = 0; nero = 0; }
00346            else
00347              { npos = 0; nneg = 1; nero = 0; }       
00348          } 
00349        else 
00350          {                                           // target is a neutron
00351            w0 = -sqr(1.+neutb)/(2.*c*c);
00352            wm = -sqr(-1.+neutb)/(2.*c*c);
00353            w0 = std::exp(w0);
00354            wm = std::exp(wm);
00355            if( G4UniformRand() < w0/(w0+wm) )
00356              { npos = 0; nneg = 0; nero = 1; }       
00357            else
00358              { npos = 0; nneg = 1; nero = 0; }       
00359          }
00360      }
00361    else
00362      {
00363        // number of total particles vs. centre of mass Energy - 2*proton mass
00364    
00365        G4double aleab = std::log(availableEnergy);
00366        G4double n     = 3.62567+aleab*(0.665843+aleab*(0.336514
00367                     + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00368    
00369        // normalization constant for kno-distribution.
00370        // calculate first the sum of all constants, check for numerical problems.   
00371        G4double test, dum, anpn = 0.0;
00372 
00373        for (nt=1; nt<=numSec; nt++) {
00374          test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00375          dum = pi*nt/(2.0*n*n);
00376          if (std::fabs(dum) < 1.0) { 
00377            if( test >= 1.0e-10 )anpn += dum*test;
00378          } else { 
00379            anpn += dum*test;
00380          }
00381        }
00382    
00383        G4double ran = G4UniformRand();
00384        G4double excs = 0.0;
00385        if( targetCode == protonCode ) 
00386          {
00387            counter = -1;
00388            for (npos=0; npos<numSec/3; npos++) {
00389              for (nneg=Imax(0,npos-1); nneg<=(npos+1); nneg++) {
00390                for (nero=0; nero<numSec/3; nero++) {
00391                  if (++counter < numMul) {
00392                    nt = npos+nneg+nero;
00393                    if ( (nt>0) && (nt<=numSec) ) {
00394                      test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00395                      dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00396                      if (std::fabs(dum) < 1.0) { 
00397                        if( test >= 1.0e-10 )excs += dum*test;
00398                      } else { 
00399                        excs += dum*test;
00400                      }
00401                      if (ran < excs) goto outOfLoop;      //----------------------->
00402                    }
00403                  }
00404                }
00405              }
00406            }
00407        
00408            // 3 previous loops continued to the end
00409            inElastic = false;                 // quasi-elastic scattering   
00410            return;
00411          }
00412        else   
00413          {                                         // target must be a neutron
00414            counter = -1;
00415            for (npos=0; npos<numSec/3; npos++) {
00416              for (nneg=npos; nneg<=(npos+2); nneg++) {
00417                for (nero=0; nero<numSec/3; nero++) {
00418                  if (++counter < numMul) {
00419                    nt = npos+nneg+nero;
00420                    if ( (nt>=1) && (nt<=numSec) ) {
00421                      test = std::exp( Amin( expxu, Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00422                      dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00423                      if (std::fabs(dum) < 1.0) { 
00424                        if( test >= 1.0e-10 )excs += dum*test;
00425                      } else { 
00426                        excs += dum*test;
00427                      }
00428                      if (ran < excs) goto outOfLoop;     // ----------------->
00429                    }
00430                  }
00431                }
00432              }
00433            }
00434            // 3 previous loops continued to the end
00435            inElastic = false;                     // quasi-elastic scattering.
00436            return;
00437          }
00438      } 
00439    outOfLoop:           //  <-----------------------------------------------   
00440     
00441    if( targetCode == protonCode)
00442      {
00443        if( npos == (1+nneg))
00444          {
00445            pv[1] = Neutron;
00446          }
00447        else if (npos == nneg)
00448          {
00449            if( G4UniformRand() < 0.75)
00450              {
00451              }
00452            else
00453              {
00454                pv[0] = PionZero;
00455                pv[1] = Neutron;
00456              }
00457          }
00458        else      
00459          {
00460            pv[0] = PionZero;
00461          } 
00462      }  
00463    else
00464      {
00465        if( npos == (nneg-1))
00466          {
00467            if( G4UniformRand() < 0.5)
00468              {
00469                pv[1] = Proton;
00470              }
00471            else
00472              {
00473                pv[0] = PionZero;
00474              }
00475          } 
00476        else if ( npos == nneg)
00477          {
00478          }
00479        else
00480          {
00481            pv[0] = PionZero;
00482            pv[1] = Proton;
00483          }
00484      }      
00485 
00486 
00487    nt = npos + nneg + nero;
00488    while ( nt > 0)
00489        {
00490          G4double ran = G4UniformRand();
00491          if ( ran < (G4double)npos/nt)
00492             { 
00493               if( npos > 0 ) 
00494                 { pv[vecLen++] = PionPlus;
00495                   npos--;
00496                 }
00497             }
00498          else if ( ran < (G4double)(npos+nneg)/nt)
00499             {   
00500               if( nneg > 0 )
00501                 { 
00502                   pv[vecLen++] = PionMinus;
00503                   nneg--;
00504                 }
00505             }
00506          else
00507             {
00508               if( nero > 0 )
00509                 { 
00510                   pv[vecLen++] = PionZero;
00511                   nero--;
00512                 }
00513             }
00514          nt = npos + nneg + nero;
00515        } 
00516    if (verboseLevel > 1)
00517       {
00518         G4cout << "Particles produced: " ;
00519         G4cout << pv[0].getName() << " " ;
00520         G4cout << pv[1].getName() << " " ;
00521         for (i=2; i < vecLen; i++)   
00522             { 
00523               G4cout << pv[i].getName() << " " ;
00524             }
00525          G4cout << G4endl;
00526       }
00527    return;
00528  }
00529 
00530 
00531 
00532 
00533 
00534 
00535 
00536 
00537 

Generated on Mon May 27 17:48:29 2013 for Geant4 by  doxygen 1.4.7