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

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