G4KokoulinMuonNuclearXS.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 // Author:      D.H. Wright
00029 // Date:        1 February 2011
00030 //
00031 // Modified:
00032 //
00033 // 19 Aug 2011, V.Ivanchenko move to new design and make x-section per element
00034 
00035 // Description: use Kokoulin's parameterized calculation of virtual 
00036 //              photon production cross section and conversion to
00037 //              real photons.
00038 
00039 #include "G4KokoulinMuonNuclearXS.hh"
00040 
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4PhysicsTable.hh"
00044 #include "G4PhysicsLogVector.hh"
00045 
00046 
00047 G4KokoulinMuonNuclearXS::G4KokoulinMuonNuclearXS()
00048  :G4VCrossSectionDataSet("KokoulinMuonNuclearXS"), theCrossSectionTable(0),
00049   LowestKineticEnergy(1*GeV), HighestKineticEnergy(1*PeV),
00050   TotBin(60), CutFixed(0.2*GeV)
00051 {
00052   BuildCrossSectionTable();
00053 }
00054 
00055 G4KokoulinMuonNuclearXS::~G4KokoulinMuonNuclearXS()
00056 {
00057   if (theCrossSectionTable) {
00058     theCrossSectionTable->clearAndDestroy();
00059     delete theCrossSectionTable;
00060   }
00061 }
00062 
00063 G4bool 
00064 G4KokoulinMuonNuclearXS::IsElementApplicable(const G4DynamicParticle*, 
00065                                              G4int, const G4Material*)
00066 {
00067   return true;
00068 }
00069 
00070 
00071 void G4KokoulinMuonNuclearXS::BuildCrossSectionTable()
00072 {
00073   G4double lowEdgeEnergy;
00074   G4PhysicsLogVector* ptrVector;
00075 
00076   G4int nEl = G4Element::GetNumberOfElements(); 
00077   const G4ElementTable* theElementTable = G4Element::GetElementTable();
00078 
00079   theCrossSectionTable = new G4PhysicsTable(nEl);
00080 
00081   G4double AtomicNumber;
00082   G4double AtomicWeight;
00083   G4double Value;
00084 
00085   for (G4int j = 0; j < nEl; j++) {
00086     ptrVector = new G4PhysicsLogVector(LowestKineticEnergy,
00087                                        HighestKineticEnergy, TotBin);
00088     AtomicNumber = (*theElementTable)[j]->GetZ();
00089     AtomicWeight = (*theElementTable)[j]->GetA();
00090 
00091     for (G4int i = 0; i <= TotBin; ++i) {
00092       lowEdgeEnergy = ptrVector->Energy(i);
00093       Value = ComputeMicroscopicCrossSection(lowEdgeEnergy,
00094                                              AtomicNumber, AtomicWeight);
00095       ptrVector->PutValue(i,Value);
00096     }
00097 
00098     theCrossSectionTable->insertAt(j, ptrVector);
00099   }
00100 
00101   // Build (Z,element) look-up table (for use in GetZandACrossSection) 
00102   for (G4int j = 0; j < nEl; j++) {
00103     G4int ZZ = G4int((*theElementTable)[j]->GetZ());
00104     zelMap[ZZ] = j;
00105   }
00106 }
00107 
00108 G4double G4KokoulinMuonNuclearXS::
00109 ComputeMicroscopicCrossSection(G4double KineticEnergy,
00110                                G4double AtomicNumber, G4double AtomicWeight)
00111 {
00112   // Calculate cross section (differential in muon incident kinetic energy) by 
00113   // integrating the double differential cross section over the energy loss
00114 
00115   const G4double xgi[] = {0.0199,0.1017,0.2372,0.4083,0.5917,0.7628,0.8983,0.9801};
00116   const G4double wgi[] = {0.0506,0.1112,0.1569,0.1813,0.1813,0.1569,0.1112,0.0506};
00117   const G4double ak1 = 6.9;
00118   const G4double ak2 = 1.0;
00119 
00120   G4double Mass = G4MuonMinus::MuonMinus()->GetPDGMass();
00121 
00122   G4double CrossSection = 0.0;
00123   if (AtomicNumber < 1.) return CrossSection;
00124   if (KineticEnergy <= CutFixed) return CrossSection; 
00125 
00126   G4double epmin = CutFixed;
00127   G4double epmax = KineticEnergy + Mass - 0.5*proton_mass_c2;
00128   if (epmax <= epmin) return CrossSection; // NaN bug correction
00129 
00130   G4double aaa = std::log(epmin) ;
00131   G4double bbb = std::log(epmax) ;
00132   G4int kkk = G4int((bbb-aaa)/ak1 +ak2) ;
00133   G4double hhh = (bbb-aaa)/kkk ;
00134   G4double epln;
00135   G4double ep;
00136   G4double x;
00137 
00138   for (G4int l = 0; l < kkk; l++) {
00139     x = aaa + hhh*l;
00140     for (G4int ll = 0; ll < 8; ll++) {
00141       epln=x+xgi[ll]*hhh;
00142       ep = std::exp(epln);
00143       CrossSection += ep*wgi[ll]*ComputeDDMicroscopicCrossSection(KineticEnergy,  
00144                                                 AtomicNumber, AtomicWeight, ep);
00145     }
00146   }
00147 
00148   CrossSection *= hhh ;
00149   if (CrossSection < 0.) CrossSection = 0.;
00150   return CrossSection;
00151 }
00152 
00153 G4double G4KokoulinMuonNuclearXS::
00154 ComputeDDMicroscopicCrossSection(G4double KineticEnergy,
00155                                  G4double, G4double AtomicWeight,
00156                                  G4double epsilon)
00157 {
00158   // Calculate the double-differential microscopic cross section (in muon
00159   // incident kinetic energy and energy loss) using the cross section formula
00160   // of R.P. Kokoulin (18/01/98)
00161 
00162   const G4double alam2 = 0.400*GeV*GeV;
00163   const G4double alam  = 0.632456*GeV;
00164   const G4double coeffn = fine_structure_const/pi;   
00165 
00166   G4double ParticleMass = G4MuonMinus::MuonMinus()->GetPDGMass();
00167   G4double TotalEnergy = KineticEnergy + ParticleMass;
00168 
00169   G4double DCrossSection = 0.;
00170 
00171   if ((epsilon >= TotalEnergy - 0.5*proton_mass_c2) ||
00172       (epsilon <= CutFixed) ) return DCrossSection;
00173 
00174   G4double ep = epsilon/GeV;
00175   G4double a = AtomicWeight/(g/mole);
00176   G4double aeff = 0.22*a+0.78*std::exp(0.89*std::log(a));       //shadowing 
00177   G4double sigph = (49.2+11.1*std::log(ep)+151.8/std::sqrt(ep))*microbarn; 
00178   
00179   G4double v = epsilon/TotalEnergy;
00180   G4double v1 = 1.-v;
00181   G4double v2 = v*v;
00182   G4double mass2 = ParticleMass*ParticleMass;
00183 
00184   G4double up = TotalEnergy*TotalEnergy*v1/mass2*(1.+mass2*v2/(alam2*v1));
00185   G4double down = 1.+epsilon/alam*(1.+alam/(2.*proton_mass_c2)+epsilon/alam);
00186 
00187   DCrossSection = coeffn*aeff*sigph/epsilon*
00188                   (-v1+(v1+0.5*v2*(1.+2.*mass2/alam2))*std::log(up/down));
00189 
00190   if (DCrossSection < 0.) DCrossSection = 0.; 
00191   return DCrossSection;
00192 }
00193 
00194 G4double G4KokoulinMuonNuclearXS::
00195 GetElementCrossSection(const G4DynamicParticle* aPart,
00196                        G4int ZZ, const G4Material*)
00197 {
00198   G4int j = zelMap[ZZ];
00199   return (*theCrossSectionTable)[j]->Value(aPart->GetKineticEnergy());
00200 }
00201 

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