Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/cross_sections/src/G4EMDissociationCrossSection.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // *                                                                  *
 21 // * Parts of this code which have been  developed by QinetiQ Ltd     *
 22 // * under contract to the European Space Agency (ESA) are the        *
 23 // * intellectual property of ESA. Rights to use, copy, modify and    *
 24 // * redistribute this software for general public use are granted    *
 25 // * in compliance with any licensing, distribution and development   *
 26 // * policy adopted by the Geant4 Collaboration. This code has been   *
 27 // * written by QinetiQ Ltd for the European Space Agency, under ESA  *
 28 // * contract 17191/03/NL/LvH (Aurora Programme).                     *
 29 // *                                                                  *
 30 // * By using,  copying,  modifying or  distributing the software (or *
 31 // * any work based  on the software)  you  agree  to acknowledge its *
 32 // * use  in  resulting  scientific  publications,  and indicate your *
 33 // * acceptance of all terms of the Geant4 Software license.          *
 34 // ********************************************************************
 35 //
 36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 37 //
 38 // MODULE:    G4EMDissociationCrossSection.cc
 39 //
 40 // Version:   B.1
 41 // Date:    15/04/04
 42 // Author:    P R Truscott
 43 // Organisation:  QinetiQ Ltd, UK
 44 // Customer:    ESA/ESTEC, NOORDWIJK
 45 // Contract:    17191/03/NL/LvH
 46 //
 47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 48 //
 49 // CHANGE HISTORY
 50 // --------------
 51 //
 52 // 17 October 2003, P R Truscott, QinetiQ Ltd, UK
 53 // Created.
 54 //
 55 // 15 March 2004, P R Truscott, QinetiQ Ltd, UK
 56 // Beta release
 57 //
 58 // 30 May 2005, J.P. Wellisch removed a compilation warning on gcc 3.4 for 
 59 //               geant4 7.1.
 60 // 09 November 2010, V.Ivanchenko make class applicable for Hydrogen but 
 61 //                   set cross section for Hydrogen to zero  
 62 //
 63 // 17 August 2011, V.Ivanchenko, provide migration to new design of cross 
 64 //                 sections considering this cross section as element-wise
 65 //
 66 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 67 //////////////////////////////////////////////////////////////////////////////
 68 //
 69 #include "G4EMDissociationCrossSection.hh"
 70 #include "G4PhysicalConstants.hh"
 71 #include "G4SystemOfUnits.hh"
 72 #include "G4ParticleTable.hh"
 73 #include "G4IonTable.hh"
 74 #include "G4HadTmpUtil.hh"
 75 #include "globals.hh"
 76 #include "G4NistManager.hh"
 77 
 78 
 79 G4EMDissociationCrossSection::G4EMDissociationCrossSection ()
 80  : G4VCrossSectionDataSet("Electromagnetic dissociation")
 81 {
 82   // This function makes use of the class which can sample the virtual photon
 83   // spectrum, G4EMDissociationSpectrum.
 84 
 85   thePhotonSpectrum = new G4EMDissociationSpectrum();
 86 
 87   // Define other constants.
 88 
 89   r0      = 1.18 * fermi;
 90   J       = 36.8 * MeV;
 91   Qprime  = 17.0 * MeV;
 92   epsilon = 0.0768;
 93   xd      = 0.25;
 94 }
 95 
 96 //////////////////////////////////////////////////////////////////////////////
 97 
 98 G4EMDissociationCrossSection::~G4EMDissociationCrossSection()
 99 {
100   delete thePhotonSpectrum;
101 }
102 /////////////////////////////////////////////////////////////////////////////
103 //
104 G4bool
105 G4EMDissociationCrossSection::IsElementApplicable(const G4DynamicParticle* part,
106               G4int /*ZZ*/, const G4Material*)
107 {
108 //
109 // The condition for the applicability of this class is that the projectile
110 // must be an ion and the target must have more than one nucleon.  In reality
111 // the value of A for either the projectile or target could be much higher,
112 // since for cases where both he projectile and target are medium to small
113 // Z, the probability of the EMD process is, I think, VERY small.
114 //
115   if (G4ParticleTable::GetParticleTable()->GetIonTable()->IsIon(part->GetDefinition())) {
116     return true;
117   } else {
118     return false;
119   }
120 }
121 
122 //////////////////////////////////////////////////////////////////////////////
123 //
124 G4double G4EMDissociationCrossSection::GetElementCrossSection
125   (const G4DynamicParticle* theDynamicParticle, G4int Z,
126    const G4Material*)
127 {
128   // VI protection for Hydrogen
129   if(1 >= Z) { return 0.0; }
130 
131   // Zero cross-section for particles with kinetic energy less than 2 MeV to prevent
132   // possible abort signal from bad arithmetic in GetCrossSectionForProjectile
133   if ( theDynamicParticle->GetKineticEnergy() < 2.0*CLHEP::MeV ) { return 0.0; }
134   
135   //
136   // Get relevant information about the projectile and target (A, Z) and
137   // velocity of the projectile.
138   //
139   const G4ParticleDefinition *definitionP = theDynamicParticle->GetDefinition();
140   G4double AP   = definitionP->GetBaryonNumber();
141   G4double ZP   = definitionP->GetPDGCharge();
142   G4double b    = theDynamicParticle->GetBeta();
143   if (b <= 0.0 && b >= 1.0) { return 0.0; }
144   
145   G4double AT   = G4NistManager::Instance()->GetAtomicMassAmu(Z);
146   G4double ZT   = (G4double)Z;
147   G4double bmin = thePhotonSpectrum->GetClosestApproach(AP, ZP, AT, ZT, b);
148   //
149   //
150   // Calculate the cross-section for the projectile and then the target.  The
151   // information is returned in a G4PhysicsFreeVector, which separates out the
152   // cross-sections for the E1 and E2 moments of the virtual photon field, and
153   // the energies (GDR and GQR).
154   //
155   G4PhysicsFreeVector *theProjectileCrossSections =
156     GetCrossSectionForProjectile (AP, ZP, AT, ZT, b, bmin);
157   G4double crossSection =
158     (*theProjectileCrossSections)[0]+(*theProjectileCrossSections)[1];
159   delete theProjectileCrossSections;
160   G4PhysicsFreeVector *theTargetCrossSections =
161     GetCrossSectionForTarget (AP, ZP, AT, ZT, b, bmin);
162   crossSection +=
163     (*theTargetCrossSections)[0]+(*theTargetCrossSections)[1];
164   delete theTargetCrossSections;
165   return crossSection;
166 }
167 ////////////////////////////////////////////////////////////////////////////////
168 //
169 G4PhysicsFreeVector *
170 G4EMDissociationCrossSection::GetCrossSectionForProjectile (G4double AP,
171   G4double ZP, G4double /* AT */, G4double ZT, G4double b, G4double bmin)
172 {
173 //
174 //
175 // Use Wilson et al's approach to calculate the cross-sections due to the E1
176 // and E2 moments of the field at the giant dipole and quadrupole resonances
177 // respectively,  Note that the algorithm is traditionally applied to the
178 // EMD break-up of the projectile in the field of the target, as is implemented
179 // here.
180 //
181 // Initialise variables and calculate the energies for the GDR and GQR.
182 //
183   G4double AProot3 = G4Pow::GetInstance()->A13(AP);
184   G4double u       = 3.0 * J / Qprime / AProot3;
185   G4double R0      = r0 * AProot3;
186   G4double E_GDR  = hbarc / std::sqrt(0.7*amu_c2*R0*R0/8.0/J*
187     (1.0 + u - (1.0 + epsilon + 3.0*u)/(1.0 + epsilon + u)*epsilon));
188   G4double E_GQR  = 63.0 * MeV / AProot3;
189 //
190 //
191 // Determine the virtual photon spectra at these energies.
192 //
193   G4double ZTsq = ZT * ZT;
194   G4double nE1 = ZTsq *
195     thePhotonSpectrum->GetGeneralE1Spectrum(E_GDR, b, bmin);
196   G4double nE2 = ZTsq *
197     thePhotonSpectrum->GetGeneralE2Spectrum(E_GQR, b, bmin);
198 //
199 //
200 // Now calculate the cross-section of the projectile for interaction with the
201 // E1 and E2 fields.
202 //
203   G4double sE1 = 60.0 * millibarn * MeV * (AP-ZP)*ZP/AP;
204   G4double sE2 = 0.22 * microbarn / MeV * ZP * AProot3 * AProot3;
205   if (AP > 100.0)     sE2 *= 0.9;
206   else if (AP > 40.0) sE2 *= 0.6;
207   else                sE2 *= 0.3;
208 //
209 //
210 // ... and multiply with the intensity of the virtual photon spectra to get
211 // the probability of interaction.
212 //
213   G4PhysicsFreeVector *theCrossSectionVector = new G4PhysicsFreeVector(2);
214   theCrossSectionVector->PutValue(0, E_GDR, sE1*nE1);
215   theCrossSectionVector->PutValue(1, E_GQR, sE2*nE2*E_GQR*E_GQR);
216 
217   return theCrossSectionVector;
218 }
219 
220 ////////////////////////////////////////////////////////////////////////////////
221 //
222 G4PhysicsFreeVector *
223 G4EMDissociationCrossSection::GetCrossSectionForTarget (G4double AP,
224   G4double ZP, G4double AT, G4double ZT, G4double b, G4double bmin)
225 {
226 //
227 // This is a cheaky little member function to calculate the probability of
228 // EMD for the target in the field of the projectile ... just by reversing the
229 // A and Z's for the participants.
230 //
231   return GetCrossSectionForProjectile (AT, ZT, AP, ZP, b, bmin);
232 }
233 
234 ////////////////////////////////////////////////////////////////////////////////
235 //
236 G4double
237 G4EMDissociationCrossSection::GetWilsonProbabilityForProtonDissociation(G4double A,
238                                                                         G4double Z)
239 {
240 //
241 // This is a simple algorithm to choose whether a proton or neutron is ejected
242 // from the nucleus in the EMD interaction.
243 //
244   G4double p = 0.0;
245   if (Z < 2.0)
246     p = 0.0;  // To avoid to remove one proton from hydrogen isotopes
247   else if (Z < 6.0)
248     p = 0.5;
249   else if (Z < 8.0)
250     p = 0.6;
251   else if (Z < 14.0)
252     p = 0.7;
253   else
254   {
255     G4double p1 = (G4double) Z / (G4double) A;
256     G4double p2 = 1.95*G4Exp(-0.075*Z);
257     if (p1 < p2) p = p1;
258     else         p = p2;
259   }
260 
261   return p;
262 }
263