G4LEAntiSigmaMinusInelastic.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: AntiSigmaMinus Inelastic Process
00029 // J.L. Chuma, TRIUMF, 19-Feb-1997
00030 // J.P. Wellisch: 25.Apr-97: counter errors removed lines 426, 447
00031 // Modified by J.L.Chuma 30-Apr-97: added originalTarget for CalculateMomenta
00032  
00033 #include "G4LEAntiSigmaMinusInelastic.hh"
00034 #include "G4PhysicalConstants.hh"
00035 #include "G4SystemOfUnits.hh"
00036 #include "Randomize.hh"
00037 
00038 void G4LEAntiSigmaMinusInelastic::ModelDescription(std::ostream& outFile) const
00039 {
00040   outFile << "G4LEAntiSigmaMinusInelastic is one of the Low Energy\n"
00041           << "Parameterized (LEP) models used to implement inelastic\n"
00042           << "antiSigma- scattering from nuclei.  It is a re-engineered\n"
00043           << "version of the GHEISHA code of H. Fesefeldt.  It divides the\n"
00044           << "initial collision products into backward- and forward-going\n"
00045           << "clusters which are then decayed into final state hadrons.  The\n"
00046           << "model does not conserve energy on an event-by-event basis.  It\n"
00047           << "may be applied to antiSigma- with initial energies between 0\n"
00048           << "and 25 GeV.\n";
00049 }
00050 
00051 
00052 G4HadFinalState*
00053 G4LEAntiSigmaMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00054                                            G4Nucleus& targetNucleus)
00055 { 
00056   const G4HadProjectile *originalIncident = &aTrack;
00057   if (originalIncident->GetKineticEnergy()<= 0.1*MeV) {
00058     theParticleChange.SetStatusChange(isAlive);
00059     theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00060     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 
00061     return &theParticleChange;      
00062   }
00063 
00064   // create the target particle
00065   G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
00066     
00067   if (verboseLevel > 1) {
00068     const G4Material *targetMaterial = aTrack.GetMaterial();
00069     G4cout << "G4LEAntiSigmaMinusInelastic::ApplyYourself called" << G4endl;
00070     G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00071     G4cout << "target material = " << targetMaterial->GetName() << ", ";
00072     G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00073            << G4endl;
00074   }
00075 
00076   // Fermi motion and evaporation
00077   // As of Geant3, the Fermi energy calculation had not been Done
00078   G4double ek = originalIncident->GetKineticEnergy()/MeV;
00079   G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00080   G4ReactionProduct modifiedOriginal;
00081   modifiedOriginal = *originalIncident;
00082     
00083   G4double tkin = targetNucleus.Cinema( ek );
00084   ek += tkin;
00085   modifiedOriginal.SetKineticEnergy( ek*MeV );
00086   G4double et = ek + amas;
00087   G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00088   G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
00089   if (pp > 0.0) {
00090     G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00091     modifiedOriginal.SetMomentum( momentum * (p/pp) );
00092   }
00093 
00094   // calculate black track energies
00095   tkin = targetNucleus.EvaporationEffects( ek );
00096   ek -= tkin;
00097   modifiedOriginal.SetKineticEnergy( ek*MeV );
00098   et = ek + amas;
00099   p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00100   pp = modifiedOriginal.GetMomentum().mag()/MeV;
00101   if (pp > 0.0) {
00102     G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00103     modifiedOriginal.SetMomentum( momentum * (p/pp) );
00104   }
00105   G4ReactionProduct currentParticle = modifiedOriginal;
00106   G4ReactionProduct targetParticle;
00107   targetParticle = *originalTarget;
00108   currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
00109   targetParticle.SetSide( -1 );  // target always goes in backward hemisphere
00110   G4bool incidentHasChanged = false;
00111   G4bool targetHasChanged = false;
00112   G4bool quasiElastic = false;
00113   G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;  // vec will contain the secondary particles
00114   G4int vecLen = 0;
00115   vec.Initialize(0);
00116     
00117   const G4double cutOff = 0.1;
00118   const G4double anni = std::min( 1.3*currentParticle.GetTotalMomentum()/GeV, 0.4 );
00119   if ((currentParticle.GetKineticEnergy()/MeV > cutOff) ||
00120       (G4UniformRand() > anni) )
00121     Cascade(vec, vecLen, originalIncident, currentParticle, targetParticle,
00122             incidentHasChanged, targetHasChanged, quasiElastic);
00123     
00124   CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
00125                    modifiedOriginal, targetNucleus, currentParticle,
00126                    targetParticle, incidentHasChanged, targetHasChanged,
00127                    quasiElastic);
00128     
00129   SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
00130 
00131   if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00132 
00133   delete originalTarget;
00134   return &theParticleChange;
00135 }
00136 
00137 
00138 void G4LEAntiSigmaMinusInelastic::Cascade(
00139    G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00140    G4int& vecLen,
00141    const G4HadProjectile *originalIncident,
00142    G4ReactionProduct &currentParticle,
00143    G4ReactionProduct &targetParticle,
00144    G4bool &incidentHasChanged,
00145    G4bool &targetHasChanged,
00146    G4bool &quasiElastic )
00147 {
00148   // derived from original FORTRAN code CASASM by H. Fesefeldt (13-Sep-1987)
00149   //
00150   // AntiSigmaMinus undergoes interaction with nucleon within a nucleus.  Check if it is
00151   // energetically possible to produce pions/kaons.  In not, assume nuclear excitation
00152   // occurs and input particle is degraded in energy. No other particles are produced.
00153   // If reaction is possible, find the correct number of pions/protons/neutrons
00154   // produced using an interpolation to multiplicity data.  Replace some pions or
00155   // protons/neutrons by kaons or strange baryons according to the average
00156   // multiplicity per Inelastic reaction.
00157 
00158   const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00159   const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
00160   const G4double pOriginal = originalIncident->GetTotalMomentum()/MeV;
00161   const G4double targetMass = targetParticle.GetMass()/MeV;
00162   G4double centerofmassEnergy = std::sqrt(mOriginal*mOriginal +
00163                                           targetMass*targetMass +
00164                                           2.0*targetMass*etOriginal);
00165   G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00166     
00167   static G4bool first = true;
00168   const G4int numMul = 1200;
00169   const G4int numMulA = 400;
00170   const G4int numSec = 60;
00171   static G4double protmul[numMul], protnorm[numSec]; // proton constants
00172   static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
00173   static G4double protmulA[numMulA], protnormA[numSec]; // proton constants
00174   static G4double neutmulA[numMulA], neutnormA[numSec]; // neutron constants
00175 
00176   // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
00177   G4int counter, nt=0, npos=0, nneg=0, nzero=0;
00178   G4double test;
00179   const G4double c = 1.25;    
00180   const G4double b[2] = { 0.7, 0.7 };
00181     if( first )       // compute normalization constants, this will only be Done once
00182     {
00183       first = false;
00184       G4int i;
00185       for( i=0; i<numMul; ++i )protmul[i] = 0.0;
00186       for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
00187       counter = -1;
00188       for( npos=0; npos<(numSec/3); ++npos )
00189       {
00190         for( nneg=std::max(0,npos-2); nneg<=npos; ++nneg )
00191         {
00192           for( nzero=0; nzero<numSec/3; ++nzero )
00193           {
00194             if( ++counter < numMul )
00195             {
00196               nt = npos+nneg+nzero;
00197               if( nt>0 && nt<=numSec )
00198               {
00199                 protmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00200                 protnorm[nt-1] += protmul[counter];
00201               }
00202             }
00203           }
00204         }
00205       }
00206       for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
00207       for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
00208       counter = -1;
00209       for( npos=0; npos<numSec/3; ++npos )
00210       {
00211         for( nneg=std::max(0,npos-1); nneg<=(npos+1); ++nneg )
00212         {
00213           for( nzero=0; nzero<numSec/3; ++nzero )
00214           {
00215             if( ++counter < numMul )
00216             {
00217               nt = npos+nneg+nzero;
00218               if( nt>0 && nt<=numSec )
00219               {
00220                 neutmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00221                 neutnorm[nt-1] += neutmul[counter];
00222               }
00223             }
00224           }
00225         }
00226       }
00227       for( i=0; i<numSec; ++i )
00228       {
00229         if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00230         if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00231       }
00232       //
00233       // do the same for annihilation channels
00234       //
00235       for( i=0; i<numMulA; ++i )protmulA[i] = 0.0;
00236       for( i=0; i<numSec; ++i )protnormA[i] = 0.0;
00237       counter = -1;
00238       for( npos=2; npos<(numSec/3); ++npos )
00239       {
00240         nneg = npos-2;
00241         for( nzero=0; nzero<numSec/3; ++nzero )
00242         {
00243           if( ++counter < numMulA )
00244           {
00245             nt = npos+nneg+nzero;
00246             if( nt>1 && nt<=numSec )
00247             {
00248               protmulA[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00249               protnormA[nt-1] += protmulA[counter];
00250             }
00251           }
00252         }
00253       }
00254       for( i=0; i<numMulA; ++i )neutmulA[i] = 0.0;
00255       for( i=0; i<numSec; ++i )neutnormA[i] = 0.0;
00256       counter = -1;
00257        for( npos=1; npos<numSec/3; ++npos )
00258        {
00259         nneg = npos-1;
00260         for( nzero=0; nzero<numSec/3; ++nzero )
00261         {
00262           if( ++counter < numMulA )
00263           {
00264             nt = npos+nneg+nzero;
00265             if( nt>1 && nt<=numSec )
00266             {
00267               neutmulA[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00268               neutnormA[nt-1] += neutmulA[counter];
00269             }
00270           }
00271         }
00272       }
00273       for( i=0; i<numSec; ++i )
00274       {
00275         if( protnormA[i] > 0.0 )protnormA[i] = 1.0/protnormA[i];
00276         if( neutnormA[i] > 0.0 )neutnormA[i] = 1.0/neutnormA[i];
00277       }
00278     }   // end of initialization
00279     const G4double expxu = 82.;           // upper bound for arg. of exp
00280     const G4double expxl = -expxu;        // lower bound for arg. of exp
00281     G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00282     G4ParticleDefinition *aProton = G4Proton::Proton();
00283     G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
00284     G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
00285     G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
00286     G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
00287     G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
00288     G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
00289     const G4double anhl[] = {1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,0.97,0.88,
00290                              0.85,0.81,0.75,0.64,0.64,0.55,0.55,0.45,0.47,0.40,
00291                              0.39,0.36,0.33,0.10,0.01};
00292     G4int iplab = G4int( pOriginal/GeV*10.0 );
00293     if( iplab >  9 )iplab = G4int( (pOriginal/GeV- 1.0)*5.0  ) + 10;
00294     if( iplab > 14 )iplab = G4int(  pOriginal/GeV- 2.0       ) + 15;
00295     if( iplab > 23 )iplab = G4int( (pOriginal/GeV-10.0)/10.0 ) + 23;
00296     if( iplab > 24 )iplab = 24;
00297     if( G4UniformRand() > anhl[iplab] )
00298     {
00299       if( availableEnergy <= aPiPlus->GetPDGMass()/MeV )
00300       {
00301         quasiElastic = true;
00302         return;
00303       }
00304       G4double n, anpn;
00305       GetNormalizationConstant( availableEnergy, n, anpn );
00306       G4double ran = G4UniformRand();
00307       G4double dum, excs = 0.0;
00308       if( targetParticle.GetDefinition() == aProton )
00309       {
00310         counter = -1;
00311         for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00312         {
00313           for( nneg=std::max(0,npos-2); nneg<=npos && ran>=excs; ++nneg )
00314           {
00315             for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00316             {
00317               if( ++counter < numMul )
00318               {
00319                 nt = npos+nneg+nzero;
00320                 if( nt>0 && nt<=numSec )
00321                 {
00322                   test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00323                   dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00324                   if( std::fabs(dum) < 1.0 )
00325                   {
00326                     if( test >= 1.0e-10 )excs += dum*test;
00327                   }
00328                   else
00329                     excs += dum*test;
00330                 }
00331               }
00332             }
00333           }
00334         }
00335         if( ran >= excs )  // 3 previous loops continued to the end
00336         {
00337           quasiElastic = true;
00338           return;
00339         }
00340         npos--; nneg--; nzero--;
00341         G4int ncht = std::min( 3, std::max( 1, npos-nneg+1 ) );
00342         switch( ncht )
00343         {
00344          case 1:
00345            break;
00346          case 2:
00347            if( G4UniformRand() < 0.5 )
00348            {
00349              targetParticle.SetDefinitionAndUpdateE( aNeutron );
00350              targetHasChanged = true;
00351            }
00352            else
00353            {
00354              if( G4UniformRand() < 0.5 )
00355                currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
00356              else
00357                currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
00358              incidentHasChanged = true;
00359            }             
00360            break;
00361          case 3:
00362            if( G4UniformRand() < 0.5 )
00363              currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
00364            else
00365              currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
00366            incidentHasChanged = true;
00367            targetParticle.SetDefinitionAndUpdateE( aNeutron );
00368            targetHasChanged = true;
00369            break;
00370         }
00371       }
00372       else  // target must be a neutron
00373       {
00374         counter = -1;
00375         for( npos=0; npos<numSec/3 && ran>=excs; ++npos )
00376         {
00377           for( nneg=std::max(0,npos-1); nneg<=(npos+1) && ran>=excs; ++nneg )
00378           {
00379             for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00380             {
00381               if( ++counter < numMul )
00382               {
00383                 nt = npos+nneg+nzero;
00384                 if( nt>0 && nt<=numSec )
00385                 {
00386                   test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00387                   dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00388                   if( std::fabs(dum) < 1.0 )
00389                   {
00390                     if( test >= 1.0e-10 )excs += dum*test;
00391                   }
00392                   else
00393                     excs += dum*test;
00394                 }
00395               }
00396             }
00397           }
00398         }
00399         if( ran >= excs )  // 3 previous loops continued to the end
00400         {
00401           quasiElastic = true;
00402           return;
00403         }
00404         npos--; nneg--; nzero--;
00405         G4int ncht = std::min( 3, std::max( 1, npos-nneg+2 ) );
00406         switch( ncht )
00407         {
00408          case 1:
00409            {
00410              targetParticle.SetDefinitionAndUpdateE( aProton );
00411              targetHasChanged = true;
00412            }
00413            break;
00414          case 2:
00415            if( G4UniformRand() < 0.5 )
00416            {
00417              if( G4UniformRand() < 0.5 )
00418              {
00419                currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
00420                incidentHasChanged = true;
00421                targetParticle.SetDefinitionAndUpdateE( aProton );
00422                targetHasChanged = true;
00423              }
00424            }
00425            else
00426            {
00427              if( G4UniformRand() < 0.5 )
00428              {
00429                currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
00430                incidentHasChanged = true;
00431                targetParticle.SetDefinitionAndUpdateE( aProton );
00432                targetHasChanged = true;
00433              }
00434            }
00435            break;
00436          case 3:
00437            if( G4UniformRand() < 0.5 )
00438              currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
00439            else
00440              currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
00441            incidentHasChanged = true;
00442            break;
00443         }
00444       }
00445     }
00446     else  // random number <= anhl[iplab]
00447     {
00448       if( centerofmassEnergy <= aPiPlus->GetPDGMass()/MeV+aKaonPlus->GetPDGMass()/MeV )
00449       {
00450         quasiElastic = true;
00451         return;
00452       }
00453       G4double n, anpn;
00454       GetNormalizationConstant( -centerofmassEnergy, n, anpn );
00455       G4double ran = G4UniformRand();
00456       G4double dum, excs = 0.0;
00457       if( targetParticle.GetDefinition() == aProton )
00458       {
00459         counter = -1;
00460         for( npos=2; npos<numSec/3 && ran>=excs; ++npos )
00461         {
00462           nneg=npos-2;
00463           for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00464           {
00465             if( ++counter < numMulA )
00466             {
00467               nt = npos+nneg+nzero;
00468               if( nt>1 && nt<=numSec )
00469               {
00470                 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00471                 dum = (pi/anpn)*nt*protmulA[counter]*protnormA[nt-1]/(2.0*n*n);
00472                 if( std::fabs(dum) < 1.0 )
00473                 {
00474                   if( test >= 1.0e-10 )excs += dum*test;
00475                 }
00476                 else
00477                   excs += dum*test;
00478               }
00479             }
00480           }
00481         }
00482         if( ran >= excs )  // 3 previous loops continued to the end
00483         {
00484           quasiElastic = true;
00485           return;
00486         }
00487         npos--; nzero--;
00488       }
00489       else  // target must be a neutron
00490       {
00491         counter = -1;
00492         for( npos=1; npos<numSec/3 && ran>=excs; ++npos )
00493         {
00494           nneg = npos-1;
00495           for( nzero=0; nzero<numSec/3 && ran>=excs; ++nzero )
00496           {
00497             if( ++counter < numMulA )
00498             {
00499               nt = npos+nneg+nzero;
00500               if( nt>1 && nt<=numSec )
00501               {
00502                 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00503                 dum = (pi/anpn)*nt*neutmulA[counter]*neutnormA[nt-1]/(2.0*n*n);
00504                 if( std::fabs(dum) < 1.0 )
00505                 {
00506                   if( test >= 1.0e-10 )excs += dum*test;
00507                 }
00508                 else
00509                   excs += dum*test;
00510               }
00511             }
00512           }
00513         }
00514         if( ran >= excs )  // 3 previous loops continued to the end
00515         {
00516           quasiElastic = true;
00517           return;
00518         }
00519         npos--; nzero--;
00520       }
00521       if( nzero > 0 )
00522       {
00523         if( nneg > 0 )
00524         {
00525           if( G4UniformRand() < 0.5 )
00526           {
00527             vec.Initialize( 1 );
00528             G4ReactionProduct *p = new G4ReactionProduct;
00529             p->SetDefinition( aKaonMinus );
00530             (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00531             vec.SetElement( vecLen++, p );
00532             --nneg;
00533           }
00534           else   // random number >= 0.5
00535           {
00536             vec.Initialize( 1 );
00537             G4ReactionProduct *p = new G4ReactionProduct;
00538             p->SetDefinition( aKaonZL );
00539             (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00540             vec.SetElement( vecLen++, p );
00541             --nzero;
00542           }
00543         }
00544         else   // nneg == 0
00545         {
00546           vec.Initialize( 1 );
00547           G4ReactionProduct *p = new G4ReactionProduct;
00548           p->SetDefinition( aKaonZL );
00549           (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00550           vec.SetElement( vecLen++, p );
00551           --nzero;
00552         }
00553       }
00554       else    //  nzero == 0
00555       {
00556         if( nneg > 0 )
00557         {
00558           vec.Initialize( 1 );
00559           G4ReactionProduct *p = new G4ReactionProduct;
00560           p->SetDefinition( aKaonMinus );
00561           (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00562           vec.SetElement( vecLen++, p );
00563           --nneg;
00564         }
00565       }
00566       currentParticle.SetMass( 0.0 );
00567       targetParticle.SetMass( 0.0 );
00568     }
00569     SetUpPions( npos, nneg, nzero, vec, vecLen );
00570     return;
00571 }
00572 
00573  /* end of file */
00574  

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