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 #include "G4DNADingfelderChargeIncreaseModel.hh"
00030 #include "G4PhysicalConstants.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "G4DNAMolecularMaterial.hh"
00033
00034
00035
00036 using namespace std;
00037
00038
00039
00040 G4DNADingfelderChargeIncreaseModel::G4DNADingfelderChargeIncreaseModel(const G4ParticleDefinition*,
00041 const G4String& nam)
00042 :G4VEmModel(nam),isInitialised(false)
00043 {
00044
00045 fpMolWaterDensity = 0;
00046
00047 numberOfPartialCrossSections[0]=0;
00048 numberOfPartialCrossSections[1]=0;
00049
00050 verboseLevel= 0;
00051
00052
00053
00054
00055
00056
00057
00058 if( verboseLevel>0 )
00059 {
00060 G4cout << "Dingfelder charge increase model is constructed " << G4endl;
00061 }
00062 fParticleChangeForGamma = 0;
00063 }
00064
00065
00066
00067 G4DNADingfelderChargeIncreaseModel::~G4DNADingfelderChargeIncreaseModel()
00068 {}
00069
00070
00071
00072 void G4DNADingfelderChargeIncreaseModel::Initialise(const G4ParticleDefinition* particle,
00073 const G4DataVector& )
00074 {
00075
00076 if (verboseLevel > 3)
00077 G4cout << "Calling G4DNADingfelderChargeIncreaseModel::Initialise()" << G4endl;
00078
00079
00080
00081 G4DNAGenericIonsManager *instance;
00082 instance = G4DNAGenericIonsManager::Instance();
00083 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
00084 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
00085 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
00086
00087 G4String hydrogen;
00088 G4String alphaPlus;
00089 G4String helium;
00090
00091
00092
00093 hydrogen = hydrogenDef->GetParticleName();
00094 lowEnergyLimit[hydrogen] = 100. * eV;
00095 highEnergyLimit[hydrogen] = 100. * MeV;
00096
00097 alphaPlus = alphaPlusDef->GetParticleName();
00098 lowEnergyLimit[alphaPlus] = 1. * keV;
00099 highEnergyLimit[alphaPlus] = 400. * MeV;
00100
00101 helium = heliumDef->GetParticleName();
00102 lowEnergyLimit[helium] = 1. * keV;
00103 highEnergyLimit[helium] = 400. * MeV;
00104
00105
00106
00107 if (particle==hydrogenDef)
00108 {
00109 SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
00110 SetHighEnergyLimit(highEnergyLimit[hydrogen]);
00111 }
00112
00113 if (particle==alphaPlusDef)
00114 {
00115 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
00116 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
00117 }
00118
00119 if (particle==heliumDef)
00120 {
00121 SetLowEnergyLimit(lowEnergyLimit[helium]);
00122 SetHighEnergyLimit(highEnergyLimit[helium]);
00123 }
00124
00125
00126
00127
00128
00129 f0[0][0]=1.;
00130 a0[0][0]=2.25;
00131 a1[0][0]=-0.75;
00132 b0[0][0]=-32.10;
00133 c0[0][0]=0.600;
00134 d0[0][0]=2.40;
00135 x0[0][0]=4.60;
00136
00137 x1[0][0]=-1.;
00138 b1[0][0]=-1.;
00139
00140 numberOfPartialCrossSections[0]=1;
00141
00142
00143
00144 f0[0][1]=1.;
00145 a0[0][1]=2.25;
00146 a1[0][1]=-0.75;
00147 b0[0][1]=-30.93;
00148 c0[0][1]=0.590;
00149 d0[0][1]=2.35;
00150 x0[0][1]=4.29;
00151
00152 f0[1][1]=1.;
00153 a0[1][1]=2.25;
00154 a1[1][1]=-0.75;
00155 b0[1][1]=-32.61;
00156 c0[1][1]=0.435;
00157 d0[1][1]=2.70;
00158 x0[1][1]=4.45;
00159
00160 x1[0][1]=-1.;
00161 b1[0][1]=-1.;
00162
00163 x1[1][1]=-1.;
00164 b1[1][1]=-1.;
00165
00166 numberOfPartialCrossSections[1]=2;
00167
00168
00169
00170 if( verboseLevel>0 )
00171 {
00172 G4cout << "Dingfelder charge increase model is initialized " << G4endl
00173 << "Energy range: "
00174 << LowEnergyLimit() / keV << " keV - "
00175 << HighEnergyLimit() / MeV << " MeV for "
00176 << particle->GetParticleName()
00177 << G4endl;
00178 }
00179
00180
00181 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00182
00183 if (isInitialised) { return; }
00184 fParticleChangeForGamma = GetParticleChangeForGamma();
00185 isInitialised = true;
00186 }
00187
00188
00189
00190 G4double G4DNADingfelderChargeIncreaseModel::CrossSectionPerVolume(const G4Material* material,
00191 const G4ParticleDefinition* particleDefinition,
00192 G4double k,
00193 G4double,
00194 G4double)
00195 {
00196 if (verboseLevel > 3)
00197 G4cout << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeIncreaseModel" << G4endl;
00198
00199
00200
00201 G4DNAGenericIonsManager *instance;
00202 instance = G4DNAGenericIonsManager::Instance();
00203
00204 if (
00205 particleDefinition != instance->GetIon("hydrogen")
00206 &&
00207 particleDefinition != instance->GetIon("alpha+")
00208 &&
00209 particleDefinition != instance->GetIon("helium")
00210 )
00211
00212 return 0;
00213
00214 G4double lowLim = 0;
00215 G4double highLim = 0;
00216 G4double totalCrossSection = 0.;
00217
00218 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
00219
00220 if(waterDensity!= 0.0)
00221
00222 {
00223 const G4String& particleName = particleDefinition->GetParticleName();
00224
00225 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
00226 pos1 = lowEnergyLimit.find(particleName);
00227
00228 if (pos1 != lowEnergyLimit.end())
00229 {
00230 lowLim = pos1->second;
00231 }
00232
00233 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
00234 pos2 = highEnergyLimit.find(particleName);
00235
00236 if (pos2 != highEnergyLimit.end())
00237 {
00238 highLim = pos2->second;
00239 }
00240
00241 if (k >= lowLim && k < highLim)
00242 {
00243
00244 if (particleDefinition == instance->GetIon("hydrogen"))
00245 {
00246 const G4double aa = 2.835;
00247 const G4double bb = 0.310;
00248 const G4double cc = 2.100;
00249 const G4double dd = 0.760;
00250 const G4double fac = 1.0e-18;
00251 const G4double rr = 13.606 * eV;
00252
00253 G4double t = k / (proton_mass_c2/electron_mass_c2);
00254 G4double x = t / rr;
00255 G4double temp = 4.0 * pi * Bohr_radius/nm * Bohr_radius/nm * fac;
00256 G4double sigmal = temp * cc * (std::pow(x,dd));
00257 G4double sigmah = temp * (aa * std::log(1.0 + x) + bb) / x;
00258 totalCrossSection = 1.0/(1.0/sigmal + 1.0/sigmah) *m*m;
00259 }
00260 else
00261 {
00262 totalCrossSection = Sum(k,particleDefinition);
00263 }
00264 }
00265
00266 if (verboseLevel > 2)
00267 {
00268 G4cout << "__________________________________" << G4endl;
00269 G4cout << "°°° G4DNADingfelderChargeIncreaseModel - XS INFO START" << G4endl;
00270 G4cout << "°°° Kinetic energy(eV)=" << k/eV << " particle : " << particleName << G4endl;
00271 G4cout << "°°° Cross section per water molecule (cm^2)=" << totalCrossSection/cm/cm << G4endl;
00272 G4cout << "°°° Cross section per water molecule (cm^-1)=" << totalCrossSection*waterDensity/(1./cm) << G4endl;
00273
00274 G4cout << "°°° G4DNADingfelderChargeIncreaseModel - XS INFO END" << G4endl;
00275 }
00276
00277 }
00278
00279 return totalCrossSection*waterDensity;
00280
00281
00282 }
00283
00284
00285
00286 void G4DNADingfelderChargeIncreaseModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00287 const G4MaterialCutsCouple* ,
00288 const G4DynamicParticle* aDynamicParticle,
00289 G4double,
00290 G4double)
00291 {
00292 if (verboseLevel > 3)
00293 G4cout << "Calling SampleSecondaries() of G4DNADingfelderChargeIncreaseModel" << G4endl;
00294
00295 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
00296 fParticleChangeForGamma->ProposeLocalEnergyDeposit(0.);
00297
00298 G4ParticleDefinition* definition = aDynamicParticle->GetDefinition();
00299
00300 G4double particleMass = definition->GetPDGMass();
00301
00302 G4double inK = aDynamicParticle->GetKineticEnergy();
00303
00304 G4int finalStateIndex = RandomSelect(inK,definition);
00305
00306 G4int n = NumberOfFinalStates(definition,finalStateIndex);
00307
00308 G4double outK = inK - IncomingParticleBindingEnergyConstant(definition,finalStateIndex);
00309
00310 G4DNAGenericIonsManager* instance;
00311 instance = G4DNAGenericIonsManager::Instance();
00312
00313 G4double electronK;
00314 if (definition == instance->GetIon("hydrogen")) electronK = inK*electron_mass_c2/proton_mass_c2;
00315 else electronK = inK*electron_mass_c2/(particleMass);
00316
00317 if (outK<0)
00318 {
00319 G4Exception("G4DNADingfelderChargeIncreaseModel::SampleSecondaries","em0004",
00320 FatalException,"Final kinetic energy is negative.");
00321 }
00322
00323 G4DynamicParticle* dp = new G4DynamicParticle
00324
00325 (OutgoingParticleDefinition(definition,finalStateIndex),
00326 aDynamicParticle->GetMomentumDirection(),
00327 outK) ;
00328
00329 fvect->push_back(dp);
00330
00331 n = n - 1;
00332
00333 while (n>0)
00334 {
00335 n--;
00336 fvect->push_back(new G4DynamicParticle
00337 (G4Electron::Electron(), aDynamicParticle->GetMomentumDirection(), electronK) );
00338 }
00339
00340 }
00341
00342
00343
00344 G4int G4DNADingfelderChargeIncreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition,
00345 G4int finalStateIndex )
00346
00347 {
00348 G4DNAGenericIonsManager* instance;
00349 instance = G4DNAGenericIonsManager::Instance();
00350 if (particleDefinition == instance->GetIon("hydrogen")) return 2;
00351 if (particleDefinition == instance->GetIon("alpha+")) return 2;
00352
00353 if (particleDefinition == instance->GetIon("helium"))
00354 { if (finalStateIndex==0) return 2;
00355 return 3;
00356 }
00357 return 0;
00358 }
00359
00360
00361
00362 G4ParticleDefinition* G4DNADingfelderChargeIncreaseModel::OutgoingParticleDefinition (G4ParticleDefinition* particleDefinition,
00363 G4int finalStateIndex)
00364 {
00365 G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
00366 if (particleDefinition == instance->GetIon("hydrogen")) return G4Proton::Proton();
00367 if (particleDefinition == instance->GetIon("alpha+")) return instance->GetIon("alpha++");
00368
00369 if (particleDefinition == instance->GetIon("helium"))
00370 {
00371 if (finalStateIndex==0) return instance->GetIon("alpha+");
00372 return instance->GetIon("alpha++");
00373 }
00374 return 0;
00375 }
00376
00377
00378
00379 G4double G4DNADingfelderChargeIncreaseModel::IncomingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
00380 G4int finalStateIndex )
00381 {
00382 G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
00383 if(particleDefinition == instance->GetIon("hydrogen")) return 13.6*eV;
00384
00385 if(particleDefinition == instance->GetIon("alpha+"))
00386 {
00387
00388
00389 return 54.509*eV;
00390 }
00391
00392 if(particleDefinition == instance->GetIon("helium"))
00393 {
00394
00395
00396
00397 if (finalStateIndex==0) return 24.587*eV;
00398 return (54.509 + 24.587)*eV;
00399 }
00400
00401 return 0;
00402 }
00403
00404
00405
00406
00407
00408 G4double G4DNADingfelderChargeIncreaseModel::PartialCrossSection(G4double k, G4int index,
00409 const G4ParticleDefinition* particleDefinition)
00410 {
00411
00412 G4int particleTypeIndex = 0;
00413 G4DNAGenericIonsManager *instance;
00414 instance = G4DNAGenericIonsManager::Instance();
00415
00416 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=0;
00417 if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=1;
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437 if (x1[index][particleTypeIndex]<x0[index][particleTypeIndex])
00438 {
00439
00440
00441
00442
00443
00444
00445
00446
00447 x1[index][particleTypeIndex]=x0[index][particleTypeIndex] + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) / (c0[index][particleTypeIndex] * d0[index][particleTypeIndex]), 1. / (d0[index][particleTypeIndex] - 1.));
00448 b1[index][particleTypeIndex]=(a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex] + b0[index][particleTypeIndex] - c0[index][particleTypeIndex] * std::pow(x1[index][particleTypeIndex] - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
00449 }
00450
00451 G4double x(std::log10(k/eV));
00452 G4double y;
00453
00454 if (x<x0[index][particleTypeIndex])
00455 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
00456 else if (x<x1[index][particleTypeIndex])
00457 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex] - c0[index][particleTypeIndex] * std::pow(x - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
00458 else
00459 y=a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
00460
00461 return f0[index][particleTypeIndex] * std::pow(10., y)*m*m;
00462
00463 }
00464
00465
00466
00467 G4int G4DNADingfelderChargeIncreaseModel::RandomSelect(G4double k,
00468 const G4ParticleDefinition* particleDefinition)
00469 {
00470
00471 G4int particleTypeIndex = 0;
00472 G4DNAGenericIonsManager *instance;
00473 instance = G4DNAGenericIonsManager::Instance();
00474
00475 if (particleDefinition == instance->GetIon("hydrogen")) return 0;
00476 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=0;
00477 if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=1;
00478
00479 const G4int n = numberOfPartialCrossSections[particleTypeIndex];
00480 G4double* values(new G4double[n]);
00481 G4double value = 0;
00482 G4int i = n;
00483
00484 while (i>0)
00485 {
00486 i--;
00487 values[i]=PartialCrossSection(k, i, particleDefinition);
00488 value+=values[i];
00489 }
00490
00491 value*=G4UniformRand();
00492
00493 i=n;
00494 while (i>0)
00495 {
00496 i--;
00497
00498 if (values[i]>value)
00499 break;
00500
00501 value-=values[i];
00502 }
00503
00504 delete[] values;
00505
00506 return i;
00507 }
00508
00509
00510
00511 G4double G4DNADingfelderChargeIncreaseModel::Sum(G4double k, const G4ParticleDefinition* particleDefinition)
00512 {
00513 G4int particleTypeIndex = 0;
00514 G4DNAGenericIonsManager *instance;
00515 instance = G4DNAGenericIonsManager::Instance();
00516
00517 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=0;
00518 if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=1;
00519
00520 G4double totalCrossSection = 0.;
00521
00522 for (G4int i=0; i<numberOfPartialCrossSections[particleTypeIndex]; i++)
00523 {
00524 totalCrossSection += PartialCrossSection(k,i,particleDefinition);
00525 }
00526 return totalCrossSection;
00527 }