G4LEDeuteronInelastic.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 // Hadronic Process: Deuteron Inelastic Process
00027 // J.L. Chuma, TRIUMF, 25-Feb-1997
00028 // J.L. Chuma, 08-May-2001: Update original incident passed back in vec[0]
00029 //                          from NuclearReaction
00030 
00031 #include "G4LEDeuteronInelastic.hh"
00032 #include "G4SystemOfUnits.hh"
00033 #include "Randomize.hh"
00034 #include "G4Electron.hh"
00035 
00036 void G4LEDeuteronInelastic::ModelDescription(std::ostream& outFile) const
00037 {
00038   outFile << "G4LEDeuteronInelastic is one of the Low Energy Parameterized\n"
00039           << "(LEP) models used to implement inelastic deuteron scattering\n"
00040           << "from nuclei.  It is a re-engineered version of the GHEISHA\n"
00041           << "code of H. Fesefeldt.  It divides the initial collision\n"
00042           << "products into backward- and forward-going clusters which are\n"
00043           << "then decayed into final state hadrons.  The model does not\n"
00044           << "conserve energy on an event-by-event basis.  It may be\n"
00045           << "applied to deuterons with initial energies between 0 and 10\n"
00046           << "TeV.\n";
00047 }
00048 
00049 
00050 G4HadFinalState*
00051 G4LEDeuteronInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00052                                      G4Nucleus& targetNucleus)
00053 { 
00054   theParticleChange.Clear();
00055   const G4HadProjectile* originalIncident = &aTrack;
00056     
00057   if (verboseLevel > 1) {
00058     const G4Material *targetMaterial = aTrack.GetMaterial();
00059     G4cout << "G4LEDeuteronInelastic::ApplyYourself called" << G4endl;
00060     G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00061     G4cout << "target material = " << targetMaterial->GetName() << ", ";
00062   }
00063     
00064   // Work-around for lack of model above 100 MeV
00065   if (originalIncident->GetKineticEnergy()/MeV > 100. ||
00066       originalIncident->GetKineticEnergy() <= 0.1*MeV) {
00067     theParticleChange.SetStatusChange(isAlive);
00068     theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00069     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 
00070     return &theParticleChange;      
00071   }
00072 
00073   G4double A = targetNucleus.GetA_asInt();
00074   G4double Z = targetNucleus.GetZ_asInt();
00075   G4double theAtomicMass = targetNucleus.AtomicMass(A, Z);
00076   G4double massVec[9];
00077   massVec[0] = targetNucleus.AtomicMass( A+2.0, Z+1.0 );
00078   massVec[1] = targetNucleus.AtomicMass( A+1.0, Z+1.0 );
00079   massVec[2] = targetNucleus.AtomicMass( A+1.0, Z     );
00080   massVec[3] = theAtomicMass;
00081   massVec[4] = 0.;
00082   if (A > 1.0 && A-1.0 > Z) 
00083     massVec[4] = targetNucleus.AtomicMass(A-1.0, Z);
00084   massVec[5] = 0.;
00085   if (A > 2.0 && Z > 1.0 && A-2.0 > Z-1.0) 
00086     massVec[5] = targetNucleus.AtomicMass(A-2.0, Z-1.0);
00087   massVec[6] = 0.;
00088   if (A > Z+1.0) 
00089     massVec[6] = targetNucleus.AtomicMass(A, Z+1.0);
00090   massVec[7] = massVec[3];
00091   massVec[8] = 0.;
00092   if (Z > 1.0) massVec[8] = targetNucleus.AtomicMass(A,Z-1.0);
00093     
00094   G4FastVector<G4ReactionProduct,4> vec;  // vec will contain the secondary particles
00095   G4int vecLen = 0;
00096   vec.Initialize( 0 );
00097     
00098   theReactionDynamics.NuclearReaction(vec, vecLen, originalIncident,
00099                                       targetNucleus, theAtomicMass, massVec);
00100 
00101   G4double p = vec[0]->GetMomentum().mag();
00102   theParticleChange.SetMomentumChange( vec[0]->GetMomentum() * (1.0/p)  );
00103   theParticleChange.SetEnergyChange( vec[0]->GetKineticEnergy() );
00104   delete vec[0];
00105 
00106   if (vecLen <= 1) 
00107     {
00108       theParticleChange.SetStatusChange(isAlive);
00109       theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00110       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00111       if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00112       return &theParticleChange;      
00113     }
00114 
00115   G4DynamicParticle* pd;
00116   for (G4int i=1; i<vecLen; ++i) {
00117     pd = new G4DynamicParticle();
00118     pd->SetDefinition( vec[i]->GetDefinition() );
00119     pd->SetMomentum( vec[i]->GetMomentum() );
00120     theParticleChange.AddSecondary( pd );
00121     delete vec[i];
00122   }
00123 
00124   if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus); 
00125   return &theParticleChange;
00126 }
00127  
00128  /* end of file */
00129  

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