G4LEAntiProtonInelastic.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: AntiProton Inelastic Process
00029 // J.L. Chuma, TRIUMF, 13-Feb-1997
00030 // J.P. Wellisch: 23-Apr-97: Bug hunting
00031 // Modified by J.L.Chuma 30-Apr-97: added originalTarget for CalculateMomenta
00032 //
00033  
00034 #include <iostream>
00035 
00036 #include "G4LEAntiProtonInelastic.hh"
00037 #include "G4PhysicalConstants.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "Randomize.hh"
00040 
00041 
00042 G4LEAntiProtonInelastic::G4LEAntiProtonInelastic(const G4String& name)
00043  :G4InelasticInteraction(name)
00044 {
00045   SetMinEnergy(0.0);
00046   SetMaxEnergy(25.*GeV);
00047   G4cout << "WARNING: model G4LEAntiProtonInelastic is being deprecated and will\n"
00048          << "disappear in Geant4 version 10.0"  << G4endl;
00049 }
00050 
00051 
00052 void G4LEAntiProtonInelastic::ModelDescription(std::ostream& outFile) const
00053 {
00054   outFile << "G4LEAntiProtonInelastic is one of the Low Energy Parameterized\n"
00055           << "(LEP) models used to implement inelastic anti-proton scattering\n"
00056           << "from nuclei.  It is a re-engineered version of the GHEISHA\n"
00057           << "code of H. Fesefeldt.  It divides the initial collision\n"
00058           << "products into backward- and forward-going clusters which are\n"
00059           << "then decayed into final state hadrons.  The model does not\n"
00060           << "conserve energy on an event-by-event basis.  It may be\n"
00061           << "applied to anti-protons with initial energies between 0 and 25\n"
00062           << "GeV.\n";
00063 }
00064 
00065 
00066 G4HadFinalState*
00067 G4LEAntiProtonInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00068                                        G4Nucleus& targetNucleus)
00069 { 
00070   const G4HadProjectile* originalIncident = &aTrack;
00071 
00072   if (originalIncident->GetKineticEnergy() <= 0.1*MeV) {
00073     theParticleChange.SetStatusChange(isAlive);
00074     theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00075     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00076     return &theParticleChange;
00077   }
00078 
00079   // create the target particle
00080 
00081   G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
00082     
00083   if (verboseLevel > 1) {
00084     const G4Material *targetMaterial = aTrack.GetMaterial();
00085     G4cout << "G4LEAntiProtonInelastic::ApplyYourself called" << G4endl;
00086     G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00087     G4cout << "target material = " << targetMaterial->GetName() << ", ";
00088     G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00089            << G4endl;
00090   }
00091 
00092   // Fermi motion and evaporation
00093   // As of Geant3, the Fermi energy calculation had not been Done
00094 
00095   G4double ek = originalIncident->GetKineticEnergy()/MeV;
00096   G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00097   G4ReactionProduct modifiedOriginal;
00098   modifiedOriginal = *originalIncident;
00099     
00100   G4double tkin = targetNucleus.Cinema( ek );
00101   ek += tkin;
00102   modifiedOriginal.SetKineticEnergy( ek*MeV );
00103   G4double et = ek + amas;
00104   G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00105   G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
00106   if (pp > 0.0) {
00107     G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00108     modifiedOriginal.SetMomentum( momentum * (p/pp) );
00109   }
00110 
00111   // calculate black track energies
00112   tkin = targetNucleus.EvaporationEffects( ek );
00113   ek -= tkin;
00114   modifiedOriginal.SetKineticEnergy( ek*MeV );
00115   et = ek + amas;
00116   p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00117   pp = modifiedOriginal.GetMomentum().mag()/MeV;
00118   if (pp > 0.0) {
00119     G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00120     modifiedOriginal.SetMomentum( momentum * (p/pp) );
00121   }
00122   G4ReactionProduct currentParticle = modifiedOriginal;
00123   G4ReactionProduct targetParticle;
00124   targetParticle = *originalTarget;
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;
00135   const G4double anni = std::min( 1.3*originalIncident->GetTotalMomentum()/GeV, 0.4 );
00136     
00137   if ((currentParticle.GetKineticEnergy()/MeV > cutOff) ||
00138       (G4UniformRand() > anni) )
00139     Cascade(vec, vecLen, originalIncident, currentParticle, targetParticle,
00140             incidentHasChanged, targetHasChanged, quasiElastic);
00141   else
00142     quasiElastic = true;
00143     
00144   CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
00145                    modifiedOriginal, targetNucleus, currentParticle,
00146                    targetParticle, incidentHasChanged, targetHasChanged,
00147                    quasiElastic);
00148     
00149   SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
00150 
00151   if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00152 
00153   delete originalTarget;
00154   return &theParticleChange;
00155 }
00156 
00157 
00158 void G4LEAntiProtonInelastic::Cascade(
00159    G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00160    G4int& vecLen,
00161    const G4HadProjectile *originalIncident,
00162    G4ReactionProduct &currentParticle,
00163    G4ReactionProduct &targetParticle,
00164    G4bool &incidentHasChanged,
00165    G4bool &targetHasChanged,
00166    G4bool &quasiElastic )
00167 {
00168   // derived from original FORTRAN code CASPB by H. Fesefeldt (13-Sep-1987)
00169   //
00170   // AntiProton undergoes interaction with nucleon within a nucleus.  Check if
00171   // it is energetically possible to produce pions/kaons.  In not, assume
00172   // nuclear excitation occurs and input particle is degraded in energy. No
00173   // other particles are produced.  If reaction is possible, find the correct
00174   // number of pions/protons/neutrons produced using an interpolation to
00175   // multiplicity data.  Replace some pions or protons/neutrons by kaons or
00176   // strange baryons according to the average  multiplicity per inelastic
00177   // reaction.
00178 
00179   const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00180   const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
00181   const G4double pOriginal = originalIncident->GetTotalMomentum()/MeV;
00182   const G4double targetMass = targetParticle.GetMass()/MeV;
00183   G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
00184                                         targetMass*targetMass +
00185                                         2.0*targetMass*etOriginal );
00186   G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00187     
00188   static G4bool first = true;
00189   const G4int numMul = 1200;
00190   const G4int numMulA = 400;
00191   const G4int numSec = 60;
00192   static G4double protmul[numMul], protnorm[numSec]; // proton constants
00193   static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
00194   static G4double protmulA[numMulA], protnormA[numSec]; // proton constants
00195   static G4double neutmulA[numMulA], neutnormA[numSec]; // neutron constants
00196 
00197   // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
00198   G4int counter, nt=0;
00199   G4int npos=0, nneg = 0, nzero=0;
00200   G4double test;
00201   const G4double c = 1.25;    
00202   const G4double b[] = { 0.7, 0.7 };
00203   if (first) {  // Computation of normalization constants will only be done once
00204     first = false;
00205     G4int i;
00206     for( i=0; i<numMul; ++i )protmul[i] = 0.0;
00207     for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
00208     counter = -1;
00209     for (npos = 0; npos < (numSec/3); ++npos) {
00210       for (nneg = std::max(0,npos-1); nneg <= (npos+1); ++nneg) {
00211         for (nzero = 0; nzero<numSec/3; ++nzero) {
00212           if ( ++counter < numMul ) {
00213             nt = npos+nneg+nzero;
00214             if (nt > 0 && nt <= numSec) {
00215               protmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00216               protnorm[nt-1] += protmul[counter];
00217             }
00218           }
00219         }
00220       }
00221     }
00222 
00223     for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
00224     for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
00225     counter = -1;
00226     for (npos = 0; npos < numSec/3; ++npos) {
00227       for (nneg = npos; nneg <= (npos+2); ++nneg) {
00228         for (nzero = 0; nzero < numSec/3; ++nzero) {
00229           if (++counter < numMul) {
00230             nt = npos+nneg+nzero;
00231             if ((nt > 0) && (nt <= numSec) ) {
00232               neutmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00233               neutnorm[nt-1] += neutmul[counter];
00234             }
00235           }
00236         }
00237       }
00238     }
00239     for (i = 0; i < numSec; ++i) {
00240       if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00241       if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00242     }
00243 
00244     // do the same for annihilation channels
00245     for (i = 0; i < numMulA; ++i) protmulA[i] = 0.0;
00246     for (i = 0; i < numSec; ++i) protnormA[i] = 0.0;
00247     counter = -1;
00248     for (npos = 0; npos < (numSec/3); ++npos) {
00249       nneg = npos;
00250       for (nzero=0; nzero<numSec/3; ++nzero) {
00251         if (++counter < numMulA) {
00252           nt = npos+nneg+nzero;
00253           if (nt > 1 && nt <= numSec) {
00254             protmulA[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00255             protnormA[nt-1] += protmulA[counter];
00256           }
00257         }
00258       }
00259     }
00260     for (i = 0; i < numMulA; ++i) neutmulA[i] = 0.0;
00261     for (i = 0; i < numSec; ++i) neutnormA[i] = 0.0;
00262     counter = -1;
00263     for (npos=0; npos < numSec/3; ++npos) {
00264       nneg = npos+1;
00265       for( nzero=0; nzero<numSec/3; ++nzero ) {
00266         if( ++counter < numMulA ) {
00267           nt = npos+nneg+nzero;
00268           if( (nt>1) && (nt<=numSec) ) {
00269             neutmulA[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00270             neutnormA[nt-1] += neutmulA[counter];
00271           }
00272         }
00273       }
00274     }
00275     for (i=0; i<numSec; ++i ) {
00276       if( protnormA[i] > 0.0 )protnormA[i] = 1.0/protnormA[i];
00277       if( neutnormA[i] > 0.0 )neutnormA[i] = 1.0/neutnormA[i];
00278     }
00279   }   // end of initialization
00280 
00281   const G4double expxu = 82.;           // upper bound for arg. of exp
00282   const G4double expxl = -expxu;        // lower bound for arg. of exp
00283     
00284   G4ParticleDefinition *aNeutron      = G4Neutron::Neutron();
00285   G4ParticleDefinition *aProton       = G4Proton::Proton();
00286   G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
00287   G4ParticleDefinition *aPiPlus       = G4PionPlus::PionPlus();
00288     
00289   const G4double anhl[] = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.90,
00290                            0.6,0.52,0.47,0.44,0.41,0.39,0.37,0.35,0.34,0.24,
00291                            0.19,0.15,0.12,0.10,0.09,0.07,0.06,0.05,0.0};
00292     
00293   G4int iplab = G4int( pOriginal/GeV*10.0 );
00294   if( iplab >  9 )iplab = G4int( pOriginal/GeV ) + 9;
00295   if( iplab > 18 )iplab = G4int( pOriginal/GeV/10.0 ) + 18;
00296   if( iplab > 27 )iplab = 28;
00297   if (G4UniformRand() > anhl[iplab]) {
00298     if (availableEnergy <= aPiPlus->GetPDGMass()/MeV) {
00299       quasiElastic = true;
00300       return;
00301     }
00302     G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
00303     const G4double supp[] = {0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98};
00304     G4double w0, wp, wt, wm;
00305     if ((availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) ) {
00306       // suppress high multiplicity events at low momentum
00307       // only one pion will be produced
00308       npos = nneg = nzero = 0;
00309       if (targetParticle.GetDefinition() == aProton) {
00310         test = std::exp(std::min(expxu,
00311                         std::max(expxl, -(1.0+b[1])*(1.0+b[1])/(2.0*c*c) ) ) );
00312         w0 = test;
00313         wp = test;        
00314         test = std::exp(std::min(expxu,
00315                         std::max(expxl, -(-1.0+b[1])*(-1.0+b[1])/(2.0*c*c) ) ) );
00316         wm = test;
00317         wt = w0+wp+wm;
00318         wp += w0;
00319         G4double ran = G4UniformRand();
00320         if( ran < w0/wt )
00321             nzero = 1;
00322           else if( ran < wp/wt )
00323             npos = 1;
00324           else
00325             nneg = 1;
00326       } else {
00327         // target is a neutron
00328         test = std::exp(std::min(expxu,
00329                         std::max(expxl, -(1.0+b[0])*(1.0+b[0])/(2.0*c*c) ) ) );
00330         w0 = test;
00331         test = std::exp(std::min(expxu,
00332                         std::max(expxl, -(-1.0+b[0])*(-1.0+b[0])/(2.0*c*c) ) ) );
00333         wm = test;
00334         G4double ran = G4UniformRand();
00335         if (ran < w0/(w0+wm) )
00336           nzero = 1;
00337         else
00338           nneg = 1;
00339       }
00340     } else {
00341       // (availableEnergy >= 2.0*GeV) || (random number < supp[ieab])
00342         G4double n, anpn;
00343         GetNormalizationConstant( availableEnergy, n, anpn );
00344         G4double ran = G4UniformRand();
00345         G4double dum, excs = 0.0;
00346         if (targetParticle.GetDefinition() == aProton) {
00347           counter = -1;
00348           for( npos=0; npos<numSec/3 && ran>=excs; ++npos ) {
00349             for( nneg=std::max(0,npos-1); nneg<=(npos+1) && ran>=excs; ++nneg) {
00350               for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00351               {
00352                 if( ++counter < numMul )
00353                 {
00354                   nt = npos+nneg+nzero;
00355                   if( (nt>0) && (nt<=numSec) )
00356                   {
00357                     test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00358                     dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00359                     if( std::fabs(dum) < 1.0 )
00360                     {
00361                       if( test >= 1.0e-10 )excs += dum*test;
00362                     }
00363                     else
00364                       excs += dum*test;
00365                   }
00366                 }
00367               }
00368             }
00369           }
00370           if( ran >= excs )  // 3 previous loops continued to the end
00371           {
00372             quasiElastic = true;
00373             return;
00374           }
00375         }
00376         else  // target must be a neutron
00377         {
00378           counter = -1;
00379           for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00380           {
00381             for (nneg = npos; nneg <= (npos+2) && ran >= excs; ++nneg) {
00382               for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00383               {
00384                 if( ++counter < numMul )
00385                 {
00386                   nt = npos+nneg+nzero;
00387                   if( (nt>0) && (nt<=numSec) )
00388                   {
00389                     test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00390                     dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00391                     if( std::fabs(dum) < 1.0 )
00392                     {
00393                       if( test >= 1.0e-10 )excs += dum*test;
00394                     }
00395                     else
00396                       excs += dum*test;
00397                   }
00398                 }
00399               }
00400             }
00401           }
00402           if( ran >= excs )  // 3 previous loops continued to the end
00403           {
00404             quasiElastic = true;
00405             return;
00406           }
00407         }
00408         npos--; nneg--; nzero--;
00409       }
00410       if( targetParticle.GetDefinition() == aProton )
00411       {
00412         switch (npos-nneg)
00413         {
00414          case 0:
00415            if( G4UniformRand() < 0.33 )
00416            {
00417              currentParticle.SetDefinitionAndUpdateE( anAntiNeutron );
00418              targetParticle.SetDefinitionAndUpdateE( aNeutron );
00419              incidentHasChanged = true;
00420              targetHasChanged = true;
00421            }
00422            break;
00423          case 1:
00424            targetParticle.SetDefinitionAndUpdateE( aNeutron );
00425            targetHasChanged = true;
00426            break;
00427          default:
00428            currentParticle.SetDefinitionAndUpdateE( anAntiNeutron );
00429            incidentHasChanged = true;
00430            break;
00431         }
00432       }
00433       else  // target must be a neutron
00434       {
00435         switch (npos-nneg)
00436         {
00437          case -1:
00438            if( G4UniformRand() < 0.5 )
00439            {
00440              targetParticle.SetDefinitionAndUpdateE( aProton );
00441              targetHasChanged = true;
00442            }
00443            else
00444            {
00445              currentParticle.SetDefinitionAndUpdateE( anAntiNeutron );
00446              incidentHasChanged = true;
00447            }
00448            break;
00449          case 0:
00450            break;
00451          default:
00452            currentParticle.SetDefinitionAndUpdateE( anAntiNeutron );
00453            targetParticle.SetDefinitionAndUpdateE( aProton );
00454            incidentHasChanged = true;
00455            targetHasChanged = true;
00456            break;
00457         }
00458       }
00459     }
00460     else   //  random number <= anhl[iplab]
00461     {
00462       if( centerofmassEnergy <= 2*aPiPlus->GetPDGMass()/MeV )
00463       {
00464         quasiElastic = true;
00465         return;
00466       }
00467       //
00468       // annihilation channels
00469       //
00470       G4double n, anpn;
00471       GetNormalizationConstant( -centerofmassEnergy, n, anpn );
00472       G4double ran = G4UniformRand();
00473       G4double dum, excs = 0.0;
00474       if( targetParticle.GetDefinition() == aProton )
00475       {
00476         counter = -1;
00477         for( npos=0; (npos<numSec/3) && (ran>=excs); ++npos )
00478         {
00479           nneg = npos;
00480           for( nzero=0; (nzero<numSec/3) && (ran>=excs); ++nzero )
00481           {
00482             if( ++counter < numMulA )
00483             {
00484               nt = npos+nneg+nzero;
00485               if( (nt>1) && (nt<=numSec) )
00486               {
00487                 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00488                 dum = (pi/anpn)*nt*protmulA[counter]*protnormA[nt-1]/(2.0*n*n);
00489                 if( std::abs(dum) < 1.0 )
00490                 {
00491                   if( test >= 1.0e-10 )excs += dum*test;
00492                 }
00493                 else
00494                   excs += dum*test;
00495               }
00496             }
00497           }
00498         }
00499       }
00500       else  // target must be a neutron
00501       {
00502         counter = -1;
00503         for( npos=0; (npos<numSec/3) && (ran>=excs); ++npos )
00504         {
00505           nneg = npos+1;
00506           for( nzero=0; (nzero<numSec/3) && (ran>=excs); ++nzero )
00507           {
00508             if( ++counter < numMulA )
00509             {
00510               nt = npos+nneg+nzero;
00511               if( (nt>1) && (nt<=numSec) )
00512               {
00513                 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00514                 dum = (pi/anpn)*nt*neutmulA[counter]*neutnormA[nt-1]/(2.0*n*n);
00515                 if( std::fabs(dum) < 1.0 )
00516                 {
00517                   if( test >= 1.0e-10 )excs += dum*test;
00518                 }
00519                 else
00520                   excs += dum*test;
00521               }
00522             }
00523           }
00524         }
00525     }
00526     if (ran >= excs) {
00527       // 3 previous loops continued to the end
00528       quasiElastic = true;
00529       return;
00530     }
00531     npos--; nzero--;
00532     currentParticle.SetMass( 0.0 );
00533     targetParticle.SetMass( 0.0 );
00534   }
00535 
00536   while(npos+nneg+nzero<3) nzero++;
00537   SetUpPions( npos, nneg, nzero, vec, vecLen );
00538   return;
00539 }

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