00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00068
00069 #include "G4EMDissociationCrossSection.hh"
00070 #include "G4PhysicalConstants.hh"
00071 #include "G4SystemOfUnits.hh"
00072 #include "G4ParticleTable.hh"
00073 #include "G4IonTable.hh"
00074 #include "G4HadTmpUtil.hh"
00075 #include "globals.hh"
00076 #include "G4NistManager.hh"
00077
00078
00079 G4EMDissociationCrossSection::G4EMDissociationCrossSection ()
00080 : G4VCrossSectionDataSet("Electromagnetic dissociation")
00081 {
00082
00083
00084
00085 thePhotonSpectrum = new G4EMDissociationSpectrum();
00086
00087
00088
00089 r0 = 1.18 * fermi;
00090 J = 36.8 * MeV;
00091 Qprime = 17.0 * MeV;
00092 epsilon = 0.0768;
00093 xd = 0.25;
00094 }
00095
00097
00098 G4EMDissociationCrossSection::~G4EMDissociationCrossSection()
00099 {
00100 delete thePhotonSpectrum;
00101 }
00103
00104 G4bool
00105 G4EMDissociationCrossSection::IsElementApplicable(const G4DynamicParticle* part,
00106 G4int , const G4Material*)
00107 {
00108
00109
00110
00111
00112
00113
00114
00115 if (G4ParticleTable::GetParticleTable()->GetIonTable()->IsIon(part->GetDefinition())) {
00116 return true;
00117 } else {
00118 return false;
00119 }
00120 }
00121
00123
00124 G4double G4EMDissociationCrossSection::GetElementCrossSection
00125 (const G4DynamicParticle* theDynamicParticle, G4int Z,
00126 const G4Material*)
00127 {
00128
00129 if(1 >= Z) { return 0.0; }
00130
00131
00132
00133
00134
00135 G4ParticleDefinition *definitionP = theDynamicParticle->GetDefinition();
00136 G4double AP = definitionP->GetBaryonNumber();
00137 G4double ZP = definitionP->GetPDGCharge();
00138 G4double b = theDynamicParticle->Get4Momentum().beta();
00139
00140 G4double AT = G4NistManager::Instance()->GetAtomicMassAmu(Z);
00141 G4double ZT = (G4double)Z;
00142 G4double bmin = thePhotonSpectrum->GetClosestApproach(AP, ZP, AT, ZT, b);
00143
00144
00145
00146
00147
00148
00149
00150 G4PhysicsFreeVector *theProjectileCrossSections =
00151 GetCrossSectionForProjectile (AP, ZP, AT, ZT, b, bmin);
00152 G4double crossSection =
00153 (*theProjectileCrossSections)[0]+(*theProjectileCrossSections)[1];
00154 delete theProjectileCrossSections;
00155 G4PhysicsFreeVector *theTargetCrossSections =
00156 GetCrossSectionForTarget (AP, ZP, AT, ZT, b, bmin);
00157 crossSection +=
00158 (*theTargetCrossSections)[0]+(*theTargetCrossSections)[1];
00159 delete theTargetCrossSections;
00160 return crossSection;
00161 }
00163
00164 G4PhysicsFreeVector *
00165 G4EMDissociationCrossSection::GetCrossSectionForProjectile (G4double AP,
00166 G4double ZP, G4double , G4double ZT, G4double b, G4double bmin)
00167 {
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178 G4double AProot3 = std::pow(AP,1.0/3.0);
00179 G4double u = 3.0 * J / Qprime / AProot3;
00180 G4double R0 = r0 * AProot3;
00181 G4double E_GDR = hbarc / std::sqrt(0.7*amu_c2*R0*R0/8.0/J*
00182 (1.0 + u - (1.0 + epsilon + 3.0*u)/(1.0 + epsilon + u)*epsilon));
00183 G4double E_GQR = 63.0 * MeV / AProot3;
00184
00185
00186
00187
00188 G4double ZTsq = ZT * ZT;
00189 G4double nE1 = ZTsq *
00190 thePhotonSpectrum->GetGeneralE1Spectrum(E_GDR, b, bmin);
00191 G4double nE2 = ZTsq *
00192 thePhotonSpectrum->GetGeneralE2Spectrum(E_GQR, b, bmin);
00193
00194
00195
00196
00197
00198 G4double sE1 = 60.0 * millibarn * MeV * (AP-ZP)*ZP/AP;
00199 G4double sE2 = 0.22 * microbarn / MeV * ZP * AProot3 * AProot3;
00200 if (AP > 100.0) sE2 *= 0.9;
00201 else if (AP > 40.0) sE2 *= 0.6;
00202 else sE2 *= 0.3;
00203
00204
00205
00206
00207
00208 G4PhysicsFreeVector *theCrossSectionVector = new G4PhysicsFreeVector(2);
00209 theCrossSectionVector->PutValue(0, E_GDR, sE1*nE1);
00210 theCrossSectionVector->PutValue(1, E_GQR, sE2*nE2*E_GQR*E_GQR);
00211
00212 return theCrossSectionVector;
00213 }
00214
00216
00217 G4PhysicsFreeVector *
00218 G4EMDissociationCrossSection::GetCrossSectionForTarget (G4double AP,
00219 G4double ZP, G4double AT, G4double ZT, G4double b, G4double bmin)
00220 {
00221
00222
00223
00224
00225
00226 return GetCrossSectionForProjectile (AT, ZT, AP, ZP, b, bmin);
00227 }
00228
00230
00231 G4double
00232 G4EMDissociationCrossSection::GetWilsonProbabilityForProtonDissociation(G4double A,
00233 G4double Z)
00234 {
00235
00236
00237
00238
00239 G4double p = 0.0;
00240 if (Z < 6.0)
00241 p = 0.5;
00242 else if (Z < 8.0)
00243 p = 0.6;
00244 else if (Z < 14.0)
00245 p = 0.7;
00246 else
00247 {
00248 G4double p1 = (G4double) Z / (G4double) A;
00249 G4double p2 = 1.95*std::exp(-0.075*Z);
00250 if (p1 < p2) p = p1;
00251 else p = p2;
00252 }
00253
00254 return p;
00255 }