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 #include "G4hImpactIonisation.hh"
00046 #include "globals.hh"
00047 #include "G4ios.hh"
00048 #include "Randomize.hh"
00049 #include "G4PhysicalConstants.hh"
00050 #include "G4SystemOfUnits.hh"
00051 #include "G4Poisson.hh"
00052 #include "G4UnitsTable.hh"
00053 #include "G4EnergyLossTables.hh"
00054 #include "G4Material.hh"
00055 #include "G4DynamicParticle.hh"
00056 #include "G4ParticleDefinition.hh"
00057 #include "G4AtomicDeexcitation.hh"
00058 #include "G4AtomicTransitionManager.hh"
00059 #include "G4PixeCrossSectionHandler.hh"
00060 #include "G4IInterpolator.hh"
00061 #include "G4LogLogInterpolator.hh"
00062 #include "G4Gamma.hh"
00063 #include "G4Electron.hh"
00064 #include "G4Proton.hh"
00065 #include "G4ProcessManager.hh"
00066 #include "G4ProductionCutsTable.hh"
00067 #include "G4PhysicsLogVector.hh"
00068 #include "G4PhysicsLinearVector.hh"
00069
00070 #include "G4VLowEnergyModel.hh"
00071 #include "G4hNuclearStoppingModel.hh"
00072 #include "G4hBetheBlochModel.hh"
00073 #include "G4hParametrisedLossModel.hh"
00074 #include "G4QAOLowEnergyLoss.hh"
00075 #include "G4hIonEffChargeSquare.hh"
00076 #include "G4IonChuFluctuationModel.hh"
00077 #include "G4IonYangFluctuationModel.hh"
00078
00079 #include "G4MaterialCutsCouple.hh"
00080 #include "G4Track.hh"
00081 #include "G4Step.hh"
00082
00083 G4hImpactIonisation::G4hImpactIonisation(const G4String& processName)
00084 : G4hRDEnergyLoss(processName),
00085 betheBlochModel(0),
00086 protonModel(0),
00087 antiprotonModel(0),
00088 theIonEffChargeModel(0),
00089 theNuclearStoppingModel(0),
00090 theIonChuFluctuationModel(0),
00091 theIonYangFluctuationModel(0),
00092 protonTable("ICRU_R49p"),
00093 antiprotonTable("ICRU_R49p"),
00094 theNuclearTable("ICRU_R49"),
00095 nStopping(true),
00096 theBarkas(true),
00097 theMeanFreePathTable(0),
00098 paramStepLimit (0.005),
00099 pixeCrossSectionHandler(0)
00100 {
00101 InitializeMe();
00102 }
00103
00104
00105
00106 void G4hImpactIonisation::InitializeMe()
00107 {
00108 LowestKineticEnergy = 10.0*eV ;
00109 HighestKineticEnergy = 100.0*GeV ;
00110 MinKineticEnergy = 10.0*eV ;
00111 TotBin = 360 ;
00112 protonLowEnergy = 1.*keV ;
00113 protonHighEnergy = 100.*MeV ;
00114 antiprotonLowEnergy = 25.*keV ;
00115 antiprotonHighEnergy = 2.*MeV ;
00116 minGammaEnergy = 100 * eV;
00117 minElectronEnergy = 250.* eV;
00118 verboseLevel = 0;
00119
00120
00121
00122 eMinPixe = 1.* keV;
00123 eMaxPixe = 200. * MeV;
00124
00125 G4String defaultPixeModel("ecpssr");
00126 modelK = defaultPixeModel;
00127 modelL = defaultPixeModel;
00128 modelM = defaultPixeModel;
00129 }
00130
00131
00132 G4hImpactIonisation::~G4hImpactIonisation()
00133 {
00134 if (theMeanFreePathTable)
00135 {
00136 theMeanFreePathTable->clearAndDestroy();
00137 delete theMeanFreePathTable;
00138 }
00139
00140 if (betheBlochModel) delete betheBlochModel;
00141 if (protonModel) delete protonModel;
00142 if (antiprotonModel) delete antiprotonModel;
00143 if (theNuclearStoppingModel) delete theNuclearStoppingModel;
00144 if (theIonEffChargeModel) delete theIonEffChargeModel;
00145 if (theIonChuFluctuationModel) delete theIonChuFluctuationModel;
00146 if (theIonYangFluctuationModel) delete theIonYangFluctuationModel;
00147
00148 delete pixeCrossSectionHandler;
00149
00150
00151
00152
00153 cutForDelta.clear();
00154 }
00155
00156
00157 void G4hImpactIonisation::SetElectronicStoppingPowerModel(const G4ParticleDefinition* particle,
00158 const G4String& dedxTable)
00159
00160 {
00161 if (particle->GetPDGCharge() > 0 )
00162 {
00163
00164 SetProtonElectronicStoppingPowerModel(dedxTable) ;
00165 }
00166 else
00167 {
00168
00169 SetAntiProtonElectronicStoppingPowerModel(dedxTable) ;
00170 }
00171 }
00172
00173
00174
00175 void G4hImpactIonisation::InitializeParametrisation()
00176
00177 {
00178
00179 betheBlochModel = new G4hBetheBlochModel("Bethe-Bloch") ;
00180 protonModel = new G4hParametrisedLossModel(protonTable) ;
00181 protonHighEnergy = std::min(protonHighEnergy,protonModel->HighEnergyLimit(0, 0));
00182 antiprotonModel = new G4QAOLowEnergyLoss(antiprotonTable) ;
00183 theNuclearStoppingModel = new G4hNuclearStoppingModel(theNuclearTable) ;
00184 theIonEffChargeModel = new G4hIonEffChargeSquare("Ziegler1988") ;
00185 theIonChuFluctuationModel = new G4IonChuFluctuationModel("Chu") ;
00186 theIonYangFluctuationModel = new G4IonYangFluctuationModel("Yang") ;
00187 }
00188
00189
00190
00191 void G4hImpactIonisation::BuildPhysicsTable(const G4ParticleDefinition& particleDef)
00192
00193
00194 {
00195
00196
00197 if(verboseLevel > 0)
00198 {
00199 G4cout << "G4hImpactIonisation::BuildPhysicsTable for "
00200 << particleDef.GetParticleName()
00201 << " mass(MeV)= " << particleDef.GetPDGMass()/MeV
00202 << " charge= " << particleDef.GetPDGCharge()/eplus
00203 << " type= " << particleDef.GetParticleType()
00204 << G4endl;
00205
00206 if(verboseLevel > 1)
00207 {
00208 G4ProcessVector* pv = particleDef.GetProcessManager()->GetProcessList();
00209
00210 G4cout << " 0: " << (*pv)[0]->GetProcessName() << " " << (*pv)[0]
00211 << " 1: " << (*pv)[1]->GetProcessName() << " " << (*pv)[1]
00212
00213 << G4endl;
00214 G4cout << "ionModel= " << theIonEffChargeModel
00215 << " MFPtable= " << theMeanFreePathTable
00216 << " iniMass= " << initialMass
00217 << G4endl;
00218 }
00219 }
00220
00221
00222 if (particleDef.GetParticleType() == "nucleus" &&
00223 particleDef.GetParticleName() != "GenericIon" &&
00224 particleDef.GetParticleSubType() == "generic")
00225 {
00226
00227 G4EnergyLossTables::Register(&particleDef,
00228 theDEDXpTable,
00229 theRangepTable,
00230 theInverseRangepTable,
00231 theLabTimepTable,
00232 theProperTimepTable,
00233 LowestKineticEnergy, HighestKineticEnergy,
00234 proton_mass_c2/particleDef.GetPDGMass(),
00235 TotBin);
00236
00237 return;
00238 }
00239
00240 if( !CutsWhereModified() && theLossTable) return;
00241
00242 InitializeParametrisation() ;
00243 G4Proton* proton = G4Proton::Proton();
00244 G4AntiProton* antiproton = G4AntiProton::AntiProton();
00245
00246 charge = particleDef.GetPDGCharge() / eplus;
00247 chargeSquare = charge*charge ;
00248
00249 const G4ProductionCutsTable* theCoupleTable=
00250 G4ProductionCutsTable::GetProductionCutsTable();
00251 size_t numOfCouples = theCoupleTable->GetTableSize();
00252
00253 cutForDelta.clear();
00254 cutForGamma.clear();
00255
00256 for (size_t j=0; j<numOfCouples; j++) {
00257
00258
00259 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
00260 const G4Material* material= couple->GetMaterial();
00261
00262
00263 G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[j];
00264 if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
00265
00266 G4double excEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
00267
00268 tCut = std::max(tCut,excEnergy);
00269 cutForDelta.push_back(tCut);
00270
00271
00272 tCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[j];
00273 if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
00274 tCut = std::max(tCut,minGammaEnergy);
00275 cutForGamma.push_back(tCut);
00276 }
00277
00278 if(verboseLevel > 0) {
00279 G4cout << "Cuts are defined " << G4endl;
00280 }
00281
00282 if(0.0 < charge)
00283 {
00284 {
00285 BuildLossTable(*proton) ;
00286
00287
00288
00289
00290
00291 RecorderOfpProcess[CounterOfpProcess] = theLossTable ;
00292 CounterOfpProcess++;
00293 }
00294 } else {
00295 {
00296 BuildLossTable(*antiproton) ;
00297
00298
00299
00300
00301
00302 RecorderOfpbarProcess[CounterOfpbarProcess] = theLossTable ;
00303 CounterOfpbarProcess++;
00304 }
00305 }
00306
00307 if(verboseLevel > 0) {
00308 G4cout << "G4hImpactIonisation::BuildPhysicsTable: "
00309 << "Loss table is built "
00310
00311 << G4endl;
00312 }
00313
00314 BuildLambdaTable(particleDef) ;
00315
00316
00317 if(verboseLevel > 1) {
00318 G4cout << (*theMeanFreePathTable) << G4endl;
00319 }
00320
00321 if(verboseLevel > 0) {
00322 G4cout << "G4hImpactIonisation::BuildPhysicsTable: "
00323 << "DEDX table will be built "
00324
00325
00326 << G4endl;
00327 }
00328
00329 BuildDEDXTable(particleDef) ;
00330
00331 if(verboseLevel > 1) {
00332 G4cout << (*theDEDXpTable) << G4endl;
00333 }
00334
00335 if((&particleDef == proton) || (&particleDef == antiproton)) PrintInfoDefinition() ;
00336
00337 if(verboseLevel > 0) {
00338 G4cout << "G4hImpactIonisation::BuildPhysicsTable: end for "
00339 << particleDef.GetParticleName() << G4endl;
00340 }
00341 }
00342
00343
00344
00345
00346
00347
00348 void G4hImpactIonisation::BuildLossTable(const G4ParticleDefinition& particleDef)
00349 {
00350
00351 G4double lowEdgeEnergy , ionloss, ionlossBB, paramB ;
00352
00353 G4double highEnergy;
00354 G4Proton* proton = G4Proton::Proton();
00355
00356 if (particleDef == *proton)
00357 {
00358
00359 highEnergy = protonHighEnergy ;
00360 charge = 1. ;
00361 }
00362 else
00363 {
00364
00365 highEnergy = antiprotonHighEnergy ;
00366 charge = -1. ;
00367 }
00368 chargeSquare = 1. ;
00369
00370 const G4ProductionCutsTable* theCoupleTable=
00371 G4ProductionCutsTable::GetProductionCutsTable();
00372 size_t numOfCouples = theCoupleTable->GetTableSize();
00373
00374 if ( theLossTable)
00375 {
00376 theLossTable->clearAndDestroy();
00377 delete theLossTable;
00378 }
00379
00380 theLossTable = new G4PhysicsTable(numOfCouples);
00381
00382
00383 for (size_t j=0; j<numOfCouples; j++) {
00384
00385
00386 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(LowestKineticEnergy,
00387 HighestKineticEnergy,
00388 TotBin);
00389
00390
00391 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
00392 const G4Material* material= couple->GetMaterial();
00393
00394 if ( charge > 0.0 ) {
00395 ionloss = ProtonParametrisedDEDX(couple,highEnergy) ;
00396 } else {
00397 ionloss = AntiProtonParametrisedDEDX(couple,highEnergy) ;
00398 }
00399
00400 ionlossBB = betheBlochModel->TheValue(&particleDef,material,highEnergy) ;
00401 ionlossBB -= DeltaRaysEnergy(couple,highEnergy,proton_mass_c2) ;
00402
00403
00404 paramB = ionloss/ionlossBB - 1.0 ;
00405
00406
00407 for (G4int i = 0 ; i < TotBin ; i++) {
00408 lowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
00409
00410
00411 if ( lowEdgeEnergy < highEnergy ) {
00412
00413 if ( charge > 0.0 ) {
00414 ionloss = ProtonParametrisedDEDX(couple,lowEdgeEnergy) ;
00415 } else {
00416 ionloss = AntiProtonParametrisedDEDX(couple,lowEdgeEnergy) ;
00417 }
00418
00419 } else {
00420
00421
00422 ionloss = betheBlochModel->TheValue(proton,material,
00423 lowEdgeEnergy) ;
00424
00425 ionloss -= DeltaRaysEnergy(couple,lowEdgeEnergy,proton_mass_c2) ;
00426
00427 ionloss *= (1.0 + paramB*highEnergy/lowEdgeEnergy) ;
00428 }
00429
00430
00431 if(verboseLevel > 1) {
00432 G4cout << "E(MeV)= " << lowEdgeEnergy/MeV
00433 << " dE/dx(MeV/mm)= " << ionloss*mm/MeV
00434 << " in " << material->GetName() << G4endl;
00435 }
00436 aVector->PutValue(i,ionloss) ;
00437 }
00438
00439 theLossTable->insert(aVector) ;
00440 }
00441 }
00442
00443
00444
00445
00446 void G4hImpactIonisation::BuildLambdaTable(const G4ParticleDefinition& particleDef)
00447
00448 {
00449
00450
00451
00452 if(verboseLevel > 1) {
00453 G4cout << "G4hImpactIonisation::BuildLambdaTable for "
00454 << particleDef.GetParticleName() << " is started" << G4endl;
00455 }
00456
00457
00458 G4double lowEdgeEnergy, value;
00459 charge = particleDef.GetPDGCharge()/eplus ;
00460 chargeSquare = charge*charge ;
00461 initialMass = particleDef.GetPDGMass();
00462
00463 const G4ProductionCutsTable* theCoupleTable=
00464 G4ProductionCutsTable::GetProductionCutsTable();
00465 size_t numOfCouples = theCoupleTable->GetTableSize();
00466
00467
00468 if (theMeanFreePathTable) {
00469 theMeanFreePathTable->clearAndDestroy();
00470 delete theMeanFreePathTable;
00471 }
00472
00473 theMeanFreePathTable = new G4PhysicsTable(numOfCouples);
00474
00475
00476
00477 for (size_t j=0 ; j < numOfCouples; j++) {
00478
00479
00480 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(LowestKineticEnergy,
00481 HighestKineticEnergy,
00482 TotBin);
00483
00484
00485 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
00486 const G4Material* material= couple->GetMaterial();
00487
00488 const G4ElementVector* theElementVector = material->GetElementVector() ;
00489 const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
00490 const G4int numberOfElements = material->GetNumberOfElements() ;
00491
00492
00493
00494
00495
00496
00497 G4double deltaCut = cutForDelta[j];
00498
00499 for ( G4int i = 0 ; i < TotBin ; i++ ) {
00500 lowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
00501 G4double sigma = 0.0 ;
00502 G4int Z;
00503
00504 for (G4int iel=0; iel<numberOfElements; iel++ )
00505 {
00506 Z = (G4int) (*theElementVector)[iel]->GetZ();
00507
00508
00509 G4double microCross = MicroscopicCrossSection( particleDef,
00510 lowEdgeEnergy,
00511 Z,
00512 deltaCut ) ;
00513
00514 sigma += theAtomicNumDensityVector[iel] * microCross ;
00515 }
00516
00517
00518
00519 value = sigma<=0 ? DBL_MAX : 1./sigma ;
00520
00521 aVector->PutValue(i, value) ;
00522 }
00523
00524 theMeanFreePathTable->insert(aVector);
00525 }
00526
00527 }
00528
00529
00530
00531 G4double G4hImpactIonisation::MicroscopicCrossSection(const G4ParticleDefinition& particleDef,
00532 G4double kineticEnergy,
00533 G4double atomicNumber,
00534 G4double deltaCutInEnergy) const
00535 {
00536
00537
00538
00539
00540
00541
00542
00543
00544 G4double totalCrossSection = 0.;
00545
00546
00547 G4double particleMass = initialMass;
00548 G4double energy = kineticEnergy + particleMass;
00549
00550
00551 G4double gamma = energy / particleMass;
00552 G4double beta2 = 1. - 1. / (gamma * gamma);
00553 G4double var = electron_mass_c2 / particleMass;
00554 G4double tMax = 2. * electron_mass_c2 * (gamma*gamma - 1.) / (1. + 2.* gamma*var + var*var);
00555
00556
00557
00558 if ( tMax > deltaCutInEnergy )
00559 {
00560 var = deltaCutInEnergy / tMax;
00561 totalCrossSection = (1. - var * (1. - beta2 * std::log(var))) / deltaCutInEnergy ;
00562
00563 G4double spin = particleDef.GetPDGSpin() ;
00564
00565
00566 if (spin == 0.5)
00567 {
00568 totalCrossSection += 0.5 * (tMax - deltaCutInEnergy) / (energy*energy);
00569 }
00570
00571 else if (spin > 0.9 )
00572 {
00573 totalCrossSection += -std::log(var) /
00574 (3. * deltaCutInEnergy) + (tMax - deltaCutInEnergy) * ( (5. + 1. /var)*0.25 / (energy*energy) -
00575 beta2 / (tMax * deltaCutInEnergy) ) / 3. ;
00576 }
00577 totalCrossSection *= twopi_mc2_rcl2 * atomicNumber / beta2 ;
00578 }
00579
00580
00581
00582
00583 return totalCrossSection ;
00584 }
00585
00586
00587
00588
00589 G4double G4hImpactIonisation::GetMeanFreePath(const G4Track& track,
00590 G4double,
00591 enum G4ForceCondition* condition)
00592 {
00593 const G4DynamicParticle* dynamicParticle = track.GetDynamicParticle();
00594 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
00595 const G4Material* material = couple->GetMaterial();
00596
00597 G4double meanFreePath = DBL_MAX;
00598
00599 G4bool isOutRange = false;
00600
00601 *condition = NotForced ;
00602
00603 G4double kineticEnergy = (dynamicParticle->GetKineticEnergy())*initialMass/(dynamicParticle->GetMass());
00604 charge = dynamicParticle->GetCharge()/eplus;
00605 chargeSquare = theIonEffChargeModel->TheValue(dynamicParticle, material);
00606
00607 if (kineticEnergy < LowestKineticEnergy)
00608 {
00609 meanFreePath = DBL_MAX;
00610 }
00611 else
00612 {
00613 if (kineticEnergy > HighestKineticEnergy) kineticEnergy = HighestKineticEnergy;
00614 meanFreePath = (((*theMeanFreePathTable)(couple->GetIndex()))->
00615 GetValue(kineticEnergy,isOutRange))/chargeSquare;
00616 }
00617
00618 return meanFreePath ;
00619 }
00620
00621
00622
00623 G4double G4hImpactIonisation::GetConstraints(const G4DynamicParticle* particle,
00624 const G4MaterialCutsCouple* couple)
00625 {
00626
00627
00628
00629
00630 const G4Material* material = couple->GetMaterial();
00631 G4Proton* proton = G4Proton::Proton();
00632 G4AntiProton* antiproton = G4AntiProton::AntiProton();
00633
00634 G4double stepLimit = 0.;
00635 G4double dx, highEnergy;
00636
00637 G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
00638 G4double kineticEnergy = particle->GetKineticEnergy() ;
00639
00640
00641
00642 G4double tScaled = kineticEnergy*massRatio ;
00643 fBarkas = 0.;
00644
00645 if (charge > 0.)
00646 {
00647 highEnergy = protonHighEnergy ;
00648 fRangeNow = G4EnergyLossTables::GetRange(proton, tScaled, couple);
00649 dx = G4EnergyLossTables::GetRange(proton, highEnergy, couple);
00650 fdEdx = G4EnergyLossTables::GetDEDX(proton, tScaled, couple)
00651 * chargeSquare ;
00652
00653
00654 if (theBarkas && tScaled > highEnergy)
00655 {
00656 fBarkas = BarkasTerm(material,tScaled)*std::sqrt(chargeSquare)*chargeSquare
00657 + BlochTerm(material,tScaled,chargeSquare);
00658 }
00659
00660 }
00661 else
00662 {
00663 highEnergy = antiprotonHighEnergy ;
00664 fRangeNow = G4EnergyLossTables::GetRange(antiproton, tScaled, couple);
00665 dx = G4EnergyLossTables::GetRange(antiproton, highEnergy, couple);
00666 fdEdx = G4EnergyLossTables::GetDEDX(antiproton, tScaled, couple) * chargeSquare ;
00667
00668 if (theBarkas && tScaled > highEnergy)
00669 {
00670 fBarkas = -BarkasTerm(material,tScaled)*std::sqrt(chargeSquare)*chargeSquare
00671 + BlochTerm(material,tScaled,chargeSquare);
00672 }
00673 }
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686 fRangeNow /= (chargeSquare*massRatio) ;
00687 dx /= (chargeSquare*massRatio) ;
00688
00689 stepLimit = fRangeNow ;
00690 G4double r = std::min(finalRange, couple->GetProductionCuts()
00691 ->GetProductionCut(idxG4ElectronCut));
00692
00693 if (fRangeNow > r)
00694 {
00695 stepLimit = dRoverRange*fRangeNow + r*(1.0 - dRoverRange)*(2.0 - r/fRangeNow);
00696 if (stepLimit > fRangeNow) stepLimit = fRangeNow;
00697 }
00698
00699 if(tScaled > highEnergy )
00700 {
00701
00702 fdEdx += fBarkas;
00703
00704 if(stepLimit > fRangeNow - dx*0.9) stepLimit = fRangeNow - dx*0.9 ;
00705
00706
00707 }
00708 else
00709 {
00710 G4double x = dx*paramStepLimit;
00711 if (stepLimit > x) stepLimit = x;
00712 }
00713 return stepLimit;
00714 }
00715
00716
00717
00718 G4VParticleChange* G4hImpactIonisation::AlongStepDoIt(const G4Track& track,
00719 const G4Step& step)
00720 {
00721
00722 G4Proton* proton = G4Proton::Proton();
00723 G4AntiProton* antiproton = G4AntiProton::AntiProton();
00724 G4double finalT = 0.;
00725
00726 aParticleChange.Initialize(track) ;
00727
00728 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
00729 const G4Material* material = couple->GetMaterial();
00730
00731
00732 const G4double stepLength = step.GetStepLength() ;
00733
00734 const G4DynamicParticle* particle = track.GetDynamicParticle() ;
00735
00736 G4double kineticEnergy = particle->GetKineticEnergy() ;
00737 G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
00738 G4double tScaled = kineticEnergy * massRatio ;
00739 G4double eLoss = 0.0 ;
00740 G4double nLoss = 0.0 ;
00741
00742
00743
00744 if (kineticEnergy < MinKineticEnergy)
00745 {
00746 eLoss = kineticEnergy ;
00747
00748 }
00749
00750 else if( kineticEnergy > HighestKineticEnergy)
00751 {
00752 eLoss = stepLength * fdEdx ;
00753
00754 }
00755 else if (stepLength >= fRangeNow )
00756 {
00757 eLoss = kineticEnergy ;
00758
00759
00760 }
00761 else
00762 {
00763
00764 if(stepLength > linLossLimit * fRangeNow)
00765 {
00766 G4double rScaled = fRangeNow * massRatio * chargeSquare ;
00767 G4double sScaled = stepLength * massRatio * chargeSquare ;
00768
00769 if(charge > 0.0)
00770 {
00771 eLoss = G4EnergyLossTables::GetPreciseEnergyFromRange(proton,rScaled, couple) -
00772 G4EnergyLossTables::GetPreciseEnergyFromRange(proton,rScaled-sScaled,couple) ;
00773
00774 }
00775 else
00776 {
00777
00778 eLoss = G4EnergyLossTables::GetPreciseEnergyFromRange(antiproton,rScaled,couple) -
00779 G4EnergyLossTables::GetPreciseEnergyFromRange(antiproton,rScaled-sScaled,couple) ;
00780 }
00781 eLoss /= massRatio ;
00782
00783
00784 eLoss += fBarkas * stepLength;
00785
00786
00787 }
00788 else
00789 {
00790 eLoss = stepLength *fdEdx ;
00791 }
00792 if (nStopping && tScaled < protonHighEnergy)
00793 {
00794 nLoss = (theNuclearStoppingModel->TheValue(particle, material)) * stepLength;
00795 }
00796 }
00797
00798 if (eLoss < 0.0) eLoss = 0.0;
00799
00800 finalT = kineticEnergy - eLoss - nLoss;
00801
00802 if ( EnlossFlucFlag && 0.0 < eLoss && finalT > MinKineticEnergy)
00803 {
00804
00805
00806 eLoss = ElectronicLossFluctuation(particle, couple, eLoss, stepLength) ;
00807 if (eLoss < 0.0) eLoss = 0.0;
00808 finalT = kineticEnergy - eLoss - nLoss;
00809 }
00810
00811
00812 if (finalT*massRatio <= MinKineticEnergy )
00813 {
00814
00815 finalT = 0.0;
00816 if (!particle->GetDefinition()->GetProcessManager()->GetAtRestProcessVector()->size())
00817 aParticleChange.ProposeTrackStatus(fStopAndKill);
00818 else
00819 aParticleChange.ProposeTrackStatus(fStopButAlive);
00820 }
00821
00822 aParticleChange.ProposeEnergy( finalT );
00823 eLoss = kineticEnergy-finalT;
00824
00825 aParticleChange.ProposeLocalEnergyDeposit(eLoss);
00826 return &aParticleChange ;
00827 }
00828
00829
00830
00831
00832 G4double G4hImpactIonisation::ProtonParametrisedDEDX(const G4MaterialCutsCouple* couple,
00833 G4double kineticEnergy) const
00834 {
00835 const G4Material* material = couple->GetMaterial();
00836 G4Proton* proton = G4Proton::Proton();
00837 G4double eLoss = 0.;
00838
00839
00840 if(kineticEnergy < protonLowEnergy) {
00841 eLoss = (protonModel->TheValue(proton, material, protonLowEnergy))
00842 * std::sqrt(kineticEnergy/protonLowEnergy) ;
00843
00844
00845 } else {
00846 eLoss = protonModel->TheValue(proton, material, kineticEnergy) ;
00847 }
00848
00849
00850 eLoss -= DeltaRaysEnergy(couple,kineticEnergy,proton_mass_c2) ;
00851
00852 if(verboseLevel > 2) {
00853 G4cout << "p E(MeV)= " << kineticEnergy/MeV
00854 << " dE/dx(MeV/mm)= " << eLoss*mm/MeV
00855 << " for " << material->GetName()
00856 << " model: " << protonModel << G4endl;
00857 }
00858
00859 if(eLoss < 0.0) eLoss = 0.0 ;
00860
00861 return eLoss ;
00862 }
00863
00864
00865
00866
00867 G4double G4hImpactIonisation::AntiProtonParametrisedDEDX(const G4MaterialCutsCouple* couple,
00868 G4double kineticEnergy) const
00869 {
00870 const G4Material* material = couple->GetMaterial();
00871 G4AntiProton* antiproton = G4AntiProton::AntiProton();
00872 G4double eLoss = 0.0 ;
00873
00874
00875 if(antiprotonModel->IsInCharge(antiproton,material)) {
00876 if(kineticEnergy < antiprotonLowEnergy) {
00877 eLoss = antiprotonModel->TheValue(antiproton,material,antiprotonLowEnergy)
00878 * std::sqrt(kineticEnergy/antiprotonLowEnergy) ;
00879
00880
00881 } else {
00882 eLoss = antiprotonModel->TheValue(antiproton,material,
00883 kineticEnergy);
00884 }
00885
00886
00887 } else {
00888 if(kineticEnergy < protonLowEnergy) {
00889 eLoss = protonModel->TheValue(G4Proton::Proton(),material,protonLowEnergy)
00890 * std::sqrt(kineticEnergy/protonLowEnergy) ;
00891
00892
00893 } else {
00894 eLoss = protonModel->TheValue(G4Proton::Proton(),material,
00895 kineticEnergy);
00896 }
00897
00898 }
00899
00900
00901 eLoss -= DeltaRaysEnergy(couple,kineticEnergy,proton_mass_c2) ;
00902
00903 if(verboseLevel > 2) {
00904 G4cout << "pbar E(MeV)= " << kineticEnergy/MeV
00905 << " dE/dx(MeV/mm)= " << eLoss*mm/MeV
00906 << " for " << material->GetName()
00907 << " model: " << protonModel << G4endl;
00908 }
00909
00910 if(eLoss < 0.0) eLoss = 0.0 ;
00911
00912 return eLoss ;
00913 }
00914
00915
00916
00917 G4double G4hImpactIonisation::DeltaRaysEnergy(const G4MaterialCutsCouple* couple,
00918 G4double kineticEnergy,
00919 G4double particleMass) const
00920 {
00921 G4double dLoss = 0.;
00922
00923 G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
00924 const G4Material* material = couple->GetMaterial();
00925 G4double electronDensity = material->GetElectronDensity();
00926 G4double excitationEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
00927
00928 G4double tau = kineticEnergy / particleMass ;
00929 G4double rateMass = electron_mass_c2/particleMass ;
00930
00931
00932
00933 G4double gamma = tau + 1.0 ;
00934 G4double bg2 = tau*(tau+2.0) ;
00935 G4double beta2 = bg2/(gamma*gamma) ;
00936 G4double tMax = 2.*electron_mass_c2*bg2/(1.0+2.0*gamma*rateMass+rateMass*rateMass) ;
00937
00938
00939 G4double deltaCut = std::max(deltaCutNow, excitationEnergy);
00940
00941 if ( deltaCut < tMax)
00942 {
00943 G4double x = deltaCut / tMax ;
00944 dLoss = ( beta2 * (x-1.) - std::log(x) ) * twopi_mc2_rcl2 * electronDensity / beta2 ;
00945 }
00946 return dLoss ;
00947 }
00948
00949
00950
00951
00952 G4VParticleChange* G4hImpactIonisation::PostStepDoIt(const G4Track& track,
00953 const G4Step& step)
00954 {
00955
00956
00957
00958
00959 aParticleChange.Initialize(track) ;
00960 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
00961 const G4DynamicParticle* aParticle = track.GetDynamicParticle() ;
00962
00963
00964
00965 G4ParticleDefinition* definition = track.GetDefinition();
00966 G4double mass = definition->GetPDGMass();
00967 G4double kineticEnergy = aParticle->GetKineticEnergy();
00968 G4double totalEnergy = kineticEnergy + mass ;
00969 G4double pSquare = kineticEnergy *( totalEnergy + mass) ;
00970 G4double eSquare = totalEnergy * totalEnergy;
00971 G4double betaSquare = pSquare / eSquare;
00972 G4ThreeVector particleDirection = aParticle->GetMomentumDirection() ;
00973
00974 G4double gamma = kineticEnergy / mass + 1.;
00975 G4double r = electron_mass_c2 / mass;
00976 G4double tMax = 2. * electron_mass_c2 *(gamma*gamma - 1.) / (1. + 2.*gamma*r + r*r);
00977
00978
00979 G4double deltaCut = cutForDelta[couple->GetIndex()];
00980
00981
00982 if (deltaCut >= tMax)
00983 return G4VContinuousDiscreteProcess::PostStepDoIt(track,step);
00984
00985 G4double xc = deltaCut / tMax;
00986 G4double rate = tMax / totalEnergy;
00987 rate = rate*rate ;
00988 G4double spin = aParticle->GetDefinition()->GetPDGSpin() ;
00989
00990
00991 G4double x = 0.;
00992 G4double gRej = 0.;
00993
00994 do {
00995 x = xc / (1. - (1. - xc) * G4UniformRand());
00996
00997 if (0.0 == spin)
00998 {
00999 gRej = 1.0 - betaSquare * x ;
01000 }
01001 else if (0.5 == spin)
01002 {
01003 gRej = (1. - betaSquare * x + 0.5 * x*x * rate) / (1. + 0.5 * rate) ;
01004 }
01005 else
01006 {
01007 gRej = (1. - betaSquare * x ) * (1. + x/(3.*xc)) +
01008 x*x * rate * (1. + 0.5*x/xc) / 3.0 /
01009 (1. + 1./(3.*xc) + rate *(1.+ 0.5/xc) / 3.);
01010 }
01011
01012 } while( G4UniformRand() > gRej );
01013
01014 G4double deltaKineticEnergy = x * tMax;
01015 G4double deltaTotalMomentum = std::sqrt(deltaKineticEnergy *
01016 (deltaKineticEnergy + 2. * electron_mass_c2 ));
01017 G4double totalMomentum = std::sqrt(pSquare) ;
01018 G4double cosTheta = deltaKineticEnergy * (totalEnergy + electron_mass_c2) / (deltaTotalMomentum*totalMomentum);
01019
01020
01021 if ( cosTheta < -1. ) cosTheta = -1.;
01022 if ( cosTheta > 1. ) cosTheta = 1.;
01023
01024
01025 G4double phi = twopi * G4UniformRand();
01026 G4double sinTheta = std::sqrt(1. - cosTheta*cosTheta);
01027 G4double dirX = sinTheta * std::cos(phi);
01028 G4double dirY = sinTheta * std::sin(phi);
01029 G4double dirZ = cosTheta;
01030
01031 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
01032 deltaDirection.rotateUz(particleDirection);
01033
01034
01035 G4DynamicParticle* deltaRay = new G4DynamicParticle;
01036 deltaRay->SetKineticEnergy( deltaKineticEnergy );
01037 deltaRay->SetMomentumDirection(deltaDirection.x(),
01038 deltaDirection.y(),
01039 deltaDirection.z());
01040 deltaRay->SetDefinition(G4Electron::Electron());
01041
01042
01043 G4double finalKineticEnergy = kineticEnergy - deltaKineticEnergy;
01044 size_t totalNumber = 1;
01045
01046
01047
01048
01049
01050 size_t nSecondaries = 0;
01051 std::vector<G4DynamicParticle*>* secondaryVector = 0;
01052
01053 if (definition == G4Proton::ProtonDefinition())
01054 {
01055 const G4Material* material = couple->GetMaterial();
01056
01057
01058 if (pixeCrossSectionHandler == 0)
01059 {
01060
01061
01062 G4IInterpolator* interpolation = new G4LogLogInterpolator();
01063 pixeCrossSectionHandler =
01064 new G4PixeCrossSectionHandler(interpolation,modelK,modelL,modelM,eMinPixe,eMaxPixe);
01065 G4String fileName("proton");
01066 pixeCrossSectionHandler->LoadShellData(fileName);
01067
01068 }
01069
01070
01071 G4int Z = pixeCrossSectionHandler->SelectRandomAtom(material,kineticEnergy);
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091 G4int shellIndex = pixeCrossSectionHandler->SelectRandomShell(Z,kineticEnergy);
01092
01093 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
01094 const G4AtomicShell* atomicShell = transitionManager->Shell(Z,shellIndex);
01095 G4double bindingEnergy = atomicShell->BindingEnergy();
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110 G4ParticleDefinition* type = 0;
01111
01112 if (finalKineticEnergy >= bindingEnergy)
01113
01114 {
01115
01116 G4int shellId = atomicShell->ShellId();
01117
01118 secondaryVector = atomicDeexcitation.GenerateParticles(Z, shellId);
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128 if (secondaryVector != 0)
01129 {
01130 nSecondaries = secondaryVector->size();
01131 for (size_t i = 0; i<nSecondaries; i++)
01132 {
01133 G4DynamicParticle* aSecondary = (*secondaryVector)[i];
01134 if (aSecondary)
01135 {
01136 G4double e = aSecondary->GetKineticEnergy();
01137 type = aSecondary->GetDefinition();
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149 if (e < finalKineticEnergy &&
01150 ((type == G4Gamma::Gamma() && e > minGammaEnergy ) ||
01151 (type == G4Electron::Electron() && e > minElectronEnergy )))
01152 {
01153
01154 finalKineticEnergy -= e;
01155 totalNumber++;
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165 }
01166 else
01167 {
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180 delete aSecondary;
01181 (*secondaryVector)[i] = 0;
01182 }
01183 }
01184 }
01185 }
01186 }
01187 }
01188
01189
01190
01191
01192 G4double eDeposit = 0.;
01193
01194 if (finalKineticEnergy > MinKineticEnergy)
01195 {
01196 G4double finalPx = totalMomentum*particleDirection.x() - deltaTotalMomentum*deltaDirection.x();
01197 G4double finalPy = totalMomentum*particleDirection.y() - deltaTotalMomentum*deltaDirection.y();
01198 G4double finalPz = totalMomentum*particleDirection.z() - deltaTotalMomentum*deltaDirection.z();
01199 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz) ;
01200 finalPx /= finalMomentum;
01201 finalPy /= finalMomentum;
01202 finalPz /= finalMomentum;
01203
01204 aParticleChange.ProposeMomentumDirection( finalPx,finalPy,finalPz );
01205 }
01206 else
01207 {
01208 eDeposit = finalKineticEnergy;
01209 finalKineticEnergy = 0.;
01210 aParticleChange.ProposeMomentumDirection(particleDirection.x(),
01211 particleDirection.y(),
01212 particleDirection.z());
01213 if(!aParticle->GetDefinition()->GetProcessManager()->
01214 GetAtRestProcessVector()->size())
01215 aParticleChange.ProposeTrackStatus(fStopAndKill);
01216 else
01217 aParticleChange.ProposeTrackStatus(fStopButAlive);
01218 }
01219
01220 aParticleChange.ProposeEnergy(finalKineticEnergy);
01221 aParticleChange.ProposeLocalEnergyDeposit (eDeposit);
01222 aParticleChange.SetNumberOfSecondaries(totalNumber);
01223 aParticleChange.AddSecondary(deltaRay);
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237 if (secondaryVector != 0)
01238 {
01239 for (size_t l = 0; l < nSecondaries; l++)
01240 {
01241 G4DynamicParticle* secondary = (*secondaryVector)[l];
01242 if (secondary) aParticleChange.AddSecondary(secondary);
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257 }
01258 delete secondaryVector;
01259 }
01260
01261 return G4VContinuousDiscreteProcess::PostStepDoIt(track,step);
01262 }
01263
01264
01265
01266 G4double G4hImpactIonisation::ComputeDEDX(const G4ParticleDefinition* aParticle,
01267 const G4MaterialCutsCouple* couple,
01268 G4double kineticEnergy)
01269 {
01270 const G4Material* material = couple->GetMaterial();
01271 G4Proton* proton = G4Proton::Proton();
01272 G4AntiProton* antiproton = G4AntiProton::AntiProton();
01273 G4double dedx = 0.;
01274
01275 G4double tScaled = kineticEnergy * proton_mass_c2 / (aParticle->GetPDGMass()) ;
01276 charge = aParticle->GetPDGCharge() ;
01277
01278 if( charge > 0.)
01279 {
01280 if (tScaled > protonHighEnergy)
01281 {
01282 dedx = G4EnergyLossTables::GetDEDX(proton,tScaled,couple) ;
01283 }
01284 else
01285 {
01286 dedx = ProtonParametrisedDEDX(couple,tScaled) ;
01287 }
01288 }
01289 else
01290 {
01291 if (tScaled > antiprotonHighEnergy)
01292 {
01293 dedx = G4EnergyLossTables::GetDEDX(antiproton,tScaled,couple);
01294 }
01295 else
01296 {
01297 dedx = AntiProtonParametrisedDEDX(couple,tScaled) ;
01298 }
01299 }
01300 dedx *= theIonEffChargeModel->TheValue(aParticle, material, kineticEnergy) ;
01301
01302 return dedx ;
01303 }
01304
01305
01306 G4double G4hImpactIonisation::BarkasTerm(const G4Material* material,
01307 G4double kineticEnergy) const
01308
01309
01310
01311
01312
01313
01314 {
01315 static double FTable[47][2] = {
01316 { 0.02, 21.5},
01317 { 0.03, 20.0},
01318 { 0.04, 18.0},
01319 { 0.05, 15.6},
01320 { 0.06, 15.0},
01321 { 0.07, 14.0},
01322 { 0.08, 13.5},
01323 { 0.09, 13.},
01324 { 0.1, 12.2},
01325 { 0.2, 9.25},
01326 { 0.3, 7.0},
01327 { 0.4, 6.0},
01328 { 0.5, 4.5},
01329 { 0.6, 3.5},
01330 { 0.7, 3.0},
01331 { 0.8, 2.5},
01332 { 0.9, 2.0},
01333 { 1.0, 1.7},
01334 { 1.2, 1.2},
01335 { 1.3, 1.0},
01336 { 1.4, 0.86},
01337 { 1.5, 0.7},
01338 { 1.6, 0.61},
01339 { 1.7, 0.52},
01340 { 1.8, 0.5},
01341 { 1.9, 0.43},
01342 { 2.0, 0.42},
01343 { 2.1, 0.3},
01344 { 2.4, 0.2},
01345 { 3.0, 0.13},
01346 { 3.08, 0.1},
01347 { 3.1, 0.09},
01348 { 3.3, 0.08},
01349 { 3.5, 0.07},
01350 { 3.8, 0.06},
01351 { 4.0, 0.051},
01352 { 4.1, 0.04},
01353 { 4.8, 0.03},
01354 { 5.0, 0.024},
01355 { 5.1, 0.02},
01356 { 6.0, 0.013},
01357 { 6.5, 0.01},
01358 { 7.0, 0.009},
01359 { 7.1, 0.008},
01360 { 8.0, 0.006},
01361 { 9.0, 0.0032},
01362 { 10.0, 0.0025} };
01363
01364
01365 G4double kinE = kineticEnergy ;
01366 if(0.5*MeV > kinE) kinE = 0.5*MeV ;
01367 G4double gamma = 1.0 + kinE / proton_mass_c2 ;
01368 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
01369 if(0.0 >= beta2) return 0.0;
01370
01371 G4double BTerm = 0.0;
01372
01373 G4double ZMaterial = 0.0;
01374 const G4ElementVector* theElementVector = material->GetElementVector();
01375 G4int numberOfElements = material->GetNumberOfElements();
01376
01377 for (G4int i = 0; i<numberOfElements; i++) {
01378
01379
01380 ZMaterial = (*theElementVector)[i]->GetZ();
01381
01382 G4double X = 137.0 * 137.0 * beta2 / ZMaterial;
01383
01384
01385 G4double Eta0Chi = 0.8;
01386 G4double EtaChi = Eta0Chi * ( 1.0 + 6.02*std::pow( ZMaterial,-1.19 ) );
01387 G4double W = ( EtaChi * std::pow( ZMaterial,1.0/6.0 ) ) / std::sqrt(X);
01388 G4double FunctionOfW = FTable[46][1]*FTable[46][0]/W ;
01389
01390 for(G4int j=0; j<47; j++) {
01391
01392 if( W < FTable[j][0] ) {
01393
01394 if(0 == j) {
01395 FunctionOfW = FTable[0][1] ;
01396
01397 } else {
01398 FunctionOfW = (FTable[j][1] - FTable[j-1][1]) * (W - FTable[j-1][0])
01399 / (FTable[j][0] - FTable[j-1][0])
01400 + FTable[j-1][1] ;
01401 }
01402
01403 break;
01404 }
01405
01406 }
01407
01408 BTerm += FunctionOfW /( std::sqrt(ZMaterial * X) * X);
01409 }
01410
01411 BTerm *= twopi_mc2_rcl2 * (material->GetElectronDensity()) / beta2 ;
01412
01413 return BTerm;
01414 }
01415
01416
01417
01418 G4double G4hImpactIonisation::BlochTerm(const G4Material* material,
01419 G4double kineticEnergy,
01420 G4double cSquare) const
01421
01422
01423
01424
01425
01426
01427 {
01428 G4double eLoss = 0.0 ;
01429 G4double gamma = 1.0 + kineticEnergy / proton_mass_c2 ;
01430 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
01431 G4double y = cSquare / (137.0*137.0*beta2) ;
01432
01433 if(y < 0.05) {
01434 eLoss = 1.202 ;
01435
01436 } else {
01437 eLoss = 1.0 / (1.0 + y) ;
01438 G4double de = eLoss ;
01439
01440 for(G4int i=2; de>eLoss*0.01; i++) {
01441 de = 1.0/( i * (i*i + y)) ;
01442 eLoss += de ;
01443 }
01444 }
01445 eLoss *= -1.0 * y * cSquare * twopi_mc2_rcl2 *
01446 (material->GetElectronDensity()) / beta2 ;
01447
01448 return eLoss;
01449 }
01450
01451
01452
01453 G4double G4hImpactIonisation::ElectronicLossFluctuation(
01454 const G4DynamicParticle* particle,
01455 const G4MaterialCutsCouple* couple,
01456 G4double meanLoss,
01457 G4double step) const
01458
01459
01460
01461 {
01462
01463
01464
01465
01466
01467 static const G4double minLoss = 1.*eV ;
01468 static const G4double kappa = 10. ;
01469 static const G4double theBohrBeta2 = 50.0 * keV/proton_mass_c2 ;
01470
01471 const G4Material* material = couple->GetMaterial();
01472 G4int imaterial = couple->GetIndex() ;
01473 G4double ipotFluct = material->GetIonisation()->GetMeanExcitationEnergy() ;
01474 G4double electronDensity = material->GetElectronDensity() ;
01475 G4double zeff = electronDensity/(material->GetTotNbOfAtomsPerVolume()) ;
01476
01477
01478 G4double tkin = particle->GetKineticEnergy();
01479 G4double particleMass = particle->GetMass() ;
01480 G4double deltaCutInKineticEnergyNow = cutForDelta[imaterial];
01481
01482
01483 if(meanLoss < minLoss) return meanLoss ;
01484
01485
01486 G4double threshold = std::max(deltaCutInKineticEnergyNow,ipotFluct);
01487 G4double loss, siga;
01488
01489 G4double rmass = electron_mass_c2/particleMass;
01490 G4double tau = tkin/particleMass;
01491 G4double tau1 = tau+1.0;
01492 G4double tau2 = tau*(tau+2.);
01493 G4double tMax = 2.*electron_mass_c2*tau2/(1.+2.*tau1*rmass+rmass*rmass);
01494
01495
01496 if(tMax > threshold) tMax = threshold;
01497 G4double beta2 = tau2/(tau1*tau1);
01498
01499
01500 if(meanLoss > kappa*tMax || tMax < kappa*ipotFluct )
01501 {
01502 siga = tMax * (1.0-0.5*beta2) * step * twopi_mc2_rcl2
01503 * electronDensity / beta2 ;
01504
01505
01506 if( beta2 > 3.0*theBohrBeta2*zeff || charge < 0.0) {
01507 siga = std::sqrt( siga * chargeSquare ) ;
01508
01509
01510
01511 } else {
01512 G4double chu = theIonChuFluctuationModel->TheValue(particle, material);
01513 G4double yang = theIonYangFluctuationModel->TheValue(particle, material);
01514 siga = std::sqrt( siga * (chargeSquare * chu + yang)) ;
01515 }
01516
01517 do {
01518 loss = G4RandGauss::shoot(meanLoss,siga);
01519 } while (loss < 0. || loss > 2.0*meanLoss);
01520 return loss;
01521 }
01522
01523
01524 static const G4double probLim = 0.01 ;
01525 static const G4double sumaLim = -std::log(probLim) ;
01526 static const G4double alim = 10.;
01527
01528 G4double suma,w1,w2,C,e0,lossc,w;
01529 G4double a1,a2,a3;
01530 G4int p1,p2,p3;
01531 G4int nb;
01532 G4double corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
01533 G4double dp3;
01534
01535 G4double f1Fluct = material->GetIonisation()->GetF1fluct();
01536 G4double f2Fluct = material->GetIonisation()->GetF2fluct();
01537 G4double e1Fluct = material->GetIonisation()->GetEnergy1fluct();
01538 G4double e2Fluct = material->GetIonisation()->GetEnergy2fluct();
01539 G4double e1LogFluct = material->GetIonisation()->GetLogEnergy1fluct();
01540 G4double e2LogFluct = material->GetIonisation()->GetLogEnergy2fluct();
01541 G4double rateFluct = material->GetIonisation()->GetRateionexcfluct();
01542 G4double ipotLogFluct= material->GetIonisation()->GetLogMeanExcEnergy();
01543
01544 w1 = tMax/ipotFluct;
01545 w2 = std::log(2.*electron_mass_c2*tau2);
01546
01547 C = meanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
01548
01549 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
01550 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
01551 a3 = rateFluct*meanLoss*(tMax-ipotFluct)/(ipotFluct*tMax*std::log(w1));
01552 if(a1 < 0.0) a1 = 0.0;
01553 if(a2 < 0.0) a2 = 0.0;
01554 if(a3 < 0.0) a3 = 0.0;
01555
01556 suma = a1+a2+a3;
01557
01558 loss = 0.;
01559
01560
01561 if(suma < sumaLim)
01562 {
01563 e0 = material->GetIonisation()->GetEnergy0fluct();
01564
01565 if(tMax == ipotFluct)
01566 {
01567 a3 = meanLoss/e0;
01568
01569 if(a3>alim)
01570 {
01571 siga=std::sqrt(a3) ;
01572 p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
01573 }
01574 else
01575 p3 = G4Poisson(a3);
01576
01577 loss = p3*e0 ;
01578
01579 if(p3 > 0)
01580 loss += (1.-2.*G4UniformRand())*e0 ;
01581
01582 }
01583 else
01584 {
01585 tMax = tMax-ipotFluct+e0 ;
01586 a3 = meanLoss*(tMax-e0)/(tMax*e0*std::log(tMax/e0));
01587
01588 if(a3>alim)
01589 {
01590 siga=std::sqrt(a3) ;
01591 p3 = std::max(0,int(G4RandGauss::shoot(a3,siga)+0.5));
01592 }
01593 else
01594 p3 = G4Poisson(a3);
01595
01596 if(p3 > 0)
01597 {
01598 w = (tMax-e0)/tMax ;
01599 if(p3 > nmaxCont2)
01600 {
01601 dp3 = G4float(p3) ;
01602 corrfac = dp3/G4float(nmaxCont2) ;
01603 p3 = nmaxCont2 ;
01604 }
01605 else
01606 corrfac = 1. ;
01607
01608 for(G4int i=0; i<p3; i++) loss += 1./(1.-w*G4UniformRand()) ;
01609 loss *= e0*corrfac ;
01610 }
01611 }
01612 }
01613
01614 else
01615 {
01616
01617 if(a1>alim)
01618 {
01619 siga=std::sqrt(a1) ;
01620 p1 = std::max(0,G4int(G4RandGauss::shoot(a1,siga)+0.5));
01621 }
01622 else
01623 p1 = G4Poisson(a1);
01624
01625
01626 if(a2>alim)
01627 {
01628 siga=std::sqrt(a2) ;
01629 p2 = std::max(0,G4int(G4RandGauss::shoot(a2,siga)+0.5));
01630 }
01631 else
01632 p2 = G4Poisson(a2);
01633
01634 loss = p1*e1Fluct+p2*e2Fluct;
01635
01636
01637 if(p2 > 0)
01638 loss += (1.-2.*G4UniformRand())*e2Fluct;
01639 else if (loss>0.)
01640 loss += (1.-2.*G4UniformRand())*e1Fluct;
01641
01642
01643 if(a3 > 0.)
01644 {
01645 if(a3>alim)
01646 {
01647 siga=std::sqrt(a3) ;
01648 p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
01649 }
01650 else
01651 p3 = G4Poisson(a3);
01652
01653 lossc = 0.;
01654 if(p3 > 0)
01655 {
01656 na = 0.;
01657 alfa = 1.;
01658 if (p3 > nmaxCont2)
01659 {
01660 dp3 = G4float(p3);
01661 rfac = dp3/(G4float(nmaxCont2)+dp3);
01662 namean = G4float(p3)*rfac;
01663 sa = G4float(nmaxCont1)*rfac;
01664 na = G4RandGauss::shoot(namean,sa);
01665 if (na > 0.)
01666 {
01667 alfa = w1*G4float(nmaxCont2+p3)/
01668 (w1*G4float(nmaxCont2)+G4float(p3));
01669 alfa1 = alfa*std::log(alfa)/(alfa-1.);
01670 ea = na*ipotFluct*alfa1;
01671 sea = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
01672 lossc += G4RandGauss::shoot(ea,sea);
01673 }
01674 }
01675
01676 nb = G4int(G4float(p3)-na);
01677 if (nb > 0)
01678 {
01679 w2 = alfa*ipotFluct;
01680 w = (tMax-w2)/tMax;
01681 for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
01682 }
01683 }
01684 loss += lossc;
01685 }
01686 }
01687
01688 return loss ;
01689 }
01690
01691
01692
01693 void G4hImpactIonisation::SetCutForSecondaryPhotons(G4double cut)
01694 {
01695 minGammaEnergy = cut;
01696 }
01697
01698
01699
01700 void G4hImpactIonisation::SetCutForAugerElectrons(G4double cut)
01701 {
01702 minElectronEnergy = cut;
01703 }
01704
01705
01706
01707 void G4hImpactIonisation::ActivateAugerElectronProduction(G4bool val)
01708 {
01709 atomicDeexcitation.ActivateAugerElectronProduction(val);
01710 }
01711
01712
01713
01714 void G4hImpactIonisation::PrintInfoDefinition() const
01715 {
01716 G4String comments = " Knock-on electron cross sections . ";
01717 comments += "\n Good description above the mean excitation energy.\n";
01718 comments += " Delta ray energy sampled from differential Xsection.";
01719
01720 G4cout << G4endl << GetProcessName() << ": " << comments
01721 << "\n PhysicsTables from " << LowestKineticEnergy / eV << " eV "
01722 << " to " << HighestKineticEnergy / TeV << " TeV "
01723 << " in " << TotBin << " bins."
01724 << "\n Electronic stopping power model is "
01725 << protonTable
01726 << "\n from " << protonLowEnergy / keV << " keV "
01727 << " to " << protonHighEnergy / MeV << " MeV " << "." << G4endl ;
01728 G4cout << "\n Parametrisation model for antiprotons is "
01729 << antiprotonTable
01730 << "\n from " << antiprotonLowEnergy / keV << " keV "
01731 << " to " << antiprotonHighEnergy / MeV << " MeV " << "." << G4endl ;
01732 if(theBarkas){
01733 G4cout << " Parametrization of the Barkas effect is switched on."
01734 << G4endl ;
01735 }
01736 if(nStopping) {
01737 G4cout << " Nuclear stopping power model is " << theNuclearTable
01738 << G4endl ;
01739 }
01740
01741 G4bool printHead = true;
01742
01743 const G4ProductionCutsTable* theCoupleTable=
01744 G4ProductionCutsTable::GetProductionCutsTable();
01745 size_t numOfCouples = theCoupleTable->GetTableSize();
01746
01747
01748
01749 for (size_t j=0 ; j < numOfCouples; j++) {
01750
01751 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
01752 const G4Material* material= couple->GetMaterial();
01753 G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
01754 G4double excitationEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
01755
01756 if(excitationEnergy > deltaCutNow) {
01757 if(printHead) {
01758 printHead = false ;
01759
01760 G4cout << " material min.delta energy(keV) " << G4endl;
01761 G4cout << G4endl;
01762 }
01763
01764 G4cout << std::setw(20) << material->GetName()
01765 << std::setw(15) << excitationEnergy/keV << G4endl;
01766 }
01767 }
01768 }