G4LEAntiOmegaMinusInelastic.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: AntiOmegaMinus 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 // NOTE:  The FORTRAN version of the cascade, CASAOM, simply called the
00033 //        routine for the OmegaMinus particle.  Hence, the Cascade function
00034 //        below is just a copy of the Cascade from the OmegaMinus particle.
00035  
00036 #include "G4LEAntiOmegaMinusInelastic.hh"
00037 #include "G4PhysicalConstants.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "Randomize.hh"
00040 
00041 void G4LEAntiOmegaMinusInelastic::ModelDescription(std::ostream& outFile) const
00042 {
00043   outFile << "G4LEAntiOmegaMinusInelastic is one of the Low Energy\n"
00044           << "Parameterized (LEP) models used to implement inelastic\n"
00045           << "antiOmega- scattering from nuclei.  It is a re-engineered\n"
00046           << "version of the GHEISHA code of H. Fesefeldt.  It divides the\n"
00047           << "initial collision products into backward- and forward-going\n"
00048           << "clusters which are then decayed into final state hadrons.  The\n"
00049           << "model does not conserve energy on an event-by-event basis.  It\n"
00050           << "may be applied to antiOmega- with initial energies between 0\n"
00051           << "and 25 GeV.\n";
00052 }
00053 
00054 G4HadFinalState*
00055 G4LEAntiOmegaMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00056                                            G4Nucleus& targetNucleus)
00057 { 
00058   const G4HadProjectile* originalIncident = &aTrack;
00059   if (originalIncident->GetKineticEnergy()<= 0.1*MeV) {
00060     theParticleChange.SetStatusChange(isAlive);
00061     theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00062     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 
00063     return &theParticleChange;      
00064   }
00065 
00066   // create the target particle
00067   G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
00068 
00069   if (verboseLevel > 1) {
00070     const G4Material *targetMaterial = aTrack.GetMaterial();
00071     G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00072     G4cout << "target material = " << targetMaterial->GetName() << ", ";
00073     G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00074            << G4endl;
00075   }
00076 
00077   // Fermi motion and evaporation
00078   // As of Geant3, the Fermi energy calculation had not been Done
00079   G4double ek = originalIncident->GetKineticEnergy()/MeV;
00080   G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00081   G4ReactionProduct modifiedOriginal;
00082   modifiedOriginal = *originalIncident;
00083     
00084   G4double tkin = targetNucleus.Cinema( ek );
00085   ek += tkin;
00086   modifiedOriginal.SetKineticEnergy( ek*MeV );
00087   G4double et = ek + amas;
00088   G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00089   G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
00090   if (pp > 0.0) {
00091     G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00092     modifiedOriginal.SetMomentum( momentum * (p/pp) );
00093   }
00094 
00095   // calculate black track energies
00096   tkin = targetNucleus.EvaporationEffects( ek );
00097   ek -= tkin;
00098   modifiedOriginal.SetKineticEnergy( ek*MeV );
00099   et = ek + amas;
00100   p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00101   pp = modifiedOriginal.GetMomentum().mag()/MeV;
00102   if (pp > 0.0) {
00103     G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00104     modifiedOriginal.SetMomentum( momentum * (p/pp) );
00105   }
00106   G4ReactionProduct currentParticle = modifiedOriginal;
00107   G4ReactionProduct targetParticle;
00108   targetParticle = *originalTarget;
00109   currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
00110   targetParticle.SetSide( -1 );  // target always goes in backward hemisphere
00111   G4bool incidentHasChanged = false;
00112   G4bool targetHasChanged = false;
00113   G4bool quasiElastic = false;
00114   G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;  // vec will contain the secondary particles
00115   G4int vecLen = 0;
00116   vec.Initialize( 0 );
00117         
00118   const G4double cutOff = 0.1;
00119   const G4double anni = std::min( 1.3*currentParticle.GetTotalMomentum()/GeV, 0.4 );
00120     
00121   if ((currentParticle.GetKineticEnergy()/MeV > cutOff) || (G4UniformRand() > anni) )
00122     Cascade(vec, vecLen, originalIncident, currentParticle, targetParticle,
00123             incidentHasChanged, targetHasChanged, quasiElastic);
00124     
00125   CalculateMomenta(vec, vecLen, originalIncident, originalTarget,
00126                    modifiedOriginal, targetNucleus, currentParticle,
00127                    targetParticle, incidentHasChanged, targetHasChanged,
00128                    quasiElastic);
00129     
00130   SetUpChange(vec, vecLen, currentParticle, targetParticle, incidentHasChanged);
00131 
00132   if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00133 
00134   delete originalTarget;
00135   return &theParticleChange;
00136 }
00137  
00138 void G4LEAntiOmegaMinusInelastic::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 CASOM by H. Fesefeldt (31-Jan-1989)
00149   //
00150   // AntiOmegaMinus 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 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   if (availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass()/MeV) {
00166     // not energetically possible to produce pion(s)
00167     quasiElastic = true;
00168     return;
00169   }
00170   static G4bool first = true;
00171   const G4int numMul = 1200;
00172   const G4int numSec = 60;
00173   static G4double protmul[numMul], protnorm[numSec]; // proton constants
00174   static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
00175 
00176   // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
00177   G4int counter, nt=0;
00178   G4int npos = 0, nneg = 0, nzero = 0;
00179   G4double test;
00180   const G4double c = 1.25;    
00181   const G4double b[] = { 0.7, 0.7 };
00182   if (first) {  // Computation of normalization constants will only be done once
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       for (nneg = std::max(0,npos-1); nneg <= (npos+1); ++nneg) {
00190         for (nzero = 0; nzero < numSec/3; ++nzero) {
00191           if (++counter < numMul) {
00192             nt = npos+nneg+nzero;
00193             if (nt > 0 && nt <= numSec) {
00194               protmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[0],c);
00195               protnorm[nt-1] += protmul[counter];
00196             }
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       for (nneg = npos; nneg <= (npos+2); ++nneg) {
00207         for (nzero = 0; nzero < numSec/3; ++nzero) {
00208           if (++counter < numMul) {
00209             nt = npos+nneg+nzero;
00210             if ( nt>0 && nt<=numSec ) {
00211               neutmul[counter] = Pmltpc(npos,nneg,nzero,nt,b[1],c);
00212               neutnorm[nt-1] += neutmul[counter];
00213             }
00214           }
00215         }
00216       }
00217     }
00218     for (i = 0; i < numSec; ++i) {
00219       if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00220       if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00221     }
00222   }   // end of initialization
00223     
00224   const G4double expxu = 82.;           // upper bound for arg. of exp
00225   const G4double expxl = -expxu;        // lower bound for arg. of exp
00226   G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00227   G4ParticleDefinition *aProton = G4Proton::Proton();
00228   G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
00229   G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
00230   G4ParticleDefinition *aXiZero = G4XiZero::XiZero();
00231   G4double n, anpn;
00232   GetNormalizationConstant( availableEnergy, n, anpn );
00233   G4double ran = G4UniformRand();
00234   G4double dum, excs = 0.0;
00235   G4int nvefix = 0;
00236   if (targetParticle.GetDefinition() == aProton) {
00237     counter = -1;
00238     for (npos = 0; npos < numSec/3 && ran>=excs; ++npos) {
00239       for (nneg = std::max(0,npos-1); nneg <= (npos+1) && ran>=excs; ++nneg) {
00240         for (nzero = 0; nzero < numSec/3 && ran>=excs; ++nzero) {
00241           if ( ++counter < numMul ) {
00242             nt = npos+nneg+nzero;
00243             if (nt > 0 && nt <= numSec) {
00244               test = std::exp(std::min(expxu,
00245                               std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00246               dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00247               if (std::fabs(dum) < 1.0) {
00248                 if( test >= 1.0e-10 )excs += dum*test;
00249               }
00250               else
00251                 excs += dum*test;
00252             }
00253           }
00254         }
00255       }
00256     }
00257     if (ran >= excs) {
00258       // 3 previous loops continued to the end
00259       quasiElastic = true;
00260       return;
00261     }
00262     npos--; nneg--; nzero--;
00263 
00264     // number of secondary mesons determined by kno distribution
00265     // check for total charge of final state mesons to determine
00266     // the kind of baryons to be produced, taking into account
00267     // charge and strangeness conservation
00268 
00269     if (npos < nneg) {
00270       if (npos+1 == nneg) {
00271         currentParticle.SetDefinitionAndUpdateE(aXiZero);
00272         incidentHasChanged = true;
00273         nvefix = 1;
00274       } else {
00275         // charge mismatch
00276         currentParticle.SetDefinitionAndUpdateE(aSigmaPlus);
00277         incidentHasChanged = true;
00278         nvefix = 2;
00279       }
00280     } else if (npos > nneg) {
00281       targetParticle.SetDefinitionAndUpdateE(aNeutron);
00282       targetHasChanged = true;
00283     }
00284   } else {
00285     // target must be a neutron
00286     counter = -1;
00287     for (npos = 0; npos < numSec/3 && ran >= excs; ++npos) {
00288       for (nneg = npos; nneg <= (npos+2) && ran>=excs; ++nneg) {
00289         for (nzero = 0; nzero < numSec/3 && ran>=excs; ++nzero) {
00290           if (++counter < numMul) {
00291             nt = npos+nneg+nzero;
00292             if (nt > 0 && nt <= numSec) {
00293               test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00294               dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00295               if( std::fabs(dum) < 1.0 )
00296                 {
00297                   if( test >= 1.0e-10 )excs += dum*test;
00298                 }
00299                 else
00300                   excs += dum*test;
00301               }
00302           }
00303         }
00304       }
00305     }
00306     if (ran >= excs) {
00307       // 3 previous loops continued to the end
00308       quasiElastic = true;
00309       return;
00310     }
00311     npos--; nneg--; nzero--;
00312     if (npos+1 < nneg) {
00313       if( npos+2 == nneg) {
00314         currentParticle.SetDefinitionAndUpdateE( aXiZero );
00315         incidentHasChanged = true;
00316         nvefix = 1;
00317       } else {
00318         // charge mismatch
00319         currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00320         incidentHasChanged = true;
00321         nvefix = 2;
00322       }
00323       targetParticle.SetDefinitionAndUpdateE( aProton );
00324       targetHasChanged = true;
00325     } else if (npos+1 == nneg) {
00326       targetParticle.SetDefinitionAndUpdateE( aProton );
00327       targetHasChanged = true;
00328     }
00329   }
00330 
00331   SetUpPions(npos, nneg, nzero, vec, vecLen);
00332   for (G4int i = 0; i < vecLen && nvefix > 0; ++i) {
00333     if (vec[i]->GetDefinition() == G4PionMinus::PionMinus() ) {
00334 
00335       // correct the strangeness by replacing a pi- by a kaon-
00336       if (nvefix >= 1) vec[i]->SetDefinitionAndUpdateE(aKaonMinus);
00337         --nvefix;
00338     }
00339   }
00340   return;
00341 }
00342 
00343  /* end of file */
00344  

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