G4LEAntiKaonZeroInelastic.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: Low Energy KaonZeroLong Inelastic Process
00029 // J.L. Chuma, TRIUMF, 11-Feb-1997
00030 // Modified by J.L.Chuma 30-Apr-97: added originalTarget for CalculateMomenta
00031  
00032 #include "G4LEAntiKaonZeroInelastic.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "Randomize.hh"
00036 #include "G4HadReentrentException.hh"
00037 
00038 
00039 void G4LEAntiKaonZeroInelastic::ModelDescription(std::ostream& outFile) const
00040 {
00041   outFile << "G4LEAntiKaonZeroInelastic is one of the Low Energy\n"
00042           << "Parameterized (LEP) models used to implement anti-K0 scattering\n"
00043           << "from nuclei.  It is a re-engineered version of the GHEISHA code\n"
00044           << "of 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           << "anti-K0s with initial energies between 0 and 25 GeV.\n";
00049 }
00050 
00051 
00052 G4HadFinalState*
00053 G4LEAntiKaonZeroInelastic::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 << "G4LEAntiKaonZeroInelastic::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 
00100   G4ReactionProduct currentParticle = modifiedOriginal;
00101   G4ReactionProduct targetParticle;
00102   targetParticle = *originalTarget;
00103   currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
00104   targetParticle.SetSide( -1 );  // target always goes in backward hemisphere
00105   G4bool incidentHasChanged = false;
00106   G4bool targetHasChanged = false;
00107   G4bool quasiElastic = false;
00108   G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;  // vec will contain the secondary particles
00109   G4int vecLen = 0;
00110   vec.Initialize( 0 );
00111     
00112   const G4double cutOff = 0.1;
00113   if (currentParticle.GetKineticEnergy()/MeV > cutOff)
00114     Cascade(vec, vecLen,
00115             originalIncident, currentParticle, targetParticle,
00116             incidentHasChanged, targetHasChanged, quasiElastic);
00117     
00118   try {
00119     CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
00120                      modifiedOriginal, targetNucleus, currentParticle,
00121                      targetParticle, incidentHasChanged, targetHasChanged,
00122                      quasiElastic);
00123   }
00124   catch(G4HadReentrentException aR)
00125   {
00126     aR.Report(G4cout);
00127     throw G4HadReentrentException(__FILE__, __LINE__, "Bailing out");
00128   }
00129   SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
00130 
00131   if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);    
00132   delete originalTarget;
00133   return &theParticleChange;
00134 }
00135 
00136  
00137 void G4LEAntiKaonZeroInelastic::Cascade(
00138    G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00139    G4int& vecLen,
00140    const G4HadProjectile *originalIncident,
00141    G4ReactionProduct &currentParticle,
00142    G4ReactionProduct &targetParticle,
00143    G4bool &incidentHasChanged,
00144    G4bool &targetHasChanged,
00145    G4bool &quasiElastic)
00146 {
00147   // derived from original FORTRAN code CASK0B by H. Fesefeldt (13-Sep-1987)
00148   //
00149   // K0Long undergoes interaction with nucleon within a nucleus.  Check if it is
00150   // energetically possible to produce pions/kaons.  In not, assume nuclear excitation
00151   // occurs and input particle is degraded in energy. No other particles are produced.
00152   // If reaction is possible, find the correct number of pions/protons/neutrons
00153   // produced using an interpolation to multiplicity data.  Replace some pions or
00154   // protons/neutrons by kaons or strange baryons according to the average
00155   // multiplicity per Inelastic reaction.
00156 
00157   const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00158   const G4double etOriginal = originalIncident->Get4Momentum().e()/MeV;
00159   const G4double pOriginal = originalIncident->Get4Momentum().vect().mag()/MeV;
00160   const G4double targetMass = targetParticle.GetMass()/MeV;
00161   G4double centerofmassEnergy = std::sqrt(mOriginal*mOriginal +
00162                                           targetMass*targetMass +
00163                                           2.0*targetMass*etOriginal );
00164   G4double availableEnergy = centerofmassEnergy - (targetMass+mOriginal);
00165 
00166   static G4bool first = true;
00167   const G4int numMul = 1200;
00168   const G4int numSec = 60;
00169   static G4double protmul[numMul], protnorm[numSec]; // proton constants
00170   static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
00171 
00172   // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
00173   G4int counter;
00174   G4int nt = 0;
00175   G4int npos = 0, nneg = 0, nzero = 0;
00176   const G4double c = 1.25;    
00177   const G4double b[] = { 0.7, 0.7 };
00178   if (first) {  // Computation of normalization constants will only be done once
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       for (nneg = std::max(0,npos - 2); nneg <= npos; ++nneg) {
00186         for (nzero = 0; nzero < numSec/3; ++nzero) {
00187           if (++counter < numMul) {
00188             nt = npos + nneg + nzero;
00189             if (nt > 0 && nt <= numSec) {
00190               protmul[counter] = Pmltpc(npos, nneg, nzero, nt, b[0], c);
00191               protnorm[nt-1] += protmul[counter];
00192             }
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 (npos = 0; npos < (numSec/3); ++npos) {
00202       for (nneg = std::max(0,npos-1); nneg <= (npos+1); ++nneg) {
00203         for (nzero = 0; nzero < numSec/3; ++nzero) {
00204           if (++counter < numMul) {
00205             nt = npos + nneg + nzero;
00206             if (nt > 0 && nt <= numSec) {
00207               neutmul[counter] = Pmltpc(npos, nneg, nzero, nt, b[1], c);
00208               neutnorm[nt-1] += neutmul[counter];
00209             }
00210           }
00211         }
00212       }
00213     }
00214 
00215     for (i = 0; i < numSec; ++i) {
00216       if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00217       if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00218     }
00219   }   // end of initialization
00220     
00221   const G4double expxu = 82.;           // upper bound for arg. of exp
00222   const G4double expxl = -expxu;        // lower bound for arg. of exp
00223   G4ParticleDefinition* aKaonMinus = G4KaonMinus::KaonMinus();
00224   G4ParticleDefinition* aKaonZS = G4KaonZeroShort::KaonZeroShort();
00225   G4ParticleDefinition* aKaonZL = G4KaonZeroLong::KaonZeroLong();
00226   G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00227   G4ParticleDefinition *aProton = G4Proton::Proton();
00228   G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
00229   G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
00230   G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
00231   G4ParticleDefinition *aLambda = G4Lambda::Lambda();
00232   G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
00233   G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
00234   G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
00235   const G4double cech[] = {1.,1.,1.,0.70,0.60,0.55,0.35,0.25,0.18,0.15};
00236   G4int iplab = G4int(std::min( 9.0, 5.0*pOriginal*MeV/GeV ));
00237   if ((pOriginal*MeV/GeV <= 2.0) && (G4UniformRand() < cech[iplab]) ) {
00238     npos = nneg = nzero = nt = 0;
00239     iplab = G4int(std::min( 19.0, pOriginal*MeV/GeV*10.0 ));
00240     const G4double cnk0[] = {0.17,0.18,0.17,0.24,0.26,0.20,0.22,0.21,0.34,0.45,
00241                              0.58,0.55,0.36,0.29,0.29,0.32,0.32,0.33,0.33,0.33};
00242     if (G4UniformRand() > cnk0[iplab]) {
00243       G4double ran = G4UniformRand();
00244       if (ran < 0.25) {
00245         // k0Long n --> pi- s+
00246         if (targetParticle.GetDefinition() == aNeutron) {
00247           currentParticle.SetDefinitionAndUpdateE(aPiMinus);
00248           targetParticle.SetDefinitionAndUpdateE(aSigmaPlus);
00249           incidentHasChanged = true;
00250           targetHasChanged = true;
00251         }
00252       } else if( ran < 0.50 ) {
00253         // k0Long p --> pi+ s0  or  k0Long n --> pi0 s0
00254           if( targetParticle.GetDefinition() == aNeutron )
00255             currentParticle.SetDefinitionAndUpdateE( aPiZero );
00256           else
00257             currentParticle.SetDefinitionAndUpdateE( aPiPlus );
00258           targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
00259           incidentHasChanged = true;
00260           targetHasChanged = true;
00261       } else if (ran < 0.75) {
00262         // k0Long n --> pi+ s-
00263           if( targetParticle.GetDefinition() == aNeutron )
00264           {
00265             currentParticle.SetDefinitionAndUpdateE( aPiPlus );
00266             targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
00267             incidentHasChanged = true;
00268             targetHasChanged = true;
00269           }
00270       } else {
00271         // k0Long p --> pi+ L  or  k0Long n --> pi0 L
00272           if( targetParticle.GetDefinition() == aNeutron )
00273             currentParticle.SetDefinitionAndUpdateE( aPiZero );
00274           else
00275             currentParticle.SetDefinitionAndUpdateE( aPiPlus );
00276           targetParticle.SetDefinitionAndUpdateE( aLambda );
00277           incidentHasChanged = true;
00278           targetHasChanged = true;
00279       }
00280     } else {
00281       // ran <= cnk0
00282       quasiElastic = true;
00283       if (targetParticle.GetDefinition() == aNeutron) {
00284           currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
00285           targetParticle.SetDefinitionAndUpdateE( aProton );
00286           incidentHasChanged = true;
00287           targetHasChanged = true;
00288       }
00289     }
00290   } else {
00291     // (pOriginal > 2.0*GeV) || (random number >= cech[iplab])
00292     if (availableEnergy < aPiPlus->GetPDGMass()/MeV) {
00293       quasiElastic = true;
00294       return;
00295     }
00296     G4double n, anpn;
00297     GetNormalizationConstant( availableEnergy, n, anpn );
00298     G4double ran = G4UniformRand();
00299     G4double dum, test, excs = 0.0;
00300     if (targetParticle.GetDefinition() == aProton) {
00301       counter = -1;
00302       for (npos = 0; (npos < numSec/3) && (ran >= excs); ++npos) {
00303         for (nneg = std::max(0,npos-2); nneg <= npos && ran >= excs; ++nneg) {
00304           for (nzero = 0; nzero < numSec/3 && ran >= excs; ++nzero) {
00305             if (++counter < numMul) {
00306               nt = npos + nneg +nzero;
00307               if (nt > 0 && nt <= numSec) {
00308                 test = std::exp(std::min(expxu,
00309                                 std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00310                 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00311                 if (std::fabs(dum) < 1.0) {
00312                   if( test >= 1.0e-10 )excs += dum*test;
00313                 } else {
00314                   excs += dum*test;
00315                 }
00316               }
00317             }
00318           }
00319         }
00320       }
00321       if (ran >= excs) {
00322         // 3 previous loops continued to the end
00323         quasiElastic = true;
00324         return;
00325       }
00326       npos--; nneg--; nzero--;
00327       switch (npos - nneg)
00328       {
00329         case 1:
00330           if (G4UniformRand() < 0.5) {
00331             currentParticle.SetDefinitionAndUpdateE(aKaonMinus);
00332             incidentHasChanged = true;
00333           } else {
00334             targetParticle.SetDefinitionAndUpdateE(aNeutron);
00335             targetHasChanged = true;
00336           }
00337         case 0:
00338           break;
00339         default:
00340           currentParticle.SetDefinitionAndUpdateE(aKaonMinus);
00341           targetParticle.SetDefinitionAndUpdateE(aNeutron);
00342           incidentHasChanged = true;
00343           targetHasChanged = true;
00344           break;
00345       }
00346     } else {
00347       // target must be a neutron
00348       counter = -1;
00349       for (npos = 0; npos < numSec/3 && ran >= excs; ++npos) {
00350         for (nneg = std::max(0,npos-1); nneg <= (npos+1) && ran >= excs; ++nneg) {
00351           for (nzero = 0; nzero < numSec/3 && ran >= excs; ++nzero) {
00352             if (++counter < numMul) {
00353               nt = npos + nneg + nzero;
00354               if (nt > 0 && nt <= numSec) {
00355                 test = std::exp(std::min(expxu,
00356                                 std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00357                 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00358                 if (std::fabs(dum) < 1.0) {
00359                   if (test >= 1.0e-10) excs += dum*test;
00360                 } else {
00361                   excs += dum*test;
00362                 }
00363               }
00364             }
00365           }
00366         }
00367       }
00368       if (ran >= excs) {
00369         // 3 previous loops continued to the end
00370         quasiElastic = true;
00371         return;
00372       }
00373       npos--; nneg--; nzero--;
00374       switch (npos - nneg)
00375       {
00376         case 0:
00377           currentParticle.SetDefinitionAndUpdateE(aKaonMinus);
00378           targetParticle.SetDefinitionAndUpdateE(aProton);
00379           incidentHasChanged = true;
00380           targetHasChanged = true;
00381           break;
00382         case 1:
00383           currentParticle.SetDefinitionAndUpdateE(aKaonMinus);
00384           incidentHasChanged = true;
00385           break;
00386         default:
00387           targetParticle.SetDefinitionAndUpdateE(aProton);
00388           targetHasChanged = true;
00389           break;
00390       }
00391     }
00392 
00393     if (G4UniformRand() >= 0.5) {
00394       if (currentParticle.GetDefinition() == aKaonMinus &&
00395           targetParticle.GetDefinition() == aNeutron) {
00396         ran = G4UniformRand();
00397         if (ran < 0.68) {
00398           currentParticle.SetDefinitionAndUpdateE(aPiMinus);
00399           targetParticle.SetDefinitionAndUpdateE(aLambda);
00400         } else if (ran < 0.84) {
00401           currentParticle.SetDefinitionAndUpdateE(aPiMinus);
00402           targetParticle.SetDefinitionAndUpdateE(aSigmaZero);
00403         } else {
00404           currentParticle.SetDefinitionAndUpdateE(aPiZero);
00405           targetParticle.SetDefinitionAndUpdateE(aSigmaMinus);
00406         }
00407       } else if ((currentParticle.GetDefinition() == aKaonZS ||
00408                   currentParticle.GetDefinition() == aKaonZL) &&
00409                   targetParticle.GetDefinition() == aProton) {
00410         ran = G4UniformRand();
00411         if (ran < 0.68) {
00412           currentParticle.SetDefinitionAndUpdateE(aPiPlus);
00413           targetParticle.SetDefinitionAndUpdateE(aLambda);
00414         } else if (ran < 0.84) {
00415           currentParticle.SetDefinitionAndUpdateE(aPiZero);
00416           targetParticle.SetDefinitionAndUpdateE(aSigmaPlus);
00417         } else {
00418           currentParticle.SetDefinitionAndUpdateE(aPiPlus);
00419           targetParticle.SetDefinitionAndUpdateE(aSigmaZero);
00420         }
00421       } else {
00422         ran = G4UniformRand();
00423         if (ran < 0.67) {
00424           currentParticle.SetDefinitionAndUpdateE(aPiZero);
00425           targetParticle.SetDefinitionAndUpdateE(aLambda);
00426         } else if (ran < 0.78) {
00427           currentParticle.SetDefinitionAndUpdateE(aPiMinus);
00428           targetParticle.SetDefinitionAndUpdateE(aSigmaPlus);
00429         } else if (ran < 0.89) {
00430           currentParticle.SetDefinitionAndUpdateE(aPiZero);
00431           targetParticle.SetDefinitionAndUpdateE(aSigmaZero);
00432         } else {
00433           currentParticle.SetDefinitionAndUpdateE(aPiPlus);
00434           targetParticle.SetDefinitionAndUpdateE(aSigmaMinus);
00435         }
00436       }
00437       incidentHasChanged = true;
00438       targetHasChanged = true;
00439     }
00440   }
00441 
00442   if (currentParticle.GetDefinition() == aKaonZL) {
00443     if (G4UniformRand() >= 0.5) {
00444       currentParticle.SetDefinitionAndUpdateE(aKaonZS);
00445       incidentHasChanged = true;
00446     }
00447   }
00448   if (targetParticle.GetDefinition() == aKaonZL) {
00449     if (G4UniformRand() >= 0.5) {
00450       targetParticle.SetDefinitionAndUpdateE(aKaonZS);
00451       targetHasChanged = true;
00452     }
00453   }
00454   SetUpPions(npos, nneg, nzero, vec, vecLen);
00455   return;
00456 }
00457 
00458  /* end of file */
00459  

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