G4HELambdaInelastic.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 // G4 Process: Gheisha High Energy Collision model.
00030 // This includes the high energy cascading model, the two-body-resonance model
00031 // and the low energy two-body model. Not included are the low energy stuff
00032 // like nuclear reactions, nuclear fission without any cascading and all
00033 // processes for particles at rest.  
00034 // First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96.  
00035 // H. Fesefeldt, RWTH-Aachen, 23-October-1996
00036  
00037 #include "G4HELambdaInelastic.hh"
00038 #include "globals.hh"
00039 #include "G4ios.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 
00043 G4HELambdaInelastic::G4HELambdaInelastic(const G4String& name)
00044  : G4HEInelastic(name)
00045 {
00046   vecLength = 0;
00047   theMinEnergy = 20*GeV;
00048   theMaxEnergy = 10*TeV;
00049   MAXPART      = 2048;
00050   verboseLevel = 0;
00051   G4cout << "WARNING: model G4HELambdaInelastic is being deprecated and will\n"
00052          << "disappear in Geant4 version 10.0"  << G4endl;  
00053 }
00054 
00055 
00056 void G4HELambdaInelastic::ModelDescription(std::ostream& outFile) const
00057 {
00058   outFile << "G4HELambdaInelastic is one of the High Energy Parameterized\n"
00059           << "(HEP) models used to implement inelastic Lambda scattering\n"
00060           << "from nuclei.  It is a re-engineered version of the GHEISHA\n"
00061           << "code of H. Fesefeldt.  It divides the initial collision\n"
00062           << "products into backward- and forward-going clusters which are\n"
00063           << "then decayed into final state hadrons.  The model does not\n"
00064           << "conserve energy on an event-by-event basis.  It may be\n"
00065           << "applied to lambdas with initial energies above 20 GeV.\n";
00066 }
00067 
00068 
00069 G4HadFinalState*
00070 G4HELambdaInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00071                                    G4Nucleus& targetNucleus)
00072 {
00073   G4HEVector* pv = new G4HEVector[MAXPART];
00074   const G4HadProjectile* aParticle = &aTrack;
00075   const G4double A = targetNucleus.GetA_asInt();
00076   const G4double Z = targetNucleus.GetZ_asInt();
00077   G4HEVector incidentParticle(aParticle);
00078      
00079   G4double atomicNumber = Z;
00080   G4double atomicWeight = A;
00081 
00082   G4int incidentCode = incidentParticle.getCode();
00083   G4double incidentMass = incidentParticle.getMass();
00084   G4double incidentTotalEnergy = incidentParticle.getEnergy();
00085 
00086   // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00087   // DHW 19 May 2011: variable set but not used
00088 
00089   G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
00090 
00091   if (incidentKineticEnergy < 1.)
00092     G4cout << "GHELambdaInelastic: incident energy < 1 GeV" << G4endl;
00093 
00094   if (verboseLevel > 1) {
00095     G4cout << "G4HELambdaInelastic::ApplyYourself" << G4endl;
00096     G4cout << "incident particle " << incidentParticle.getName()
00097            << "mass "              << incidentMass
00098            << "kinetic energy "    << incidentKineticEnergy
00099            << G4endl;
00100     G4cout << "target material with (A,Z) = (" 
00101            << atomicWeight << "," << atomicNumber << ")" << G4endl;
00102   }
00103 
00104   G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 
00105                                               atomicWeight, atomicNumber);
00106   if (verboseLevel > 1)
00107     G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
00108     
00109   incidentKineticEnergy -= inelasticity;
00110     
00111   G4double excitationEnergyGNP = 0.;
00112   G4double excitationEnergyDTA = 0.; 
00113 
00114   G4double excitation = NuclearExcitation(incidentKineticEnergy,
00115                                           atomicWeight, atomicNumber,
00116                                           excitationEnergyGNP,
00117                                           excitationEnergyDTA);
00118   if (verboseLevel > 1)
00119     G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 
00120            << excitationEnergyDTA << G4endl;             
00121 
00122   incidentKineticEnergy -= excitation;
00123   incidentTotalEnergy = incidentKineticEnergy + incidentMass;
00124   // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)                    
00125   //                                   *(incidentTotalEnergy+incidentMass));
00126   // DHW 19 May 2011: variable set but not used
00127 
00128   G4HEVector targetParticle;
00129   if (G4UniformRand() < atomicNumber/atomicWeight) { 
00130     targetParticle.setDefinition("Proton");
00131   } else { 
00132     targetParticle.setDefinition("Neutron");
00133   }
00134 
00135   G4double targetMass = targetParticle.getMass();
00136   G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
00137                                         + targetMass*targetMass
00138                                         + 2.0*targetMass*incidentTotalEnergy);
00139   G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
00140 
00141   G4bool inElastic = true;
00142   vecLength = 0;           
00143         
00144   if (verboseLevel > 1)
00145     G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
00146            << incidentCode << G4endl;
00147 
00148   G4bool successful = false; 
00149     
00150   FirstIntInCasLambda(inElastic, availableEnergy, pv, vecLength,
00151                       incidentParticle, targetParticle, atomicWeight);
00152 
00153   if (verboseLevel > 1)
00154     G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
00155 
00156   if ((vecLength > 0) && (availableEnergy > 1.)) 
00157     StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
00158                                   pv, vecLength,
00159                                   incidentParticle, targetParticle);
00160 
00161   HighEnergyCascading(successful, pv, vecLength,
00162                       excitationEnergyGNP, excitationEnergyDTA,
00163                       incidentParticle, targetParticle,
00164                       atomicWeight, atomicNumber);
00165   if (!successful)
00166     HighEnergyClusterProduction(successful, pv, vecLength,
00167                                 excitationEnergyGNP, excitationEnergyDTA,
00168                                 incidentParticle, targetParticle,
00169                                 atomicWeight, atomicNumber);
00170   if (!successful) 
00171     MediumEnergyCascading(successful, pv, vecLength, 
00172                           excitationEnergyGNP, excitationEnergyDTA, 
00173                           incidentParticle, targetParticle,
00174                           atomicWeight, atomicNumber);
00175 
00176   if (!successful)
00177     MediumEnergyClusterProduction(successful, pv, vecLength,
00178                                   excitationEnergyGNP, excitationEnergyDTA,       
00179                                   incidentParticle, targetParticle,
00180                                   atomicWeight, atomicNumber);
00181   if (!successful)
00182     QuasiElasticScattering(successful, pv, vecLength,
00183                            excitationEnergyGNP, excitationEnergyDTA,
00184                            incidentParticle, targetParticle, 
00185                            atomicWeight, atomicNumber);
00186   if (!successful)
00187     ElasticScattering(successful, pv, vecLength,
00188                       incidentParticle,    
00189                       atomicWeight, atomicNumber);
00190 
00191   if (!successful)
00192     G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
00193            << G4endl;
00194   FillParticleChange(pv, vecLength);
00195   delete [] pv;
00196   theParticleChange.SetStatusChange(stopAndKill);
00197   return &theParticleChange;
00198 }
00199 
00200 
00201 void
00202 G4HELambdaInelastic::FirstIntInCasLambda(G4bool& inElastic,
00203                                          const G4double availableEnergy,
00204                                          G4HEVector pv[],
00205                                          G4int& vecLen,
00206                                          const G4HEVector& incidentParticle,
00207                                          const G4HEVector& targetParticle,
00208                                          const G4double atomicWeight)
00209 
00210 // Lambda undergoes interaction with nucleon within a nucleus.  Check if it is
00211 // energetically possible to produce pions/kaons.  In not, assume nuclear 
00212 // excitation occurs and input particle is degraded in energy. No other 
00213 // particles are produced.  If reaction is possible, find the correct number 
00214 // of pions/protons/neutrons produced using an interpolation to multiplicity 
00215 // data.  Replace some pions or protons/neutrons by kaons or strange baryons 
00216 // according to the average multiplicity per inelastic reaction.
00217 {
00218   static const G4double expxu = 82.;     // upper bound for arg. of exp
00219   static const G4double expxl = -expxu;  // lower bound for arg. of exp
00220 
00221   static const G4double protb = 0.7;
00222   static const G4double neutb = 0.7;              
00223   static const G4double c = 1.25;
00224 
00225   static const G4int numMul = 1200;
00226   static const G4int numSec = 60;
00227 
00228   G4int protonCode  = Proton.getCode();
00229   G4int targetCode = targetParticle.getCode();
00230   G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
00231 
00232   static G4bool first = true;
00233   static G4double protmul[numMul], protnorm[numSec];  // proton constants
00234   static G4double neutmul[numMul], neutnorm[numSec];  // neutron constants
00235 
00236   //  misc. local variables
00237   //  npos = number of pi+,  nneg = number of pi-,  nzero = number of pi0
00238 
00239   G4int i, counter, nt, npos, nneg, nzero;
00240 
00241   if (first) {  // Computation of normalization constants will only be done once
00242     first = false;
00243     for( i=0; i<numMul; i++ )protmul[i]  = 0.0;
00244     for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
00245     counter = -1;
00246     for (npos = 0; npos < (numSec/3); npos++) {
00247       for (nneg = std::max(0,npos-2); nneg <= (npos+1); nneg++) {
00248         for (nzero = 0; nzero < numSec/3; nzero++) {
00249           if (++counter < numMul) {
00250             nt = npos+nneg+nzero;
00251             if ((nt>0) && (nt<=numSec) ) {
00252               protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
00253               protnorm[nt-1] += protmul[counter];
00254             }
00255           }
00256         }
00257       }
00258     }
00259 
00260     for( i=0; i<numMul; i++ )neutmul[i]  = 0.0;
00261     for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
00262     counter = -1;
00263     for (npos = 0; npos < numSec/3; npos++) {
00264             for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ ) 
00265                {
00266                  for( nzero=0; nzero<numSec/3; nzero++ ) 
00267                     {
00268                       if( ++counter < numMul ) 
00269                         {
00270                           nt = npos+nneg+nzero;
00271                           if( (nt>0) && (nt<=numSec) ) 
00272                             {
00273                                neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
00274                                neutnorm[nt-1] += neutmul[counter];
00275                             }
00276                         }
00277                     }
00278                }
00279           }
00280        for( i=0; i<numSec; i++ ) 
00281           {
00282             if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00283             if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00284           }
00285   }   // end of initialization
00286          
00287   pv[0] = incidentParticle;  // initialize the first two places
00288   pv[1] = targetParticle;    // the same as beam and target
00289   vecLen = 2;
00290 
00291   if (!inElastic) {  // quasi-elastic scattering, no pions produced
00292     G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
00293     G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*2.5 ) );
00294     if (G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) {
00295       G4double ran = G4UniformRand();
00296       if (targetCode == protonCode) { 
00297                if( ran < 0.2)
00298                  {
00299                    pv[0] = SigmaPlus;
00300                    pv[1] = Neutron;
00301                  }
00302                else if(ran < 0.4)
00303                  {
00304                    pv[0] = SigmaZero;
00305                  }
00306                else if(ran < 0.6)
00307                  {
00308                    pv[0] = Proton;
00309                    pv[1] = Lambda;
00310                  }
00311                else if(ran < 0.8)
00312                  {
00313                    pv[0] = Proton;
00314                    pv[1] = SigmaZero;
00315                  }
00316                else
00317                  {
00318                    pv[0] = Neutron;
00319                    pv[1] = SigmaPlus;
00320                  }
00321       } else { 
00322                if(ran < 0.2)
00323                  {
00324                    pv[0] = SigmaZero;
00325                  }
00326                else if(ran < 0.4)
00327                  {
00328                    pv[0] = SigmaMinus;
00329                    pv[1] = Proton;
00330                  }
00331                else if(ran < 0.6)
00332                  {
00333                    pv[0] = Neutron;
00334                    pv[1] = Lambda;
00335                  } 
00336                else if(ran < 0.8)
00337                  {
00338                    pv[0] = Neutron;
00339                    pv[1] = SigmaZero;
00340                  }
00341                else
00342                  {
00343                    pv[0] = Proton;
00344                    pv[1] = SigmaMinus;
00345                  }  
00346       }
00347     }
00348     return;
00349   }
00350   else if (availableEnergy <= PionPlus.getMass())
00351     return;
00352 
00353   // inelastic scattering
00354   npos = 0; nneg = 0; nzero = 0;
00355 
00356   // number of total particles vs. centre of mass Energy - 2*proton mass
00357    
00358   G4double aleab = std::log(availableEnergy);
00359   G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
00360                 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
00361    
00362   // normalization constant for kno-distribution.
00363   // calculate first the sum of all constants, check for numerical problems.  
00364   G4double test, dum, anpn = 0.0;
00365 
00366   for (nt=1; nt<=numSec; nt++) {
00367      test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00368      dum = pi*nt/(2.0*n*n);
00369      if (std::fabs(dum) < 1.0) {
00370        if( test >= 1.0e-10 )anpn += dum*test;
00371      } else { 
00372        anpn += dum*test;
00373      }
00374   }
00375    
00376   G4double ran = G4UniformRand();
00377   G4double excs = 0.0;
00378   if (targetCode == protonCode) {
00379     counter = -1;
00380     for (npos = 0; npos < numSec/3; npos++) {
00381       for (nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++) {
00382            for (nzero=0; nzero<numSec/3; nzero++) {
00383              if (++counter < numMul) {
00384                nt = npos+nneg+nzero;
00385                if ( (nt>0) && (nt<=numSec) ) {
00386                  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00387                  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00388                  if (std::fabs(dum) < 1.0) { 
00389                    if( test >= 1.0e-10 )excs += dum*test;
00390                  } else { 
00391                    excs += dum*test;
00392                  }
00393                  if (ran < excs) goto outOfLoop;      //--------------------->
00394                }   
00395              }
00396            }
00397       }
00398     }
00399     // 3 previous loops continued to the end
00400 
00401     inElastic = false;                 // quasi-elastic scattering   
00402     return;
00403 
00404   } else {            // target must be a neutron
00405     counter = -1;
00406     for (npos=0; npos<numSec/3; npos++) {
00407       for (nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++) {
00408         for (nzero=0; nzero<numSec/3; nzero++) {
00409           if (++counter < numMul) {
00410             nt = npos+nneg+nzero;
00411             if ( (nt>=1) && (nt<=numSec) ) {
00412               test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00413               dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00414               if (std::fabs(dum) < 1.0) { 
00415                 if( test >= 1.0e-10 )excs += dum*test;
00416               } else { 
00417                 excs += dum*test;
00418               }
00419               if (ran < excs) goto outOfLoop;       // ------------->
00420             }
00421           }
00422         }
00423       }
00424     }
00425                            // 3 previous loops continued to the end
00426     inElastic = false;     // quasi-elastic scattering.
00427     return;
00428   }
00429  
00430   outOfLoop:           //  <--------------------------------------------   
00431     
00432   ran = G4UniformRand();
00433   if (targetCode == protonCode) {
00434     if (npos == nneg) { 
00435            if (ran < 0.25)
00436              { 
00437              }
00438            else if(ran < 0.5)
00439              {
00440                pv[0] = SigmaZero;
00441              }  
00442            else 
00443              {
00444                pv[0] = SigmaPlus;   
00445                pv[1] = Neutron;
00446              }
00447     } else if (npos == (nneg+1)) {
00448            if( G4UniformRand() < 0.25)
00449              {
00450                pv[1] = Neutron;
00451              }
00452            else if(ran < 0.5)
00453              {
00454                pv[0] = SigmaZero;
00455                pv[1] = Neutron; 
00456              }
00457            else
00458              {
00459                pv[0] = SigmaMinus;
00460              }
00461     } else if (npos == (nneg-1)) {
00462       pv[0] = SigmaPlus;
00463     } else {
00464       pv[0] = SigmaMinus;
00465       pv[1] = Neutron;
00466     }
00467  
00468   } else {
00469     if (npos == nneg) {
00470            if(ran < 0.5) 
00471              {
00472              }
00473            else
00474              {
00475                pv[0] = SigmaMinus;
00476                pv[1] = Proton;
00477              }
00478     } else if (npos == (nneg-1)) {  
00479            if( ran < 0.25)
00480              {
00481                pv[1] = Proton;
00482              }
00483            else if(ran < 0.5)               
00484              {
00485                pv[0] = SigmaZero;
00486                pv[1] = Proton;
00487              }
00488            else
00489              {
00490                pv[1] = SigmaPlus;
00491              }
00492     } else if (npos == (1+nneg)) {
00493       pv[0] = SigmaMinus;
00494     } else {
00495       pv[0] = SigmaPlus; 
00496       pv[1] = Proton;
00497     }
00498   }      
00499 
00500   nt = npos + nneg + nzero;
00501   while (nt > 0) {
00502     G4double rnd = G4UniformRand();
00503     if (rnd < (G4double)npos/nt) { 
00504       if (npos > 0) {
00505         pv[vecLen++] = PionPlus;
00506         npos--;
00507       }
00508     } else if (rnd < (G4double)(npos+nneg)/nt) {   
00509       if (nneg > 0) { 
00510         pv[vecLen++] = PionMinus;
00511         nneg--;
00512       }
00513     } else {
00514       if (nzero > 0) { 
00515         pv[vecLen++] = PionZero;
00516         nzero--;
00517       }
00518     }
00519     nt = npos + nneg + nzero;
00520   }
00521 
00522   if (verboseLevel > 1) {
00523     G4cout << "Particles produced: " ;
00524     G4cout << pv[0].getName() << " " ;
00525     G4cout << pv[1].getName() << " " ;
00526     for (i = 2; i < vecLen; i++) G4cout << pv[i].getName() << " " ;
00527     G4cout << G4endl;
00528   }
00529   return;
00530 }
00531 

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