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

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