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 #include "G4Nucleus.hh"
00043 #include "G4NucleiProperties.hh"
00044 #include "G4PhysicalConstants.hh"
00045 #include "G4SystemOfUnits.hh"
00046 #include "Randomize.hh"
00047 #include "G4HadronicException.hh"
00048
00049 G4Nucleus::G4Nucleus()
00050 : theA(0), theZ(0), aEff(0.0), zEff(0)
00051 {
00052 pnBlackTrackEnergy = 0.0;
00053 dtaBlackTrackEnergy = 0.0;
00054 pnBlackTrackEnergyfromAnnihilation = 0.0;
00055 dtaBlackTrackEnergyfromAnnihilation = 0.0;
00056 excitationEnergy = 0.0;
00057 momentum = G4ThreeVector(0.,0.,0.);
00058 fermiMomentum = 1.52*hbarc/fermi;
00059 theTemp = 293.16*kelvin;
00060 fIsotope = 0;
00061 }
00062
00063 G4Nucleus::G4Nucleus( const G4double A, const G4double Z )
00064 {
00065 SetParameters( A, Z );
00066 pnBlackTrackEnergy = 0.0;
00067 dtaBlackTrackEnergy = 0.0;
00068 pnBlackTrackEnergyfromAnnihilation = 0.0;
00069 dtaBlackTrackEnergyfromAnnihilation = 0.0;
00070 excitationEnergy = 0.0;
00071 momentum = G4ThreeVector(0.,0.,0.);
00072 fermiMomentum = 1.52*hbarc/fermi;
00073 theTemp = 293.16*kelvin;
00074 fIsotope = 0;
00075 }
00076
00077 G4Nucleus::G4Nucleus( const G4int A, const G4int Z )
00078 {
00079 SetParameters( A, Z );
00080 pnBlackTrackEnergy = 0.0;
00081 dtaBlackTrackEnergy = 0.0;
00082 pnBlackTrackEnergyfromAnnihilation = 0.0;
00083 dtaBlackTrackEnergyfromAnnihilation = 0.0;
00084 excitationEnergy = 0.0;
00085 momentum = G4ThreeVector(0.,0.,0.);
00086 fermiMomentum = 1.52*hbarc/fermi;
00087 theTemp = 293.16*kelvin;
00088 fIsotope = 0;
00089 }
00090
00091 G4Nucleus::G4Nucleus( const G4Material *aMaterial )
00092 {
00093 ChooseParameters( aMaterial );
00094 pnBlackTrackEnergy = 0.0;
00095 dtaBlackTrackEnergy = 0.0;
00096 pnBlackTrackEnergyfromAnnihilation = 0.0;
00097 dtaBlackTrackEnergyfromAnnihilation = 0.0;
00098 excitationEnergy = 0.0;
00099 momentum = G4ThreeVector(0.,0.,0.);
00100 fermiMomentum = 1.52*hbarc/fermi;
00101 theTemp = aMaterial->GetTemperature();
00102 fIsotope = 0;
00103 }
00104
00105 G4Nucleus::~G4Nucleus() {}
00106
00107 G4ReactionProduct G4Nucleus::
00108 GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp) const
00109 {
00110 G4double velMag = aVelocity.mag();
00111 G4ReactionProduct result;
00112 G4double value = 0;
00113 G4double random = 1;
00114 G4double norm = 3.*std::sqrt(k_Boltzmann*temp*aMass*G4Neutron::Neutron()->GetPDGMass());
00115 norm /= G4Neutron::Neutron()->GetPDGMass();
00116 norm *= 5.;
00117 norm += velMag;
00118 norm /= velMag;
00119 while(value/norm<random)
00120 {
00121 result = GetThermalNucleus(aMass, temp);
00122 G4ThreeVector targetVelocity = 1./result.GetMass()*result.GetMomentum();
00123 value = (targetVelocity+aVelocity).mag()/velMag;
00124 random = G4UniformRand();
00125 }
00126 return result;
00127 }
00128
00129 G4ReactionProduct
00130 G4Nucleus::GetThermalNucleus(G4double targetMass, G4double temp) const
00131 {
00132 G4double currentTemp = temp;
00133 if(currentTemp < 0) currentTemp = theTemp;
00134 G4ReactionProduct theTarget;
00135 theTarget.SetMass(targetMass*G4Neutron::Neutron()->GetPDGMass());
00136 G4double px, py, pz;
00137 px = GetThermalPz(theTarget.GetMass(), currentTemp);
00138 py = GetThermalPz(theTarget.GetMass(), currentTemp);
00139 pz = GetThermalPz(theTarget.GetMass(), currentTemp);
00140 theTarget.SetMomentum(px, py, pz);
00141 G4double tMom = std::sqrt(px*px+py*py+pz*pz);
00142 G4double tEtot = std::sqrt((tMom+theTarget.GetMass())*
00143 (tMom+theTarget.GetMass())-
00144 2.*tMom*theTarget.GetMass());
00145 if(1-tEtot/theTarget.GetMass()>0.001)
00146 {
00147 theTarget.SetTotalEnergy(tEtot);
00148 }
00149 else
00150 {
00151 theTarget.SetKineticEnergy(tMom*tMom/(2.*theTarget.GetMass()));
00152 }
00153 return theTarget;
00154 }
00155
00156
00157 void
00158 G4Nucleus::ChooseParameters(const G4Material* aMaterial)
00159 {
00160 G4double random = G4UniformRand();
00161 G4double sum = aMaterial->GetTotNbOfAtomsPerVolume();
00162 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
00163 G4double running(0);
00164
00165 G4Element* element = (*theElementVector)[aMaterial->GetNumberOfElements()-1];
00166
00167 for (unsigned int i = 0; i < aMaterial->GetNumberOfElements(); ++i) {
00168 running += aMaterial->GetVecNbOfAtomsPerVolume()[i];
00169 if (running > random*sum) {
00170 element = (*theElementVector)[i];
00171 break;
00172 }
00173 }
00174
00175 if (element->GetNumberOfIsotopes() > 0) {
00176 G4double randomAbundance = G4UniformRand();
00177 G4double sumAbundance = element->GetRelativeAbundanceVector()[0];
00178 unsigned int iso=0;
00179 while (iso < element->GetNumberOfIsotopes() &&
00180 sumAbundance < randomAbundance) {
00181 ++iso;
00182 sumAbundance += element->GetRelativeAbundanceVector()[iso];
00183 }
00184 theA=element->GetIsotope(iso)->GetN();
00185 theZ=element->GetIsotope(iso)->GetZ();
00186 aEff=theA;
00187 zEff=theZ;
00188 } else {
00189 aEff = element->GetN();
00190 zEff = element->GetZ();
00191 theZ = G4int(zEff + 0.5);
00192 theA = G4int(aEff + 0.5);
00193 }
00194 }
00195
00196
00197 void
00198 G4Nucleus::SetParameters(G4double A, G4double Z)
00199 {
00200 theZ = G4lrint(Z);
00201 theA = G4lrint(A);
00202 if (theA<1 || theZ<0 || theZ>theA) {
00203 throw G4HadronicException(__FILE__, __LINE__,
00204 "G4Nucleus::SetParameters called with non-physical parameters");
00205 }
00206 aEff = A;
00207 zEff = Z;
00208 fIsotope = 0;
00209 }
00210
00211 void
00212 G4Nucleus::SetParameters(G4int A, const G4int Z )
00213 {
00214 theZ = Z;
00215 theA = A;
00216 if( theA<1 || theZ<0 || theZ>theA )
00217 {
00218 throw G4HadronicException(__FILE__, __LINE__,
00219 "G4Nucleus::SetParameters called with non-physical parameters");
00220 }
00221 aEff = A;
00222 zEff = Z;
00223 fIsotope = 0;
00224 }
00225
00226 G4DynamicParticle *
00227 G4Nucleus::ReturnTargetParticle() const
00228 {
00229
00230
00231 G4DynamicParticle *targetParticle = new G4DynamicParticle;
00232 if( G4UniformRand() < zEff/aEff )
00233 targetParticle->SetDefinition( G4Proton::Proton() );
00234 else
00235 targetParticle->SetDefinition( G4Neutron::Neutron() );
00236 return targetParticle;
00237 }
00238
00239 G4double
00240 G4Nucleus::AtomicMass( const G4double A, const G4double Z ) const
00241 {
00242
00243 return G4NucleiProperties::GetNuclearMass(A, Z);
00244 }
00245
00246 G4double
00247 G4Nucleus::AtomicMass( const G4int A, const G4int Z ) const
00248 {
00249
00250 return G4NucleiProperties::GetNuclearMass(A, Z);
00251 }
00252
00253 G4double
00254 G4Nucleus::GetThermalPz( const G4double mass, const G4double temp ) const
00255 {
00256 G4double result = G4RandGauss::shoot();
00257 result *= std::sqrt(k_Boltzmann*temp*mass);
00258
00259
00260 return result;
00261 }
00262
00263 G4double
00264 G4Nucleus::EvaporationEffects( G4double kineticEnergy )
00265 {
00266
00267
00268
00269
00270
00271
00272
00273 if( aEff < 1.5 )
00274 {
00275 pnBlackTrackEnergy = dtaBlackTrackEnergy = 0.0;
00276 return 0.0;
00277 }
00278 G4double ek = kineticEnergy/GeV;
00279 G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
00280 const G4float atno = std::min( 120., aEff );
00281 const G4float gfa = 2.0*((aEff-1.0)/70.)*std::exp(-(aEff-1.0)/70.);
00282
00283
00284
00285
00286 G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*std::log(ekin) );
00287 G4float exnu = 7.716 * cfa * std::exp(-cfa)
00288 * ((atno-1.0)/120.)*std::exp(-(atno-1.0)/120.);
00289 G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin );
00290
00291
00292
00293
00294
00295
00296 pnBlackTrackEnergy = exnu*fpdiv;
00297 dtaBlackTrackEnergy = exnu*(1.0-fpdiv);
00298
00299 if( G4int(zEff+0.1) != 82 )
00300 {
00301 G4double ran1 = -6.0;
00302 G4double ran2 = -6.0;
00303 for( G4int i=0; i<12; ++i )
00304 {
00305 ran1 += G4UniformRand();
00306 ran2 += G4UniformRand();
00307 }
00308 pnBlackTrackEnergy *= 1.0 + ran1*gfa;
00309 dtaBlackTrackEnergy *= 1.0 + ran2*gfa;
00310 }
00311 pnBlackTrackEnergy = std::max( 0.0, pnBlackTrackEnergy );
00312 dtaBlackTrackEnergy = std::max( 0.0, dtaBlackTrackEnergy );
00313 while( pnBlackTrackEnergy+dtaBlackTrackEnergy >= ek )
00314 {
00315 pnBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
00316 dtaBlackTrackEnergy *= 1.0 - 0.5*G4UniformRand();
00317 }
00318
00319
00320 return (pnBlackTrackEnergy+dtaBlackTrackEnergy)*GeV;
00321 }
00322
00323 G4double G4Nucleus::AnnihilationEvaporationEffects(G4double kineticEnergy, G4double ekOrg)
00324 {
00325
00326
00327
00328 if( aEff < 1.5 || ekOrg < 0.)
00329 {
00330 pnBlackTrackEnergyfromAnnihilation = 0.0;
00331 dtaBlackTrackEnergyfromAnnihilation = 0.0;
00332 return 0.0;
00333 }
00334 G4double ek = kineticEnergy/GeV;
00335 G4float ekin = std::min( 4.0, std::max( 0.1, ek ) );
00336 const G4float atno = std::min( 120., aEff );
00337 const G4float gfa = 2.0*((aEff-1.0)/70.)*std::exp(-(aEff-1.0)/70.);
00338
00339 G4float cfa = std::max( 0.15, 0.35 + ((0.35-0.05)/2.3)*std::log(ekin) );
00340 G4float exnu = 7.716 * cfa * std::exp(-cfa)
00341 * ((atno-1.0)/120.)*std::exp(-(atno-1.0)/120.);
00342 G4float fpdiv = std::max( 0.5, 1.0-0.25*ekin*ekin );
00343
00344 pnBlackTrackEnergyfromAnnihilation = exnu*fpdiv;
00345 dtaBlackTrackEnergyfromAnnihilation = exnu*(1.0-fpdiv);
00346
00347 G4double ran1 = -6.0;
00348 G4double ran2 = -6.0;
00349 for( G4int i=0; i<12; ++i ) {
00350 ran1 += G4UniformRand();
00351 ran2 += G4UniformRand();
00352 }
00353 pnBlackTrackEnergyfromAnnihilation *= 1.0 + ran1*gfa;
00354 dtaBlackTrackEnergyfromAnnihilation *= 1.0 + ran2*gfa;
00355
00356 pnBlackTrackEnergyfromAnnihilation = std::max( 0.0, pnBlackTrackEnergyfromAnnihilation);
00357 dtaBlackTrackEnergyfromAnnihilation = std::max( 0.0, dtaBlackTrackEnergyfromAnnihilation);
00358 G4double blackSum = pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation;
00359 if (blackSum >= ekOrg/GeV) {
00360 pnBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
00361 dtaBlackTrackEnergyfromAnnihilation *= ekOrg/GeV/blackSum;
00362 }
00363
00364 return (pnBlackTrackEnergyfromAnnihilation+dtaBlackTrackEnergyfromAnnihilation)*GeV;
00365 }
00366
00367 G4double
00368 G4Nucleus::Cinema( G4double kineticEnergy )
00369 {
00370
00371
00372
00373
00374
00375 static const G4double expxu = 82.;
00376 static const G4double expxl = -expxu;
00377
00378 G4double ek = kineticEnergy/GeV;
00379 G4double ekLog = std::log( ek );
00380 G4double aLog = std::log( aEff );
00381 G4double em = std::min( 1.0, 0.2390 + 0.0408*aLog*aLog );
00382 G4double temp1 = -ek * std::min( 0.15, 0.0019*aLog*aLog*aLog );
00383 G4double temp2 = std::exp( std::max( expxl, std::min( expxu, -(ekLog-em)*(ekLog-em)*2.0 ) ) );
00384 G4double result = 0.0;
00385 if( std::abs( temp1 ) < 1.0 )
00386 {
00387 if( temp2 > 1.0e-10 )result = temp1*temp2;
00388 }
00389 else result = temp1*temp2;
00390 if( result < -ek )result = -ek;
00391 return result*GeV;
00392 }
00393
00394
00395
00396
00397
00398 G4ThreeVector G4Nucleus::GetFermiMomentum()
00399 {
00400
00401
00402
00403 G4double ranflat1=
00404 CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
00405 G4double ranflat2=
00406 CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
00407 G4double ranflat3=
00408 CLHEP::RandFlat::shoot((G4double)0.,(G4double)fermiMomentum);
00409 G4double ranmax = (ranflat1>ranflat2? ranflat1: ranflat2);
00410 ranmax = (ranmax>ranflat3? ranmax : ranflat3);
00411
00412
00413 G4double costheta = 2.*G4UniformRand() - 1.0;
00414 G4double sintheta = std::sqrt(1.0 - costheta*costheta);
00415 G4double phi = 2.0*pi*G4UniformRand();
00416
00417 G4double pz=costheta*ranmax;
00418 G4double px=sintheta*std::cos(phi)*ranmax;
00419 G4double py=sintheta*std::sin(phi)*ranmax;
00420 G4ThreeVector p(px,py,pz);
00421 return p;
00422 }
00423
00424 G4ReactionProductVector* G4Nucleus::Fragmentate()
00425 {
00426
00427 return NULL;
00428 }
00429
00430 void G4Nucleus::AddMomentum(const G4ThreeVector aMomentum)
00431 {
00432 momentum+=(aMomentum);
00433 }
00434
00435 void G4Nucleus::AddExcitationEnergy( G4double anEnergy )
00436 {
00437 excitationEnergy+=anEnergy;
00438 }
00439
00440
00441