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

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