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

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