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 #include "G4DNAMillerGreenExcitationModel.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "G4DNAChemistryManager.hh"
00033 #include "G4DNAMolecularMaterial.hh"
00034
00035
00036
00037 using namespace std;
00038
00039
00040
00041 G4DNAMillerGreenExcitationModel::G4DNAMillerGreenExcitationModel(const G4ParticleDefinition*,
00042 const G4String& nam)
00043 :G4VEmModel(nam),isInitialised(false)
00044 {
00045
00046 fpMolWaterDensity = 0;
00047
00048 nLevels=0;
00049 kineticEnergyCorrection[0]=0.;
00050 kineticEnergyCorrection[1]=0.;
00051 kineticEnergyCorrection[2]=0.;
00052 kineticEnergyCorrection[3]=0.;
00053
00054 verboseLevel= 0;
00055
00056
00057
00058
00059
00060
00061
00062 if( verboseLevel>0 )
00063 {
00064 G4cout << "Miller & Green excitation model is constructed " << G4endl;
00065 }
00066 fParticleChangeForGamma = 0;
00067 }
00068
00069
00070
00071 G4DNAMillerGreenExcitationModel::~G4DNAMillerGreenExcitationModel()
00072 {}
00073
00074
00075
00076 void G4DNAMillerGreenExcitationModel::Initialise(const G4ParticleDefinition* particle,
00077 const G4DataVector& )
00078 {
00079
00080 if (verboseLevel > 3)
00081 G4cout << "Calling G4DNAMillerGreenExcitationModel::Initialise()" << G4endl;
00082
00083
00084
00085 G4DNAGenericIonsManager *instance;
00086 instance = G4DNAGenericIonsManager::Instance();
00087 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
00088 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
00089 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
00090 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
00091 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
00092
00093 G4String proton;
00094 G4String hydrogen;
00095 G4String alphaPlusPlus;
00096 G4String alphaPlus;
00097 G4String helium;
00098
00099
00100
00101 proton = protonDef->GetParticleName();
00102 lowEnergyLimit[proton] = 10. * eV;
00103 highEnergyLimit[proton] = 500. * keV;
00104
00105 kineticEnergyCorrection[0] = 1.;
00106 slaterEffectiveCharge[0][0] = 0.;
00107 slaterEffectiveCharge[1][0] = 0.;
00108 slaterEffectiveCharge[2][0] = 0.;
00109 sCoefficient[0][0] = 0.;
00110 sCoefficient[1][0] = 0.;
00111 sCoefficient[2][0] = 0.;
00112
00113 hydrogen = hydrogenDef->GetParticleName();
00114 lowEnergyLimit[hydrogen] = 10. * eV;
00115 highEnergyLimit[hydrogen] = 500. * keV;
00116
00117 kineticEnergyCorrection[0] = 1.;
00118 slaterEffectiveCharge[0][0] = 0.;
00119 slaterEffectiveCharge[1][0] = 0.;
00120 slaterEffectiveCharge[2][0] = 0.;
00121 sCoefficient[0][0] = 0.;
00122 sCoefficient[1][0] = 0.;
00123 sCoefficient[2][0] = 0.;
00124
00125 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
00126 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
00127 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
00128
00129 kineticEnergyCorrection[1] = 0.9382723/3.727417;
00130 slaterEffectiveCharge[0][1]=0.;
00131 slaterEffectiveCharge[1][1]=0.;
00132 slaterEffectiveCharge[2][1]=0.;
00133 sCoefficient[0][1]=0.;
00134 sCoefficient[1][1]=0.;
00135 sCoefficient[2][1]=0.;
00136
00137 alphaPlus = alphaPlusDef->GetParticleName();
00138 lowEnergyLimit[alphaPlus] = 1. * keV;
00139 highEnergyLimit[alphaPlus] = 400. * MeV;
00140
00141 kineticEnergyCorrection[2] = 0.9382723/3.727417;
00142 slaterEffectiveCharge[0][2]=2.0;
00143
00144
00145 slaterEffectiveCharge[1][2]=2.00;
00146 slaterEffectiveCharge[2][2]=2.00;
00147
00148 sCoefficient[0][2]=0.7;
00149 sCoefficient[1][2]=0.15;
00150 sCoefficient[2][2]=0.15;
00151
00152 helium = heliumDef->GetParticleName();
00153 lowEnergyLimit[helium] = 1. * keV;
00154 highEnergyLimit[helium] = 400. * MeV;
00155
00156 kineticEnergyCorrection[3] = 0.9382723/3.727417;
00157 slaterEffectiveCharge[0][3]=1.7;
00158 slaterEffectiveCharge[1][3]=1.15;
00159 slaterEffectiveCharge[2][3]=1.15;
00160 sCoefficient[0][3]=0.5;
00161 sCoefficient[1][3]=0.25;
00162 sCoefficient[2][3]=0.25;
00163
00164
00165
00166 if (particle==protonDef)
00167 {
00168 SetLowEnergyLimit(lowEnergyLimit[proton]);
00169 SetHighEnergyLimit(highEnergyLimit[proton]);
00170 }
00171
00172 if (particle==hydrogenDef)
00173 {
00174 SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
00175 SetHighEnergyLimit(highEnergyLimit[hydrogen]);
00176 }
00177
00178 if (particle==alphaPlusPlusDef)
00179 {
00180 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
00181 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
00182 }
00183
00184 if (particle==alphaPlusDef)
00185 {
00186 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
00187 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
00188 }
00189
00190 if (particle==heliumDef)
00191 {
00192 SetLowEnergyLimit(lowEnergyLimit[helium]);
00193 SetHighEnergyLimit(highEnergyLimit[helium]);
00194 }
00195
00196
00197
00198 nLevels = waterExcitation.NumberOfLevels();
00199
00200
00201 if( verboseLevel>0 )
00202 {
00203 G4cout << "Miller & Green excitation model is initialized " << G4endl
00204 << "Energy range: "
00205 << LowEnergyLimit() / eV << " eV - "
00206 << HighEnergyLimit() / keV << " keV for "
00207 << particle->GetParticleName()
00208 << G4endl;
00209 }
00210
00211
00212 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00213
00214 if (isInitialised) { return; }
00215 fParticleChangeForGamma = GetParticleChangeForGamma();
00216 isInitialised = true;
00217
00218 }
00219
00220
00221
00222 G4double G4DNAMillerGreenExcitationModel::CrossSectionPerVolume(const G4Material* material,
00223 const G4ParticleDefinition* particleDefinition,
00224 G4double k,
00225 G4double,
00226 G4double)
00227 {
00228 if (verboseLevel > 3)
00229 G4cout << "Calling CrossSectionPerVolume() of G4DNAMillerGreenExcitationModel" << G4endl;
00230
00231
00232
00233 G4DNAGenericIonsManager *instance;
00234 instance = G4DNAGenericIonsManager::Instance();
00235
00236 if (
00237 particleDefinition != G4Proton::ProtonDefinition()
00238 &&
00239 particleDefinition != instance->GetIon("hydrogen")
00240 &&
00241 particleDefinition != instance->GetIon("alpha++")
00242 &&
00243 particleDefinition != instance->GetIon("alpha+")
00244 &&
00245 particleDefinition != instance->GetIon("helium")
00246 )
00247
00248 return 0;
00249
00250 G4double lowLim = 0;
00251 G4double highLim = 0;
00252 G4double crossSection = 0.;
00253
00254 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
00255
00256 if(waterDensity!= 0.0)
00257
00258 {
00259
00260
00261 const G4String& particleName = particleDefinition->GetParticleName();
00262
00263 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
00264 pos1 = lowEnergyLimit.find(particleName);
00265
00266 if (pos1 != lowEnergyLimit.end())
00267 {
00268 lowLim = pos1->second;
00269 }
00270
00271 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
00272 pos2 = highEnergyLimit.find(particleName);
00273
00274 if (pos2 != highEnergyLimit.end())
00275 {
00276 highLim = pos2->second;
00277 }
00278
00279 if (k >= lowLim && k < highLim)
00280 {
00281 crossSection = Sum(k,particleDefinition);
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332 }
00333
00334 if (verboseLevel > 2)
00335 {
00336 G4cout << "__________________________________" << G4endl;
00337 G4cout << "°°° G4DNAMillerGreenExcitationModel - XS INFO START" << G4endl;
00338 G4cout << "°°° Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
00339 G4cout << "°°° Cross section per water molecule (cm^2)=" << crossSection/cm/cm << G4endl;
00340 G4cout << "°°° Cross section per water molecule (cm^-1)=" << crossSection*waterDensity/(1./cm) << G4endl;
00341
00342 G4cout << "°°° G4DNAMillerGreenExcitationModel - XS INFO END" << G4endl;
00343 }
00344 }
00345 else
00346 {
00347 if (verboseLevel > 2)
00348 {
00349 G4cout << "MillerGreenExcitationModel : WARNING Water density is NULL" << G4endl;
00350 }
00351 }
00352
00353 return crossSection*waterDensity;
00354
00355
00356 }
00357
00358
00359
00360 void G4DNAMillerGreenExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* ,
00361 const G4MaterialCutsCouple* ,
00362 const G4DynamicParticle* aDynamicParticle,
00363 G4double,
00364 G4double)
00365 {
00366
00367 if (verboseLevel > 3)
00368 G4cout << "Calling SampleSecondaries() of G4DNAMillerGreenExcitationModel" << G4endl;
00369
00370 G4double particleEnergy0 = aDynamicParticle->GetKineticEnergy();
00371
00372 G4int level = RandomSelect(particleEnergy0,aDynamicParticle->GetDefinition());
00373
00374
00375
00376
00377 const G4double excitation[]={ 8.17*eV, 10.13*eV, 11.31*eV, 12.91*eV, 14.50*eV};
00378 G4double excitationEnergy = excitation[level];
00379
00380 G4double newEnergy = particleEnergy0 - excitationEnergy;
00381
00382 if (newEnergy>0)
00383 {
00384 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
00385 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
00386 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
00387
00388 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
00389 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
00390 level,
00391 theIncomingTrack);
00392
00393 }
00394
00395 }
00396
00397
00398
00399 G4double G4DNAMillerGreenExcitationModel::PartialCrossSection(G4double k, G4int excitationLevel,
00400 const G4ParticleDefinition* particleDefinition)
00401 {
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416 const G4double sigma0(1.E+8 * barn);
00417 const G4double nu(1.);
00418 const G4double aj[]={876.*eV, 2084.* eV, 1373.*eV, 692.*eV, 900.*eV};
00419 const G4double jj[]={19820.*eV, 23490.*eV, 27770.*eV, 30830.*eV, 33080.*eV};
00420 const G4double omegaj[]={0.85, 0.88, 0.88, 0.78, 0.78};
00421
00422
00423 const G4double Eliq[5]={ 8.17*eV, 10.13*eV, 11.31*eV, 12.91*eV, 14.50*eV};
00424
00425 G4int particleTypeIndex = 0;
00426 G4DNAGenericIonsManager* instance;
00427 instance = G4DNAGenericIonsManager::Instance();
00428
00429 if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0;
00430 if (particleDefinition == instance->GetIon("hydrogen")) particleTypeIndex=0;
00431 if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex=1;
00432 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2;
00433 if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=3;
00434
00435 G4double tCorrected;
00436 tCorrected = k * kineticEnergyCorrection[particleTypeIndex];
00437
00438
00439 if (tCorrected < Eliq[excitationLevel]) return 0;
00440
00441
00442 G4int z = 10;
00443
00444 G4double numerator;
00445 numerator = std::pow(z * aj[excitationLevel], omegaj[excitationLevel]) *
00446 std::pow(tCorrected - Eliq[excitationLevel], nu);
00447
00448
00449
00450 if (particleDefinition == instance->GetIon("hydrogen"))
00451 numerator = std::pow(z * 0.75*aj[excitationLevel], omegaj[excitationLevel]) *
00452 std::pow(tCorrected - Eliq[excitationLevel], nu);
00453
00454
00455 G4double power;
00456 power = omegaj[excitationLevel] + nu;
00457
00458 G4double denominator;
00459 denominator = std::pow(jj[excitationLevel], power) + std::pow(tCorrected, power);
00460
00461 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
00462
00463 zEff -= ( sCoefficient[0][particleTypeIndex] * S_1s(k, Eliq[excitationLevel], slaterEffectiveCharge[0][particleTypeIndex], 1.) +
00464 sCoefficient[1][particleTypeIndex] * S_2s(k, Eliq[excitationLevel], slaterEffectiveCharge[1][particleTypeIndex], 2.) +
00465 sCoefficient[2][particleTypeIndex] * S_2p(k, Eliq[excitationLevel], slaterEffectiveCharge[2][particleTypeIndex], 2.) );
00466
00467 if (particleDefinition == instance->GetIon("hydrogen")) zEff = 1.;
00468
00469 G4double cross = sigma0 * zEff * zEff * numerator / denominator;
00470
00471
00472 return cross;
00473 }
00474
00475
00476
00477 G4int G4DNAMillerGreenExcitationModel::RandomSelect(G4double k,const G4ParticleDefinition* particle)
00478 {
00479 G4int i = nLevels;
00480 G4double value = 0.;
00481 std::deque<double> values;
00482
00483 G4DNAGenericIonsManager *instance;
00484 instance = G4DNAGenericIonsManager::Instance();
00485
00486 if ( particle == instance->GetIon("alpha++") ||
00487 particle == G4Proton::ProtonDefinition()||
00488 particle == instance->GetIon("hydrogen") ||
00489 particle == instance->GetIon("alpha+") ||
00490 particle == instance->GetIon("helium")
00491 )
00492 {
00493 while (i > 0)
00494 {
00495 i--;
00496 G4double partial = PartialCrossSection(k,i,particle);
00497 values.push_front(partial);
00498 value += partial;
00499 }
00500
00501 value *= G4UniformRand();
00502
00503 i = nLevels;
00504
00505 while (i > 0)
00506 {
00507 i--;
00508 if (values[i] > value) return i;
00509 value -= values[i];
00510 }
00511 }
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556 return 0;
00557 }
00558
00559
00560
00561 G4double G4DNAMillerGreenExcitationModel::Sum(G4double k, const G4ParticleDefinition* particle)
00562 {
00563 G4double totalCrossSection = 0.;
00564
00565 for (G4int i=0; i<nLevels; i++)
00566 {
00567 totalCrossSection += PartialCrossSection(k,i,particle);
00568 }
00569 return totalCrossSection;
00570 }
00571
00572
00573
00574 G4double G4DNAMillerGreenExcitationModel::S_1s(G4double t,
00575 G4double energyTransferred,
00576 G4double _slaterEffectiveCharge,
00577 G4double shellNumber)
00578 {
00579
00580
00581
00582 G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
00583 G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
00584
00585 return value;
00586 }
00587
00588
00589
00590
00591 G4double G4DNAMillerGreenExcitationModel::S_2s(G4double t,
00592 G4double energyTransferred,
00593 G4double _slaterEffectiveCharge,
00594 G4double shellNumber)
00595 {
00596
00597
00598
00599 G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
00600 G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
00601
00602 return value;
00603
00604 }
00605
00606
00607
00608 G4double G4DNAMillerGreenExcitationModel::S_2p(G4double t,
00609 G4double energyTransferred,
00610 G4double _slaterEffectiveCharge,
00611 G4double shellNumber)
00612 {
00613
00614
00615
00616 G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
00617 G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
00618
00619 return value;
00620 }
00621
00622
00623
00624 G4double G4DNAMillerGreenExcitationModel::R(G4double t,
00625 G4double energyTransferred,
00626 G4double _slaterEffectiveCharge,
00627 G4double shellNumber)
00628 {
00629
00630
00631
00632 G4double tElectron = 0.511/3728. * t;
00633
00634
00635 G4double H = 2.*13.60569172 * eV;
00636 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (_slaterEffectiveCharge/shellNumber);
00637
00638 return value;
00639 }
00640