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