G4LEPionPlusInelastic.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 // Hadronic Process: PionPlus Inelastic Process
00029 // J.L. Chuma, TRIUMF, 19-Nov-1996
00030 
00031 // Modified by J.L.Chuma 30-Apr-97: added originalTarget for CalculateMomenta
00032 // fixing charge exchange - HPW Sep 2002.
00033  
00034 #include <iostream>
00035 
00036 #include "G4LEPionPlusInelastic.hh"
00037 #include "G4PhysicalConstants.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "Randomize.hh"
00040 
00041 G4LEPionPlusInelastic::G4LEPionPlusInelastic(const G4String& name)
00042  :G4InelasticInteraction(name)
00043 {
00044   SetMinEnergy(0.0);
00045   SetMaxEnergy(55.*GeV);
00046   G4cout << "WARNING: model G4LEPionPlusInelastic is being deprecated and will\n"
00047          << "disappear in Geant4 version 10.0"  << G4endl;
00048 }
00049 
00050 
00051 void G4LEPionPlusInelastic::ModelDescription(std::ostream& outFile) const
00052 {
00053   outFile << "G4LEPionPlusInelastic is one of the Low Energy Parameterized\n"
00054           << "(LEP) models used to implement inelastic pi+ scattering\n"
00055           << "from nuclei.  It is a re-engineered version of the GHEISHA\n"
00056           << "code of H. Fesefeldt.  It divides the initial collision\n"
00057           << "products into backward- and forward-going clusters which are\n"
00058           << "then decayed into final state hadrons.  The model does not\n"
00059           << "conserve energy on an event-by-event basis.  It may be\n"
00060           << "applied to pions with initial energies between 0 and 25\n"
00061           << "GeV.\n";
00062 }
00063 
00064 
00065 G4HadFinalState*
00066 G4LEPionPlusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00067                                      G4Nucleus& targetNucleus)
00068 {
00069   const G4HadProjectile *originalIncident = &aTrack;
00070   if (originalIncident->GetKineticEnergy()<= 0.1*MeV) {
00071     theParticleChange.SetStatusChange(isAlive);
00072     theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00073     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 
00074     return &theParticleChange;      
00075   }
00076 
00077   // create the target particle
00078     
00079   G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
00080   G4ReactionProduct targetParticle( originalTarget->GetDefinition() );
00081     
00082   if (verboseLevel > 1) {
00083     const G4Material* targetMaterial = aTrack.GetMaterial();
00084     G4cout << "G4LEPionPlusInelastic::ApplyYourself called" << G4endl;
00085     G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy() << "MeV, ";
00086     G4cout << "target material = " << targetMaterial->GetName() << ", ";
00087     G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00088            << G4endl;
00089   }
00090   G4ReactionProduct currentParticle( 
00091   const_cast<G4ParticleDefinition *>(originalIncident->GetDefinition() ) );
00092   currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
00093   currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
00094     
00095   // Fermi motion and evaporation
00096   // As of Geant3, the Fermi energy calculation had not been Done  
00097   G4double ek = originalIncident->GetKineticEnergy();
00098   G4double amas = originalIncident->GetDefinition()->GetPDGMass();
00099     
00100   G4double tkin = targetNucleus.Cinema(ek);
00101   ek += tkin;
00102   currentParticle.SetKineticEnergy(ek);
00103   G4double et = ek + amas;
00104   G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00105   G4double pp = currentParticle.GetMomentum().mag();
00106   if (pp > 0.0) {
00107     G4ThreeVector momentum = currentParticle.GetMomentum();
00108     currentParticle.SetMomentum( momentum * (p/pp) );
00109   }
00110     
00111   // calculate black track energies  
00112   tkin = targetNucleus.EvaporationEffects(ek);
00113   ek -= tkin;
00114   currentParticle.SetKineticEnergy(ek);
00115   et = ek + amas;
00116   p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00117   pp = currentParticle.GetMomentum().mag();
00118   if (pp > 0.0) {
00119     G4ThreeVector momentum = currentParticle.GetMomentum();
00120     currentParticle.SetMomentum( momentum * (p/pp) );
00121   }
00122 
00123   G4ReactionProduct modifiedOriginal = currentParticle;
00124 
00125   currentParticle.SetSide(1); // incident always goes in forward hemisphere
00126   targetParticle.SetSide(-1);  // target always goes in backward hemisphere
00127   G4bool incidentHasChanged = false;
00128   G4bool targetHasChanged = false;
00129   G4bool quasiElastic = false;
00130   G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;  // vec will contain the secondary particles
00131   G4int vecLen = 0;
00132   vec.Initialize(0);
00133     
00134   const G4double cutOff = 0.1*MeV;
00135   if (currentParticle.GetKineticEnergy() > cutOff)
00136     Cascade(vec, vecLen, originalIncident, currentParticle, targetParticle,
00137             incidentHasChanged, targetHasChanged, quasiElastic);
00138     
00139   CalculateMomenta(vec, vecLen,
00140                    originalIncident, originalTarget, modifiedOriginal,
00141                    targetNucleus, currentParticle, targetParticle,
00142                    incidentHasChanged, targetHasChanged, quasiElastic);
00143     
00144   SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
00145 
00146   if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00147 
00148   delete originalTarget;
00149   return &theParticleChange;
00150 }
00151 
00152  
00153 void G4LEPionPlusInelastic::Cascade(
00154    G4FastVector<G4ReactionProduct,GHADLISTSIZE>& vec,
00155    G4int& vecLen,
00156    const G4HadProjectile* originalIncident,
00157    G4ReactionProduct& currentParticle,
00158    G4ReactionProduct& targetParticle,
00159    G4bool& incidentHasChanged,
00160    G4bool& targetHasChanged,
00161    G4bool& quasiElastic)
00162 {
00163   // derived from original FORTRAN code CASPIP by H. Fesefeldt (18-Sep-1987)
00164   //
00165   // pi+  undergoes interaction with nucleon within nucleus.
00166   // Check if energetically possible to produce pions/kaons.
00167   // If not assume nuclear excitation occurs and input particle
00168   // is degraded in energy.  No other particles produced.
00169   // If reaction is possible find correct number of pions/protons/neutrons
00170   // produced using an interpolation to multiplicity data.
00171   // Replace some pions or protons/neutrons by kaons or strange baryons
00172   // according to average multiplicity per inelastic reactions.
00173 
00174   const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass();
00175   const G4double etOriginal = originalIncident->GetTotalEnergy();
00176   const G4double pOriginal = originalIncident->GetTotalMomentum();
00177   const G4double targetMass = targetParticle.GetMass();
00178   G4double centerofmassEnergy = std::sqrt(mOriginal*mOriginal +
00179                                           targetMass*targetMass +
00180                                           2.0*targetMass*etOriginal);
00181   G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00182   static G4bool first = true;
00183   const G4int numMul = 1200;
00184   const G4int numSec = 60;
00185   static G4double protmul[numMul], protnorm[numSec]; // proton constants
00186   static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
00187 
00188   // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
00189   G4int counter, nt=0, npos=0, nneg=0, nzero=0;
00190   const G4double c = 1.25;
00191   const G4double b[] = { 0.70, 0.70 };
00192     if( first ) {      // compute normalization constants, this will only be Done once
00193       first = false;
00194       G4int i;
00195       for( i=0; i<numMul; ++i )protmul[i] = 0.0;
00196       for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
00197       counter = -1;
00198       for( npos=0; npos<(numSec/3); ++npos ) {
00199         for( nneg=std::max(0,npos-2); nneg<=npos; ++nneg ) {
00200           for( nzero=0; nzero<numSec/3; ++nzero ) {
00201             if( ++counter < numMul ) {
00202               nt = npos+nneg+nzero;
00203               if( nt > 0 ) {
00204                 protmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00205                 protnorm[nt-1] += protmul[counter];
00206               }
00207             }
00208           }
00209         }
00210       }
00211       for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
00212       for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
00213       counter = -1;
00214       for( npos=0; npos<numSec/3; ++npos ) {
00215         for( nneg=std::max(0,npos-1); nneg<=(npos+1); ++nneg ) {
00216           for( nzero=0; nzero<numSec/3; ++nzero ) {
00217             if( ++counter < numMul ) {
00218               nt = npos+nneg+nzero;
00219               if( (nt>0) && (nt<=numSec) ) {
00220                 neutmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00221                 neutnorm[nt-1] += neutmul[counter];
00222               }
00223             }
00224           }
00225         }
00226       }
00227       for( i=0; i<numSec; ++i ) {
00228         if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00229         if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00230       }
00231     }   // end of initialization
00232     
00233     const G4double expxu = 82.;           // upper bound for arg. of exp
00234     const G4double expxl = -expxu;        // lower bound for arg. of exp
00235     G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00236     G4ParticleDefinition *aProton = G4Proton::Proton();
00237     G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
00238     G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
00239     const G4double supp[] = {0.,0.2,0.45,0.55,0.65,0.75,0.85,0.90,0.94,0.98};
00240     G4double test, w0, wp, wt, wm;
00241     if( (availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) )
00242     {
00243       // suppress high multiplicity events at low momentum
00244       // only one pion will be produced
00245       // charge exchange reaction is included in inelastic cross section
00246       
00247       const G4double cech[] = {1.,0.95,0.79,0.32,0.19,0.16,0.14,0.12,0.10,0.08};
00248       G4int iplab = G4int(std::min( 9.0, pOriginal/GeV*5.0 ));
00249       if( G4UniformRand() <= cech[iplab] )
00250       {
00251         if( targetParticle.GetDefinition() == aNeutron )
00252         {
00253           currentParticle.SetDefinitionAndUpdateE( aPiZero );  // charge exchange
00254           targetParticle.SetDefinitionAndUpdateE( aProton );
00255           incidentHasChanged = true;
00256           targetHasChanged = true;
00257         }
00258       }
00259       
00260       if( availableEnergy <= G4PionMinus::PionMinus()->GetPDGMass() )
00261       {
00262         quasiElastic = true;
00263         return;
00264       }
00265       
00266       nneg = npos = nzero = 0;
00267       if( targetParticle.GetDefinition() == aProton ) {
00268         test = std::exp( std::min( expxu, std::max( expxl, -sqr(1.0+b[0])/(2.0*c*c) ) ) );
00269         w0 = test;
00270         wp = test;        
00271         if( G4UniformRand() < w0/(w0+wp) )
00272           nzero =1;
00273         else
00274           npos = 1;
00275       } else { // target is a neutron
00276         test = std::exp( std::min( expxu, std::max( expxl, -sqr(1.0+b[1])/(2.0*c*c) ) ) );
00277         w0 = test;
00278         wp = test;        
00279         test = std::exp( std::min( expxu, std::max( expxl, -sqr(-1.0+b[1])/(2.0*c*c) ) ) );
00280         wm = test;
00281         wt = w0+wp+wm;
00282         wp = w0+wp;
00283         G4double ran = G4UniformRand();
00284         if( ran < w0/wt )
00285           nzero = 1;
00286         else if( ran < wp/wt )
00287           npos = 1;
00288         else
00289           nneg = 1;
00290       }
00291     } else {
00292       if( availableEnergy <= G4PionMinus::PionMinus()->GetPDGMass() )
00293       {
00294         quasiElastic = true;
00295         return;
00296       }
00297       G4double n, anpn;
00298       GetNormalizationConstant( availableEnergy, n, anpn );
00299       G4double ran = G4UniformRand();
00300       G4double dum, excs = 0.0;
00301       if( targetParticle.GetDefinition() == aProton ) {
00302         counter = -1;
00303         for( npos=0; (npos<numSec/3) && (ran>=excs); ++npos ) {
00304           for( nneg=std::max(0,npos-2); (nneg<=npos) && (ran>=excs); ++nneg ) {
00305             for( nzero=0; (nzero<numSec/3) && (ran>=excs); ++nzero ) {
00306               if( ++counter < numMul ) {
00307                 nt = npos+nneg+nzero;
00308                 if( nt > 0 ) {
00309                   test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00310                   dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00311                   if( std::fabs(dum) < 1.0 ) {
00312                     if( test >= 1.0e-10 )excs += dum*test;
00313                   } else {
00314                     excs += dum*test;
00315                   }
00316                 }
00317               }
00318             }
00319           }
00320         }
00321         if( ran >= excs )
00322         {
00323           quasiElastic = true;
00324           return;  // 3 previous loops continued to the end
00325         }
00326         npos--; nneg--; nzero--;
00327       } else { // target must be a neutron
00328         counter = -1;
00329         for( npos=0; (npos<numSec/3) && (ran>=excs); ++npos ) {
00330           for( nneg=std::max(0,npos-1); (nneg<=(npos+1)) && (ran>=excs); ++nneg ) {
00331             for( nzero=0; (nzero<numSec/3) && (ran>=excs); ++nzero ) {
00332               if( ++counter < numMul ) {
00333                 nt = npos+nneg+nzero;
00334                 if( (nt>=1) && (nt<=numSec) ) {
00335                   test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00336                   dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00337                   if( std::fabs(dum) < 1.0 ) {
00338                     if( test >= 1.0e-10 )excs += dum*test;
00339                   } else {
00340                     excs += dum*test;
00341                   }
00342                 }
00343               }
00344             }
00345           }
00346         }
00347         if( ran >= excs )  // 3 previous loops continued to the end
00348         {
00349           quasiElastic = true;
00350           return;  // 3 previous loops continued to the end
00351         }
00352         npos--; nneg--; nzero--;
00353       }
00354     }
00355     if( targetParticle.GetDefinition() == aProton ) {
00356       switch( npos-nneg ) {
00357        case 1:
00358          if( G4UniformRand() < 0.5 ) {
00359            currentParticle.SetDefinitionAndUpdateE( aPiZero );
00360            incidentHasChanged = true;
00361          } else {
00362            targetParticle.SetDefinitionAndUpdateE( aNeutron );
00363            targetHasChanged = true;
00364          }
00365          break;
00366        case 2:
00367          currentParticle.SetDefinitionAndUpdateE( aPiZero );
00368          targetParticle.SetDefinitionAndUpdateE( aNeutron );
00369          incidentHasChanged = true;
00370          targetHasChanged = true;
00371          break;
00372        default:
00373          break;
00374       }
00375     } else {
00376       switch( npos-nneg ) {
00377        case 0:
00378          if( G4UniformRand() < 0.25 ) {
00379            currentParticle.SetDefinitionAndUpdateE( aPiZero );
00380            targetParticle.SetDefinitionAndUpdateE( aProton );
00381            incidentHasChanged = true;
00382            targetHasChanged = true;
00383          }
00384          break;
00385        case 1:
00386          currentParticle.SetDefinitionAndUpdateE( aPiZero );
00387          incidentHasChanged = true;
00388          break;
00389        default:
00390          targetParticle.SetDefinitionAndUpdateE( aProton );
00391          targetHasChanged = true;
00392          break;
00393       }
00394     }
00395   SetUpPions( npos, nneg, nzero, vec, vecLen );
00396   return;
00397 }

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