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

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