Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 28 // ------------------------------------------- 29 // G4RDHadronIonisation 30 // 31 // 32 // Author: Maria Grazia Pia (MariaGrazia.Pia@g 33 // 34 // 08 Sep 2008 - MGP - Created (initially base 35 // Added PIXE capabilities 36 // Partial clean-up of the 37 // Calculation of Microsco 38 // M.G. Pia et al., PIXE Simulation With Geant 39 // IEEE Trans. Nucl. Sci., vol. 56, no. 6, pp. 40 41 // 42 // ------------------------------------------- 43 44 #include "G4hImpactIonisation.hh" 45 #include "globals.hh" 46 #include "G4ios.hh" 47 #include "Randomize.hh" 48 #include "G4PhysicalConstants.hh" 49 #include "G4SystemOfUnits.hh" 50 #include "G4Poisson.hh" 51 #include "G4UnitsTable.hh" 52 #include "G4EnergyLossTables.hh" 53 #include "G4Material.hh" 54 #include "G4DynamicParticle.hh" 55 #include "G4ParticleDefinition.hh" 56 #include "G4AtomicDeexcitation.hh" 57 #include "G4AtomicTransitionManager.hh" 58 #include "G4PixeCrossSectionHandler.hh" 59 #include "G4IInterpolator.hh" 60 #include "G4LogLogInterpolator.hh" 61 #include "G4Gamma.hh" 62 #include "G4Electron.hh" 63 #include "G4Proton.hh" 64 #include "G4ProcessManager.hh" 65 #include "G4ProductionCutsTable.hh" 66 #include "G4PhysicsLogVector.hh" 67 #include "G4PhysicsLinearVector.hh" 68 69 #include "G4VLowEnergyModel.hh" 70 #include "G4hNuclearStoppingModel.hh" 71 #include "G4hBetheBlochModel.hh" 72 #include "G4hParametrisedLossModel.hh" 73 #include "G4QAOLowEnergyLoss.hh" 74 #include "G4hIonEffChargeSquare.hh" 75 #include "G4IonChuFluctuationModel.hh" 76 #include "G4IonYangFluctuationModel.hh" 77 78 #include "G4MaterialCutsCouple.hh" 79 #include "G4Track.hh" 80 #include "G4Step.hh" 81 82 G4hImpactIonisation::G4hImpactIonisation(const 83 : G4hRDEnergyLoss(processName), 84 betheBlochModel(0), 85 protonModel(0), 86 antiprotonModel(0), 87 theIonEffChargeModel(0), 88 theNuclearStoppingModel(0), 89 theIonChuFluctuationModel(0), 90 theIonYangFluctuationModel(0), 91 protonTable("ICRU_R49p"), 92 antiprotonTable("ICRU_R49p"), 93 theNuclearTable("ICRU_R49"), 94 nStopping(true), 95 theBarkas(true), 96 theMeanFreePathTable(0), 97 paramStepLimit (0.005), 98 pixeCrossSectionHandler(0) 99 { 100 InitializeMe(); 101 } 102 103 104 105 void G4hImpactIonisation::InitializeMe() 106 { 107 LowestKineticEnergy = 10.0*eV ; 108 HighestKineticEnergy = 100.0*GeV ; 109 MinKineticEnergy = 10.0*eV ; 110 TotBin = 360 ; 111 protonLowEnergy = 1.*keV ; 112 protonHighEnergy = 100.*MeV ; 113 antiprotonLowEnergy = 25.*keV ; 114 antiprotonHighEnergy = 2.*MeV ; 115 minGammaEnergy = 100 * eV; 116 minElectronEnergy = 250.* eV; 117 verboseLevel = 0; 118 119 // Min and max energy of incident particle f 120 // for PIXE generation 121 eMinPixe = 1.* keV; 122 eMaxPixe = 200. * MeV; 123 124 G4String defaultPixeModel("ecpssr"); 125 modelK = defaultPixeModel; 126 modelL = defaultPixeModel; 127 modelM = defaultPixeModel; 128 129 // Calculated on the fly, but should be init 130 fdEdx = 0.0; 131 fRangeNow = 0.0; 132 charge = 0.0; 133 chargeSquare = 0.0; 134 initialMass = 0.0; 135 fBarkas = 0.0; 136 } 137 138 139 G4hImpactIonisation::~G4hImpactIonisation() 140 { 141 if (theMeanFreePathTable) 142 { 143 theMeanFreePathTable->clearAndDestroy(); 144 delete theMeanFreePathTable; 145 } 146 147 if (betheBlochModel) delete betheBlochModel; 148 if (protonModel) delete protonModel; 149 if (antiprotonModel) delete antiprotonModel; 150 if (theNuclearStoppingModel) delete theNucle 151 if (theIonEffChargeModel) delete theIonEffCh 152 if (theIonChuFluctuationModel) delete theIon 153 if (theIonYangFluctuationModel) delete theIo 154 155 delete pixeCrossSectionHandler; 156 157 // ---- MGP ---- The following is to be chec 158 // if (shellVacancy) delete shellVacancy; 159 160 cutForDelta.clear(); 161 } 162 163 // ------------------------------------------- 164 void G4hImpactIonisation::SetElectronicStoppin 165 const G4String& dedxTable) 166 // This method defines the ionisation parametr 167 { 168 if (particle->GetPDGCharge() > 0 ) 169 { 170 // Positive charge 171 SetProtonElectronicStoppingPowerModel(de 172 } 173 else 174 { 175 // Antiprotons 176 SetAntiProtonElectronicStoppingPowerMode 177 } 178 } 179 180 181 // ------------------------------------------- 182 void G4hImpactIonisation::InitializeParametris 183 184 { 185 // Define models for parametrisation of elec 186 betheBlochModel = new G4hBetheBlochModel("Be 187 protonModel = new G4hParametrisedLossModel(p 188 protonHighEnergy = std::min(protonHighEnergy 189 antiprotonModel = new G4QAOLowEnergyLoss(ant 190 theNuclearStoppingModel = new G4hNuclearStop 191 theIonEffChargeModel = new G4hIonEffChargeSq 192 theIonChuFluctuationModel = new G4IonChuFluc 193 theIonYangFluctuationModel = new G4IonYangFl 194 } 195 196 197 // ------------------------------------------- 198 void G4hImpactIonisation::BuildPhysicsTable(co 199 200 // just call BuildLossTable+BuildLambdaTable 201 { 202 203 // Verbose print-out 204 if(verboseLevel > 0) 205 { 206 G4cout << "G4hImpactIonisation::BuildPhy 207 << particleDef.GetParticleName() 208 << " mass(MeV)= " << particleDef.GetPDG 209 << " charge= " << particleDef.GetPDGCha 210 << " type= " << particleDef.GetParticle 211 << G4endl; 212 213 if(verboseLevel > 1) 214 { 215 G4ProcessVector* pv = particleDef.GetProce 216 217 G4cout << " 0: " << (*pv)[0]->GetProcessNa 218 << " 1: " << (*pv)[1]->GetProcessName() < 219 // << " 2: " << (*pv)[2]->GetProc 220 << G4endl; 221 G4cout << "ionModel= " << theIonEffChargeM 222 << " MFPtable= " << theMeanFreePathTable 223 << " iniMass= " << initialMass 224 << G4endl; 225 } 226 } 227 // End of verbose print-out 228 229 if (particleDef.GetParticleType() == "nucleu 230 particleDef.GetParticleName() != "Generi 231 particleDef.GetParticleSubType() == "gen 232 { 233 234 G4EnergyLossTables::Register(&particleDe 235 theDEDXpTable, 236 theRangepTable, 237 theInverseRangepTable, 238 theLabTimepTable, 239 theProperTimepTable, 240 LowestKineticEnergy, HighestKinetic 241 proton_mass_c2/particleDef.GetPDGMa 242 TotBin); 243 244 return; 245 } 246 247 if( !CutsWhereModified() && theLossTable) re 248 249 InitializeParametrisation() ; 250 G4Proton* proton = G4Proton::Proton(); 251 G4AntiProton* antiproton = G4AntiProton::Ant 252 253 charge = particleDef.GetPDGCharge() / eplus; 254 chargeSquare = charge*charge ; 255 256 const G4ProductionCutsTable* theCoupleTable= 257 G4ProductionCutsTable::GetProductionCutsTa 258 G4int numOfCouples = (G4int)theCoupleTable-> 259 260 cutForDelta.clear(); 261 cutForGamma.clear(); 262 263 for (G4int j=0; j<numOfCouples; ++j) { 264 265 // get material parameters needed for the 266 const G4MaterialCutsCouple* couple = theCo 267 const G4Material* material= couple->GetMat 268 269 // the cut cannot be below lowest limit 270 G4double tCut = (*(theCoupleTable->GetEner 271 if(tCut > HighestKineticEnergy) tCut = Hig 272 273 G4double excEnergy = material->GetIonisati 274 275 tCut = std::max(tCut,excEnergy); 276 cutForDelta.push_back(tCut); 277 278 // the cut cannot be below lowest limit 279 tCut = (*(theCoupleTable->GetEnergyCutsVec 280 if(tCut > HighestKineticEnergy) tCut = Hig 281 tCut = std::max(tCut,minGammaEnergy); 282 cutForGamma.push_back(tCut); 283 } 284 285 if(verboseLevel > 0) { 286 G4cout << "Cuts are defined " << G4endl; 287 } 288 289 if(0.0 < charge) 290 { 291 { 292 BuildLossTable(*proton) ; 293 294 // The following vector has a fixed dim 295 // It happended in the past that caused 296 // G4cout << "[NOTE]: __LINE__=" << _ 297 298 RecorderOfpProcess[CounterOfpProcess] 299 CounterOfpProcess++; 300 } 301 } else { 302 { 303 BuildLossTable(*antiproton) ; 304 305 // The following vector has a fixed 306 // It happended in the past that ca 307 // G4cout << "[NOTE]: __LINE__=" 308 309 RecorderOfpbarProcess[CounterOfpbarProce 310 CounterOfpbarProcess++; 311 } 312 } 313 314 if(verboseLevel > 0) { 315 G4cout << "G4hImpactIonisation::BuildPhysi 316 << "Loss table is built " 317 // << theLossTable 318 << G4endl; 319 } 320 321 BuildLambdaTable(particleDef) ; 322 // BuildDataForFluorescence(particleDef); 323 324 if(verboseLevel > 1) { 325 G4cout << (*theMeanFreePathTable) << G4end 326 } 327 328 if(verboseLevel > 0) { 329 G4cout << "G4hImpactIonisation::BuildPhysi 330 << "DEDX table will be built " 331 // << theDEDXpTable << " " << theDED 332 // << " " << theRangepTable << " " < 333 << G4endl; 334 } 335 336 BuildDEDXTable(particleDef) ; 337 338 if(verboseLevel > 1) { 339 G4cout << (*theDEDXpTable) << G4endl; 340 } 341 342 if((&particleDef == proton) || (&particleDe 343 344 if(verboseLevel > 0) { 345 G4cout << "G4hImpactIonisation::BuildPhysi 346 << particleDef.GetParticleName() << 347 } 348 } 349 350 351 352 353 354 // ------------------------------------------- 355 void G4hImpactIonisation::BuildLossTable(const 356 { 357 // Initialisation 358 G4double lowEdgeEnergy , ionloss, ionlossBB, 359 //G4double lowEnergy, highEnergy; 360 G4double highEnergy; 361 G4Proton* proton = G4Proton::Proton(); 362 363 if (particleDef == *proton) 364 { 365 //lowEnergy = protonLowEnergy ; 366 highEnergy = protonHighEnergy ; 367 charge = 1. ; 368 } 369 else 370 { 371 //lowEnergy = antiprotonLowEnergy ; 372 highEnergy = antiprotonHighEnergy ; 373 charge = -1. ; 374 } 375 chargeSquare = 1. ; 376 377 const G4ProductionCutsTable* theCoupleTable= 378 G4ProductionCutsTable::GetProductionCutsTa 379 G4int numOfCouples = (G4int)theCoupleTable-> 380 381 if ( theLossTable) 382 { 383 theLossTable->clearAndDestroy(); 384 delete theLossTable; 385 } 386 387 theLossTable = new G4PhysicsTable(numOfCoupl 388 389 // loop for materials 390 for (G4int j=0; j<numOfCouples; ++j) { 391 392 // create physics vector and fill it 393 G4PhysicsLogVector* aVector = new G4Physic 394 395 396 397 // get material parameters needed for the 398 const G4MaterialCutsCouple* couple = theCo 399 const G4Material* material= couple->GetMat 400 401 if ( charge > 0.0 ) { 402 ionloss = ProtonParametrisedDEDX(couple, 403 } else { 404 ionloss = AntiProtonParametrisedDEDX(cou 405 } 406 407 ionlossBB = betheBlochModel->TheValue(&par 408 ionlossBB -= DeltaRaysEnergy(couple,highEn 409 410 411 paramB = ionloss/ionlossBB - 1.0 ; 412 413 // now comes the loop for the kinetic ener 414 for (G4int i = 0 ; i < TotBin ; i++) { 415 lowEdgeEnergy = aVector->GetLowEdgeEnerg 416 417 // low energy part for this material, pa 418 if ( lowEdgeEnergy < highEnergy ) { 419 420 if ( charge > 0.0 ) { 421 ionloss = ProtonParametrisedDEDX(cou 422 } else { 423 ionloss = AntiProtonParametrisedDEDX 424 } 425 426 } else { 427 428 // high energy part for this material, 429 ionloss = betheBlochModel->TheValue(pr 430 lowEdgeEnergy) ; 431 432 ionloss -= DeltaRaysEnergy(couple,lowE 433 434 ionloss *= (1.0 + paramB*highEnergy/lowEdgeE 435 } 436 437 // now put the loss into the vector 438 if(verboseLevel > 1) { 439 G4cout << "E(MeV)= " << lowEdgeEnergy/ 440 << " dE/dx(MeV/mm)= " << ionlo 441 << " in " << material->GetName( 442 } 443 aVector->PutValue(i,ionloss) ; 444 } 445 // Insert vector for this material into th 446 theLossTable->insert(aVector) ; 447 } 448 } 449 450 451 452 // ------------------------------------------- 453 void G4hImpactIonisation::BuildLambdaTable(con 454 455 { 456 // Build mean free path tables for the delta 457 // tables are built for MATERIALS 458 459 if(verboseLevel > 1) { 460 G4cout << "G4hImpactIonisation::BuildLambd 461 << particleDef.GetParticleName() << 462 } 463 464 465 G4double lowEdgeEnergy, value; 466 charge = particleDef.GetPDGCharge()/eplus ; 467 chargeSquare = charge*charge ; 468 initialMass = particleDef.GetPDGMass(); 469 470 const G4ProductionCutsTable* theCoupleTable= 471 G4ProductionCutsTable::GetProductionCutsTa 472 G4int numOfCouples = (G4int)theCoupleTable-> 473 474 475 if (theMeanFreePathTable) { 476 theMeanFreePathTable->clearAndDestroy(); 477 delete theMeanFreePathTable; 478 } 479 480 theMeanFreePathTable = new G4PhysicsTable(nu 481 482 // loop for materials 483 484 for (G4int j=0 ; j < numOfCouples; ++j) { 485 486 //create physics vector then fill it .... 487 G4PhysicsLogVector* aVector = new G4Physic 488 489 490 491 // compute the (macroscopic) cross section 492 const G4MaterialCutsCouple* couple = theCo 493 const G4Material* material= couple->GetMat 494 495 const G4ElementVector* theElementVector = 496 const G4double* theAtomicNumDensityVector 497 const G4int numberOfElements = (G4int)mate 498 499 // get the electron kinetic energy cut for 500 // it will be used in ComputeMicroscopicC 501 // ( it is the SAME for ALL the ELEMENTS i 502 // ------------------------------------- 503 504 G4double deltaCut = cutForDelta[j]; 505 506 for ( G4int i = 0 ; i < TotBin ; i++ ) { 507 lowEdgeEnergy = aVector->GetLowEdgeEnerg 508 G4double sigma = 0.0 ; 509 G4int Z; 510 511 for (G4int iel=0; iel<numberOfElements; 512 { 513 Z = (G4int) (*theElementVector)[iel]->GetZ 514 // ---- MGP --- Corrected duplicated cross 515 // G4hLowEnergyIonisation original 516 G4double microCross = MicroscopicCrossSect 517 lowEdgeEnergy, 518 Z, 519 deltaCut ) ; 520 //totalCrossSectionMap [Z] = microCross; 521 sigma += theAtomicNumDensityVector[iel] * 522 } 523 524 // mean free path = 1./macroscopic cross 525 526 value = sigma<=0 ? DBL_MAX : 1./sigma ; 527 528 aVector->PutValue(i, value) ; 529 } 530 531 theMeanFreePathTable->insert(aVector); 532 } 533 534 } 535 536 537 // ------------------------------------------- 538 G4double G4hImpactIonisation::MicroscopicCross 539 G4double kineticEnergy, 540 G4double atomicNumber, 541 G4double deltaCutInEnergy) c 542 { 543 //****************************************** 544 // cross section formula is OK for spin=0, 545 // *************************************** 546 547 // Calculates the microscopic cross sectio 548 // Formula documented in Geant4 Phys. Ref. 549 // ( it is called for elements, AtomicNumb 550 551 G4double totalCrossSection = 0.; 552 553 // Particle mass and energy 554 G4double particleMass = initialMass; 555 G4double energy = kineticEnergy + particleMa 556 557 // Some kinematics 558 G4double gamma = energy / particleMass; 559 G4double beta2 = 1. - 1. / (gamma * gamma); 560 G4double var = electron_mass_c2 / particleMa 561 G4double tMax = 2. * electron_mass_c2 * (g 562 563 // Calculate the total cross section 564 565 if ( tMax > deltaCutInEnergy ) 566 { 567 var = deltaCutInEnergy / tMax; 568 totalCrossSection = (1. - var * (1. - be 569 570 G4double spin = particleDef.GetPDGSpin() 571 572 // +term for spin=1/2 particle 573 if (spin == 0.5) 574 { 575 totalCrossSection += 0.5 * (tMax - deltaC 576 } 577 // +term for spin=1 particle 578 else if (spin > 0.9 ) 579 { 580 totalCrossSection += -std::log(var) / 581 (3. * deltaCutInEnergy) + (tMax - deltaC 582 beta2 / (tMax * deltaCutIn 583 } 584 totalCrossSection *= twopi_mc2_rcl2 * at 585 } 586 587 //std::cout << "Microscopic = " << totalCros 588 // << ", e = " << kineticEnergy/MeV <<std 589 590 return totalCrossSection ; 591 } 592 593 594 595 // ------------------------------------------- 596 G4double G4hImpactIonisation::GetMeanFreePath( 597 G4double, // previousStepSize 598 enum G4ForceCondition* conditi 599 { 600 const G4DynamicParticle* dynamicParticle = t 601 const G4MaterialCutsCouple* couple = track.G 602 const G4Material* material = couple->GetMate 603 604 G4double meanFreePath = DBL_MAX; 605 // ---- MGP ---- What is the meaning of the 606 G4bool isOutRange = false; 607 608 *condition = NotForced ; 609 610 G4double kineticEnergy = (dynamicParticle->G 611 charge = dynamicParticle->GetCharge()/eplus; 612 chargeSquare = theIonEffChargeModel->TheValu 613 614 if (kineticEnergy < LowestKineticEnergy) 615 { 616 meanFreePath = DBL_MAX; 617 } 618 else 619 { 620 if (kineticEnergy > HighestKineticEnergy 621 meanFreePath = (((*theMeanFreePathTable) 622 GetValue(kineticEnergy,isOutRange))/ 623 } 624 625 return meanFreePath ; 626 } 627 628 629 // ------------------------------------------- 630 G4double G4hImpactIonisation::GetConstraints(c 631 const G4MaterialCutsCouple* cou 632 { 633 // returns the Step limit 634 // dEdx is calculated as well as the range 635 // based on Effective Charge Approach 636 637 const G4Material* material = couple->GetMate 638 G4Proton* proton = G4Proton::Proton(); 639 G4AntiProton* antiproton = G4AntiProton::Ant 640 641 G4double stepLimit = 0.; 642 G4double dx, highEnergy; 643 644 G4double massRatio = proton_mass_c2/(particl 645 G4double kineticEnergy = particle->GetKineti 646 647 // Scale the kinetic energy 648 649 G4double tScaled = kineticEnergy*massRatio ; 650 fBarkas = 0.; 651 652 if (charge > 0.) 653 { 654 highEnergy = protonHighEnergy ; 655 fRangeNow = G4EnergyLossTables::GetRange 656 dx = G4EnergyLossTables::GetRange(proton 657 fdEdx = G4EnergyLossTables::GetDEDX(prot 658 * chargeSquare ; 659 660 // Correction for positive ions 661 if (theBarkas && tScaled > highEnergy) 662 { 663 fBarkas = BarkasTerm(material,tScaled)*std 664 + BlochTerm(material,tScaled,chargeSquar 665 } 666 // Antiprotons and negative hadrons 667 } 668 else 669 { 670 highEnergy = antiprotonHighEnergy ; 671 fRangeNow = G4EnergyLossTables::GetRange 672 dx = G4EnergyLossTables::GetRange(antipr 673 fdEdx = G4EnergyLossTables::GetDEDX(anti 674 675 if (theBarkas && tScaled > highEnergy) 676 { 677 fBarkas = -BarkasTerm(material,tScaled)*st 678 + BlochTerm(material,tScaled,chargeSquar 679 } 680 } 681 /* 682 const G4Material* mat = couple->GetMateria 683 G4double fac = gram/(MeV*cm2*mat->GetDensi 684 G4cout << particle->GetDefinition()->GetPa 685 << " in " << mat->GetName() 686 << " E(MeV)= " << kineticEnergy/MeV 687 << " dedx(MeV*cm^2/g)= " << fdEdx*fac 688 << " barcas(MeV*cm^2/gram)= " << fBarkas*f 689 << " Q^2= " << chargeSquare 690 << G4endl; 691 */ 692 // scaling back 693 fRangeNow /= (chargeSquare*massRatio) ; 694 dx /= (chargeSquare*massRatio) ; 695 696 stepLimit = fRangeNow ; 697 G4double r = std::min(finalRange, couple->Ge 698 ->GetProductionCut(idxG4ElectronCut)); 699 700 if (fRangeNow > r) 701 { 702 stepLimit = dRoverRange*fRangeNow + r*(1 703 if (stepLimit > fRangeNow) stepLimit = f 704 } 705 // compute the (random) Step limit in stan 706 if(tScaled > highEnergy ) 707 { 708 // add Barkas correction directly to ded 709 fdEdx += fBarkas; 710 711 if(stepLimit > fRangeNow - dx*0.9) stepL 712 713 // Step limit in low energy range 714 } 715 else 716 { 717 G4double x = dx*paramStepLimit; 718 if (stepLimit > x) stepLimit = x; 719 } 720 return stepLimit; 721 } 722 723 724 // ------------------------------------------- 725 G4VParticleChange* G4hImpactIonisation::AlongS 726 const G4Step& step) 727 { 728 // compute the energy loss after a step 729 G4Proton* proton = G4Proton::Proton(); 730 G4AntiProton* antiproton = G4AntiProton::Ant 731 G4double finalT = 0.; 732 733 aParticleChange.Initialize(track) ; 734 735 const G4MaterialCutsCouple* couple = track.G 736 const G4Material* material = couple->GetMate 737 738 // get the actual (true) Step length from st 739 const G4double stepLength = step.GetStepLeng 740 741 const G4DynamicParticle* particle = track.Ge 742 743 G4double kineticEnergy = particle->GetKineti 744 G4double massRatio = proton_mass_c2/(particl 745 G4double tScaled = kineticEnergy * massRatio 746 G4double eLoss = 0.0 ; 747 G4double nLoss = 0.0 ; 748 749 750 // very small particle energy 751 if (kineticEnergy < MinKineticEnergy) 752 { 753 eLoss = kineticEnergy ; 754 // particle energy outside tabulated ene 755 } 756 757 else if( kineticEnergy > HighestKineticEnerg 758 { 759 eLoss = stepLength * fdEdx ; 760 // big step 761 } 762 else if (stepLength >= fRangeNow ) 763 { 764 eLoss = kineticEnergy ; 765 766 // tabulated range 767 } 768 else 769 { 770 // step longer than linear step limit 771 if(stepLength > linLossLimit * fRangeNow 772 { 773 G4double rScaled = fRangeNow * massRatio * 774 G4double sScaled = stepLength * massRatio 775 776 if(charge > 0.0) 777 { 778 eLoss = G4EnergyLossTables::GetPrecise 779 G4EnergyLossTables::GetPreciseEnergyFromRa 780 781 } 782 else 783 { 784 // Antiproton 785 eLoss = G4EnergyLossTables::GetPrecise 786 G4EnergyLossTables::GetPreciseEnergyFromRa 787 } 788 eLoss /= massRatio ; 789 790 // Barkas correction at big step 791 eLoss += fBarkas * stepLength; 792 793 // step shorter than linear step limit 794 } 795 else 796 { 797 eLoss = stepLength *fdEdx ; 798 } 799 if (nStopping && tScaled < protonHighEne 800 { 801 nLoss = (theNuclearStoppingModel->TheValue 802 } 803 } 804 805 if (eLoss < 0.0) eLoss = 0.0; 806 807 finalT = kineticEnergy - eLoss - nLoss; 808 809 if ( EnlossFlucFlag && 0.0 < eLoss && finalT 810 { 811 812 // now the electron loss with fluctuati 813 eLoss = ElectronicLossFluctuation(partic 814 if (eLoss < 0.0) eLoss = 0.0; 815 finalT = kineticEnergy - eLoss - nLoss; 816 } 817 818 // stop particle if the kinetic energy <= M 819 if (finalT*massRatio <= MinKineticEnergy ) 820 { 821 822 finalT = 0.0; 823 if (!particle->GetDefinition()->GetProce 824 aParticleChange.ProposeTrackStatus(fSt 825 else 826 aParticleChange.ProposeTrackStatus(fSt 827 } 828 829 aParticleChange.ProposeEnergy( finalT ); 830 eLoss = kineticEnergy-finalT; 831 832 aParticleChange.ProposeLocalEnergyDeposit(eL 833 return &aParticleChange ; 834 } 835 836 837 838 // ------------------------------------------- 839 G4double G4hImpactIonisation::ProtonParametris 840 G4double kineticEnergy) const 841 { 842 const G4Material* material = couple->GetMate 843 G4Proton* proton = G4Proton::Proton(); 844 G4double eLoss = 0.; 845 846 // Free Electron Gas Model 847 if(kineticEnergy < protonLowEnergy) { 848 eLoss = (protonModel->TheValue(proton, mat 849 * std::sqrt(kineticEnergy/protonLowEnerg 850 851 // Parametrisation 852 } else { 853 eLoss = protonModel->TheValue(proton, mate 854 } 855 856 // Delta rays energy 857 eLoss -= DeltaRaysEnergy(couple,kineticEnerg 858 859 if(verboseLevel > 2) { 860 G4cout << "p E(MeV)= " << kineticEnergy/Me 861 << " dE/dx(MeV/mm)= " << eLoss*mm/M 862 << " for " << material->GetName() 863 << " model: " << protonModel << G4e 864 } 865 866 if(eLoss < 0.0) eLoss = 0.0 ; 867 868 return eLoss ; 869 } 870 871 872 873 // ------------------------------------------- 874 G4double G4hImpactIonisation::AntiProtonParame 875 G4double kineticEnergy) const 876 { 877 const G4Material* material = couple->GetMate 878 G4AntiProton* antiproton = G4AntiProton::Ant 879 G4double eLoss = 0.0 ; 880 881 // Antiproton model is used 882 if(antiprotonModel->IsInCharge(antiproton,ma 883 if(kineticEnergy < antiprotonLowEnergy) { 884 eLoss = antiprotonModel->TheValue(antipr 885 * std::sqrt(kineticEnergy/antiprotonLowEnerg 886 887 // Parametrisation 888 } else { 889 eLoss = antiprotonModel->TheValue(antipr 890 kineticEnergy); 891 } 892 893 // The proton model is used + Barkas corre 894 } else { 895 if(kineticEnergy < protonLowEnergy) { 896 eLoss = protonModel->TheValue(G4Proton:: 897 * std::sqrt(kineticEnergy/protonLowEnergy) ; 898 899 // Parametrisation 900 } else { 901 eLoss = protonModel->TheValue(G4Proton:: 902 kineticEnergy); 903 } 904 //if(theBarkas) eLoss -= 2.0*BarkasTerm(ma 905 } 906 907 // Delta rays energy 908 eLoss -= DeltaRaysEnergy(couple,kineticEnerg 909 910 if(verboseLevel > 2) { 911 G4cout << "pbar E(MeV)= " << kineticEnergy 912 << " dE/dx(MeV/mm)= " << eLoss*mm/M 913 << " for " << material->GetName() 914 << " model: " << protonModel << G4e 915 } 916 917 if(eLoss < 0.0) eLoss = 0.0 ; 918 919 return eLoss ; 920 } 921 922 923 // ------------------------------------------- 924 G4double G4hImpactIonisation::DeltaRaysEnergy( 925 G4double kineticEnergy, 926 G4double particleMass) const 927 { 928 G4double dLoss = 0.; 929 930 G4double deltaCutNow = cutForDelta[(couple-> 931 const G4Material* material = couple->GetMate 932 G4double electronDensity = material->GetElec 933 G4double excitationEnergy = material->GetIon 934 935 G4double tau = kineticEnergy / particleMass 936 G4double rateMass = electron_mass_c2/particl 937 938 // some local variables 939 940 G4double gamma = tau + 1.0 ; 941 G4double bg2 = tau*(tau+2.0) ; 942 G4double beta2 = bg2/(gamma*gamma) ; 943 G4double tMax = 2.*electron_mass_c2*bg2/(1.0 944 945 // Validity range for delta electron cross s 946 G4double deltaCut = std::max(deltaCutNow, ex 947 948 if ( deltaCut < tMax) 949 { 950 G4double x = deltaCut / tMax ; 951 dLoss = ( beta2 * (x-1.) - std::log(x) ) 952 } 953 return dLoss ; 954 } 955 956 957 // ------------------------------------------- 958 959 G4VParticleChange* G4hImpactIonisation::PostSt 960 const G4Step& step) 961 { 962 // Units are expressed in GEANT4 internal un 963 964 // std::cout << "----- Calling PostStepDoIt 965 966 aParticleChange.Initialize(track) ; 967 const G4MaterialCutsCouple* couple = track.G 968 const G4DynamicParticle* aParticle = track.G 969 970 // Some kinematics 971 972 G4ParticleDefinition* definition = track.Get 973 G4double mass = definition->GetPDGMass(); 974 G4double kineticEnergy = aParticle->GetKinet 975 G4double totalEnergy = kineticEnergy + mass 976 G4double pSquare = kineticEnergy *( totalEne 977 G4double eSquare = totalEnergy * totalEnergy 978 G4double betaSquare = pSquare / eSquare; 979 G4ThreeVector particleDirection = aParticle- 980 981 G4double gamma = kineticEnergy / mass + 1.; 982 G4double r = electron_mass_c2 / mass; 983 G4double tMax = 2. * electron_mass_c2 *(gamm 984 985 // Validity range for delta electron cross s 986 G4double deltaCut = cutForDelta[couple->GetI 987 988 // This should not be a case 989 if (deltaCut >= tMax) 990 return G4VContinuousDiscreteProcess::PostS 991 992 G4double xc = deltaCut / tMax; 993 G4double rate = tMax / totalEnergy; 994 rate = rate*rate ; 995 G4double spin = aParticle->GetDefinition()-> 996 997 // Sampling follows ... 998 G4double x = 0.; 999 G4double gRej = 0.; 1000 1001 do { 1002 x = xc / (1. - (1. - xc) * G4UniformRand( 1003 1004 if (0.0 == spin) 1005 { 1006 gRej = 1.0 - betaSquare * x ; 1007 } 1008 else if (0.5 == spin) 1009 { 1010 gRej = (1. - betaSquare * x + 0.5 * x*x * r 1011 } 1012 else 1013 { 1014 gRej = (1. - betaSquare * x ) * (1. + x/(3. 1015 x*x * rate * (1. + 0.5*x/xc) / 3.0 / 1016 (1. + 1./(3.*xc) + rate *(1.+ 0.5/xc) / 3 1017 } 1018 1019 } while( G4UniformRand() > gRej ); 1020 1021 G4double deltaKineticEnergy = x * tMax; 1022 G4double deltaTotalMomentum = std::sqrt(del 1023 (deltaKineticEnergy + 2. * electr 1024 G4double totalMomentum = std::sqrt(pSquare) 1025 G4double cosTheta = deltaKineticEnergy * (t 1026 1027 // protection against cosTheta > 1 or < -1 1028 if ( cosTheta < -1. ) cosTheta = -1.; 1029 if ( cosTheta > 1. ) cosTheta = 1.; 1030 1031 // direction of the delta electron 1032 G4double phi = twopi * G4UniformRand(); 1033 G4double sinTheta = std::sqrt(1. - cosTheta 1034 G4double dirX = sinTheta * std::cos(phi); 1035 G4double dirY = sinTheta * std::sin(phi); 1036 G4double dirZ = cosTheta; 1037 1038 G4ThreeVector deltaDirection(dirX,dirY,dirZ 1039 deltaDirection.rotateUz(particleDirection); 1040 1041 // create G4DynamicParticle object for delt 1042 G4DynamicParticle* deltaRay = new G4Dynamic 1043 deltaRay->SetKineticEnergy( deltaKineticEne 1044 deltaRay->SetMomentumDirection(deltaDirecti 1045 deltaDirection.y(), 1046 deltaDirection.z()); 1047 deltaRay->SetDefinition(G4Electron::Electro 1048 1049 // fill aParticleChange 1050 G4double finalKineticEnergy = kineticEnergy 1051 std::size_t totalNumber = 1; 1052 1053 // Atomic relaxation 1054 1055 // ---- MGP ---- Temporary limitation: curr 1056 1057 std::size_t nSecondaries = 0; 1058 std::vector<G4DynamicParticle*>* secondaryV 1059 1060 if (definition == G4Proton::ProtonDefinitio 1061 { 1062 const G4Material* material = couple->Ge 1063 1064 // Lazy initialization of pixeCrossSect 1065 if (pixeCrossSectionHandler == 0) 1066 { 1067 // Instantiate pixeCrossSectionHandler wi 1068 // Ownership of interpolation is transfer 1069 G4IInterpolator* interpolation = new G4Lo 1070 pixeCrossSectionHandler = 1071 new G4PixeCrossSectionHandler(interpola 1072 G4String fileName("proton"); 1073 pixeCrossSectionHandler->LoadShellData(fi 1074 // pixeCrossSectionHandler->PrintData( 1075 } 1076 1077 // Select an atom in the current materi 1078 G4int Z = pixeCrossSectionHandler->Sele 1079 // std::cout << "G4hImpactIonisati 1080 1081 // G4double microscopicCross = Mic 1082 // kineticEnergy 1083 // Z, deltaCut); 1084 // G4double crossFromShells = pixeCrossS 1085 1086 //std::cout << "G4hImpactIonisation: Z= 1087 // << Z 1088 // << ", energy = " 1089 // << kineticEnergy/MeV 1090 // <<" MeV, microscopic = " 1091 // << microscopicCross/barn 1092 // << " barn, from shells = " 1093 // << crossFromShells/barn 1094 // << " barn" 1095 // << std::endl; 1096 1097 // Select a shell in the target atom ba 1098 G4int shellIndex = pixeCrossSectionHand 1099 1100 G4AtomicTransitionManager* transitionMa 1101 const G4AtomicShell* atomicShell = tran 1102 G4double bindingEnergy = atomicShell->B 1103 1104 // if (verboseLevel > 1) 1105 // { 1106 // G4cout << "G4hImpactIonisation::Po 1107 // << Z 1108 // << ", shell = " 1109 // << shellIndex 1110 // << ", bindingE (keV) = " 1111 // << bindingEnergy/keV 1112 // << G4endl; 1113 // } 1114 1115 // Generate PIXE if binding energy larg 1116 1117 G4ParticleDefinition* type = 0; 1118 1119 if (finalKineticEnergy >= bindingEnergy 1120 // && (bindingEnergy >= minGammaEnergy | 1121 { 1122 // Vacancy in subshell shellIndex; shellI 1123 G4int shellId = atomicShell->ShellId(); 1124 // Atomic relaxation: generate secondarie 1125 secondaryVector = atomicDeexcitation.Gene 1126 1127 // ---- Debug ---- 1128 //std::cout << "ShellId = " 1129 // <<shellId << " ---- Atomic relaxati 1130 // << secondaryVector->size() 1131 // << std::endl; 1132 1133 // ---- End debug --- 1134 1135 if (secondaryVector != 0) 1136 { 1137 nSecondaries = secondaryVector->size( 1138 for (std::size_t i = 0; i<nSecondarie 1139 { 1140 G4DynamicParticle* aSecondary = (*secon 1141 if (aSecondary) 1142 { 1143 G4double e = aSecondary->GetKinetic 1144 type = aSecondary->GetDefinition(); 1145 1146 // ---- Debug ---- 1147 //if (type == G4Gamma::GammaDefinit 1148 // { 1149 // std::cout << "Z = " << Z 1150 // << ", shell: " << shellId 1151 // << ", PIXE photon energy 1152 // << std::endl; 1153 // } 1154 // ---- End debug --- 1155 1156 if (e < finalKineticEnergy && 1157 ((type == G4Gamma::Gamma() && e > min 1158 (type == G4Electron::Electron() && e 1159 { 1160 // Subtract the energy of the emitted 1161 finalKineticEnergy -= e; 1162 totalNumber++; 1163 // ---- Debug ---- 1164 //if (type == G4Gamma::GammaDefinitio 1165 // { 1166 // std::cout << "Z = " << Z 1167 // << ", shell: " << shellId 1168 // << ", PIXE photon energy (k 1169 // << std::endl; 1170 // } 1171 // ---- End debug --- 1172 } 1173 else 1174 { 1175 // The atomic relaxation product has 1176 // ---- Debug ---- 1177 // if (type == G4Gamma::GammaDefiniti 1178 // { 1179 // std::cout << "Z = " << Z 1180 // 1181 // << ", PIXE photon energy = 1182 // << " keV below threshold 1183 // << std::endl; 1184 // } 1185 // ---- End debug --- 1186 1187 delete aSecondary; 1188 (*secondaryVector)[i] = 0; 1189 } 1190 } 1191 } 1192 } 1193 } 1194 } 1195 1196 1197 // Save delta-electrons 1198 1199 G4double eDeposit = 0.; 1200 1201 if (finalKineticEnergy > MinKineticEnergy) 1202 { 1203 G4double finalPx = totalMomentum*partic 1204 G4double finalPy = totalMomentum*partic 1205 G4double finalPz = totalMomentum*partic 1206 G4double finalMomentum = std::sqrt(fina 1207 finalPx /= finalMomentum; 1208 finalPy /= finalMomentum; 1209 finalPz /= finalMomentum; 1210 1211 aParticleChange.ProposeMomentumDirectio 1212 } 1213 else 1214 { 1215 eDeposit = finalKineticEnergy; 1216 finalKineticEnergy = 0.; 1217 aParticleChange.ProposeMomentumDirectio 1218 particleDirection.y(), 1219 particleDirection.z()); 1220 if(!aParticle->GetDefinition()->GetProc 1221 GetAtRestProcessVector()->size()) 1222 aParticleChange.ProposeTrackStatus(fS 1223 else 1224 aParticleChange.ProposeTrackStatus(fS 1225 } 1226 1227 aParticleChange.ProposeEnergy(finalKineticE 1228 aParticleChange.ProposeLocalEnergyDeposit ( 1229 aParticleChange.SetNumberOfSecondaries((G4i 1230 aParticleChange.AddSecondary(deltaRay); 1231 1232 // ---- Debug ---- 1233 // std::cout << "RDHadronIonisation - fina 1234 // << finalKineticEnergy/MeV 1235 // << ", delta KineticEnergy (keV) = " 1236 // << deltaKineticEnergy/keV 1237 // << ", energy deposit (MeV) = " 1238 // << eDeposit/MeV 1239 // << std::endl; 1240 // ---- End debug --- 1241 1242 // Save Fluorescence and Auger 1243 1244 if (secondaryVector != 0) 1245 { 1246 for (std::size_t l = 0; l < nSecondarie 1247 { 1248 G4DynamicParticle* secondary = (*secondar 1249 if (secondary) aParticleChange.AddSeconda 1250 1251 // ---- Debug ---- 1252 //if (secondary != 0) 1253 // { 1254 // if (secondary->GetDefinition() == G4 1255 // { 1256 // G4double eX = secondary->GetKinetic 1257 // std::cout << " PIXE photon of energ 1258 // << " keV added to ParticleChang 1259 // << std::endl; 1260 // } 1261 //} 1262 // ---- End debug --- 1263 1264 } 1265 delete secondaryVector; 1266 } 1267 1268 return G4VContinuousDiscreteProcess::PostSt 1269 } 1270 1271 // ------------------------------------------ 1272 1273 G4double G4hImpactIonisation::ComputeDEDX(con 1274 const G4MaterialCutsCouple* coupl 1275 G4double kineticEnergy) 1276 { 1277 const G4Material* material = couple->GetMat 1278 G4Proton* proton = G4Proton::Proton(); 1279 G4AntiProton* antiproton = G4AntiProton::An 1280 G4double dedx = 0.; 1281 1282 G4double tScaled = kineticEnergy * proton_m 1283 charge = aParticle->GetPDGCharge() ; 1284 1285 if( charge > 0.) 1286 { 1287 if (tScaled > protonHighEnergy) 1288 { 1289 dedx = G4EnergyLossTables::GetDEDX(proton 1290 } 1291 else 1292 { 1293 dedx = ProtonParametrisedDEDX(couple,tSca 1294 } 1295 } 1296 else 1297 { 1298 if (tScaled > antiprotonHighEnergy) 1299 { 1300 dedx = G4EnergyLossTables::GetDEDX(antipr 1301 } 1302 else 1303 { 1304 dedx = AntiProtonParametrisedDEDX(couple, 1305 } 1306 } 1307 dedx *= theIonEffChargeModel->TheValue(aPar 1308 1309 return dedx ; 1310 } 1311 1312 1313 G4double G4hImpactIonisation::BarkasTerm(cons 1314 G4double kineticEnergy) const 1315 //Function to compute the Barkas term for pro 1316 // 1317 //Ref. Z_1^3 effect in the stopping power of 1318 // J.C Ashley and R.H.Ritchie 1319 // Physical review B Vol.5 No.7 1 April 1 1320 // 1321 { 1322 static G4ThreadLocal G4double FTable[47][2] 1323 { 0.02, 21.5}, 1324 { 0.03, 20.0}, 1325 { 0.04, 18.0}, 1326 { 0.05, 15.6}, 1327 { 0.06, 15.0}, 1328 { 0.07, 14.0}, 1329 { 0.08, 13.5}, 1330 { 0.09, 13.}, 1331 { 0.1, 12.2}, 1332 { 0.2, 9.25}, 1333 { 0.3, 7.0}, 1334 { 0.4, 6.0}, 1335 { 0.5, 4.5}, 1336 { 0.6, 3.5}, 1337 { 0.7, 3.0}, 1338 { 0.8, 2.5}, 1339 { 0.9, 2.0}, 1340 { 1.0, 1.7}, 1341 { 1.2, 1.2}, 1342 { 1.3, 1.0}, 1343 { 1.4, 0.86}, 1344 { 1.5, 0.7}, 1345 { 1.6, 0.61}, 1346 { 1.7, 0.52}, 1347 { 1.8, 0.5}, 1348 { 1.9, 0.43}, 1349 { 2.0, 0.42}, 1350 { 2.1, 0.3}, 1351 { 2.4, 0.2}, 1352 { 3.0, 0.13}, 1353 { 3.08, 0.1}, 1354 { 3.1, 0.09}, 1355 { 3.3, 0.08}, 1356 { 3.5, 0.07}, 1357 { 3.8, 0.06}, 1358 { 4.0, 0.051}, 1359 { 4.1, 0.04}, 1360 { 4.8, 0.03}, 1361 { 5.0, 0.024}, 1362 { 5.1, 0.02}, 1363 { 6.0, 0.013}, 1364 { 6.5, 0.01}, 1365 { 7.0, 0.009}, 1366 { 7.1, 0.008}, 1367 { 8.0, 0.006}, 1368 { 9.0, 0.0032}, 1369 { 10.0, 0.0025} }; 1370 1371 // Information on particle and material 1372 G4double kinE = kineticEnergy ; 1373 if(0.5*MeV > kinE) kinE = 0.5*MeV ; 1374 G4double gamma = 1.0 + kinE / proton_mass_c 1375 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ; 1376 if(0.0 >= beta2) return 0.0; 1377 1378 G4double BTerm = 0.0; 1379 //G4double AMaterial = 0.0; 1380 G4double ZMaterial = 0.0; 1381 const G4ElementVector* theElementVector = m 1382 G4int numberOfElements = (G4int)material->G 1383 1384 for (G4int i = 0; i<numberOfElements; i++) 1385 1386 //AMaterial = (*theElementVector)[i]->Get 1387 ZMaterial = (*theElementVector)[i]->GetZ( 1388 1389 G4double X = 137.0 * 137.0 * beta2 / ZMat 1390 1391 // Variables to compute L_1 1392 G4double Eta0Chi = 0.8; 1393 G4double EtaChi = Eta0Chi * ( 1.0 + 6.02* 1394 G4double W = ( EtaChi * std::pow( ZMateri 1395 G4double FunctionOfW = FTable[46][1]*FTab 1396 1397 for(G4int j=0; j<47; j++) { 1398 1399 if( W < FTable[j][0] ) { 1400 1401 if(0 == j) { 1402 FunctionOfW = FTable[0][1] ; 1403 1404 } else { 1405 FunctionOfW = (FTable[j][1] - FTabl 1406 / (FTable[j][0] - FTable[j-1][0]) 1407 + FTable[j-1][1] ; 1408 } 1409 1410 break; 1411 } 1412 1413 } 1414 1415 BTerm += FunctionOfW /( std::sqrt(ZMateri 1416 } 1417 1418 BTerm *= twopi_mc2_rcl2 * (material->GetEle 1419 1420 return BTerm; 1421 } 1422 1423 1424 1425 G4double G4hImpactIonisation::BlochTerm(const 1426 G4double kineticEnergy, 1427 G4double cSquare) const 1428 //Function to compute the Bloch term for prot 1429 // 1430 //Ref. Z_1^3 effect in the stopping power of 1431 // J.C Ashley and R.H.Ritchie 1432 // Physical review B Vol.5 No.7 1 April 1 1433 // 1434 { 1435 G4double eLoss = 0.0 ; 1436 G4double gamma = 1.0 + kineticEnergy / prot 1437 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ; 1438 G4double y = cSquare / (137.0*137.0*beta2) 1439 1440 if(y < 0.05) { 1441 eLoss = 1.202 ; 1442 1443 } else { 1444 eLoss = 1.0 / (1.0 + y) ; 1445 G4double de = eLoss ; 1446 1447 for(G4int i=2; de>eLoss*0.01; i++) { 1448 de = 1.0/( i * (i*i + y)) ; 1449 eLoss += de ; 1450 } 1451 } 1452 eLoss *= -1.0 * y * cSquare * twopi_mc2_rcl 1453 (material->GetElectronDensity()) / beta2 1454 1455 return eLoss; 1456 } 1457 1458 1459 1460 G4double G4hImpactIonisation::ElectronicLossF 1461 const G4DynamicParticle* partic 1462 const G4MaterialCutsCouple* cou 1463 G4double meanLoss, 1464 G4double step) const 1465 // calculate actual loss from the mean loss 1466 // The model used to get the fluctuation is 1467 // as in Glandz in Geant3. 1468 { 1469 // data members to speed up the fluctuation 1470 // G4int imat ; 1471 // G4double f1Fluct,f2Fluct,e1Fluct,e2Fluc 1472 // G4double e1LogFluct,e2LogFluct,ipotLogF 1473 1474 static const G4double minLoss = 1.*eV ; 1475 static const G4double kappa = 10. ; 1476 static const G4double theBohrBeta2 = 50.0 * 1477 1478 const G4Material* material = couple->GetMat 1479 G4int imaterial = couple->GetIndex() ; 1480 G4double ipotFluct = material->GetIonisat 1481 G4double electronDensity = material->GetEle 1482 G4double zeff = electronDensity/(material-> 1483 1484 // get particle data 1485 G4double tkin = particle->GetKineticEnerg 1486 G4double particleMass = particle->GetMass() 1487 G4double deltaCutInKineticEnergyNow = cutFo 1488 1489 // shortcut for very very small loss 1490 if(meanLoss < minLoss) return meanLoss ; 1491 1492 // Validity range for delta electron cross 1493 G4double threshold = std::max(deltaCutInKin 1494 G4double loss, siga; 1495 1496 G4double rmass = electron_mass_c2/particleM 1497 G4double tau = tkin/particleMass; 1498 G4double tau1 = tau+1.0; 1499 G4double tau2 = tau*(tau+2.); 1500 G4double tMax = 2.*electron_mass_c2*tau2 1501 1502 1503 if(tMax > threshold) tMax = threshold; 1504 G4double beta2 = tau2/(tau1*tau1); 1505 1506 // Gaussian fluctuation 1507 if(meanLoss > kappa*tMax || tMax < kappa*ip 1508 { 1509 siga = tMax * (1.0-0.5*beta2) * step * 1510 * electronDensity / beta2 ; 1511 1512 // High velocity or negatively charged 1513 if( beta2 > 3.0*theBohrBeta2*zeff || ch 1514 siga = std::sqrt( siga * chargeSquare ) ; 1515 1516 // Low velocity - additional ion charge flu 1517 // Q.Yang et al., NIM B61(1991)149-155. 1518 } else { 1519 G4double chu = theIonChuFluctuationModel->T 1520 G4double yang = theIonYangFluctuationModel- 1521 siga = std::sqrt( siga * (chargeSquare * ch 1522 } 1523 1524 do { 1525 loss = G4RandGauss::shoot(meanLoss,si 1526 } while (loss < 0. || loss > 2.0*meanLo 1527 return loss; 1528 } 1529 1530 // Non Gaussian fluctuation 1531 static const G4double probLim = 0.01 ; 1532 static const G4double sumaLim = -std::log(p 1533 static const G4double alim = 10.; 1534 1535 G4double suma,w1,w2,C,e0,lossc,w; 1536 G4double a1,a2,a3; 1537 G4int p1,p2,p3; 1538 G4int nb; 1539 G4double corrfac, na,alfa,rfac,namean,sa,al 1540 G4double dp3; 1541 1542 G4double f1Fluct = material->GetIonisat 1543 G4double f2Fluct = material->GetIonisat 1544 G4double e1Fluct = material->GetIonisat 1545 G4double e2Fluct = material->GetIonisat 1546 G4double e1LogFluct = material->GetIonisat 1547 G4double e2LogFluct = material->GetIonisat 1548 G4double rateFluct = material->GetIonisat 1549 G4double ipotLogFluct= material->GetIonisat 1550 1551 w1 = tMax/ipotFluct; 1552 w2 = std::log(2.*electron_mass_c2*tau2); 1553 1554 C = meanLoss*(1.-rateFluct)/(w2-ipotLogFluc 1555 1556 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluc 1557 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluc 1558 a3 = rateFluct*meanLoss*(tMax-ipotFluct)/(i 1559 if(a1 < 0.0) a1 = 0.0; 1560 if(a2 < 0.0) a2 = 0.0; 1561 if(a3 < 0.0) a3 = 0.0; 1562 1563 suma = a1+a2+a3; 1564 1565 loss = 0.; 1566 1567 1568 if(suma < sumaLim) // very smal 1569 { 1570 e0 = material->GetIonisation()->GetEner 1571 1572 if(tMax == ipotFluct) 1573 { 1574 a3 = meanLoss/e0; 1575 1576 if(a3>alim) 1577 { 1578 siga=std::sqrt(a3) ; 1579 p3 = std::max(0,G4int(G4RandGauss::sh 1580 } 1581 else 1582 p3 = (G4int)G4Poisson(a3); 1583 1584 loss = p3*e0 ; 1585 1586 if(p3 > 0) 1587 loss += (1.-2.*G4UniformRand())*e0 ; 1588 1589 } 1590 else 1591 { 1592 tMax = tMax-ipotFluct+e0 ; 1593 a3 = meanLoss*(tMax-e0)/(tMax*e0*std::log 1594 1595 if(a3>alim) 1596 { 1597 siga=std::sqrt(a3) ; 1598 p3 = std::max(0,int(G4RandGauss::shoo 1599 } 1600 else 1601 p3 = (G4int)G4Poisson(a3); 1602 1603 if(p3 > 0) 1604 { 1605 w = (tMax-e0)/tMax ; 1606 if(p3 > nmaxCont2) 1607 { 1608 dp3 = G4float(p3) ; 1609 corrfac = dp3/G4float(nmaxCont2) ; 1610 p3 = G4int(nmaxCont2) ; 1611 } 1612 else 1613 corrfac = 1. ; 1614 1615 for(G4int i=0; i<p3; i++) loss += 1./ 1616 loss *= e0*corrfac ; 1617 } 1618 } 1619 } 1620 1621 else // not so 1622 { 1623 // excitation type 1 1624 if(a1>alim) 1625 { 1626 siga=std::sqrt(a1) ; 1627 p1 = std::max(0,G4int(G4RandGauss::shoot( 1628 } 1629 else 1630 p1 = (G4int)G4Poisson(a1); 1631 1632 // excitation type 2 1633 if(a2>alim) 1634 { 1635 siga=std::sqrt(a2) ; 1636 p2 = std::max(0,G4int(G4RandGauss::shoot( 1637 } 1638 else 1639 p2 = (G4int)G4Poisson(a2); 1640 1641 loss = p1*e1Fluct+p2*e2Fluct; 1642 1643 // smearing to avoid unphysical peaks 1644 if(p2 > 0) 1645 loss += (1.-2.*G4UniformRand())*e2Flu 1646 else if (loss>0.) 1647 loss += (1.-2.*G4UniformRand())*e1Flu 1648 1649 // ionisation ......................... 1650 if(a3 > 0.) 1651 { 1652 if(a3>alim) 1653 { 1654 siga=std::sqrt(a3) ; 1655 p3 = std::max(0,G4int(G4RandGauss::sh 1656 } 1657 else 1658 p3 = (G4int)G4Poisson(a3); 1659 1660 lossc = 0.; 1661 if(p3 > 0) 1662 { 1663 na = 0.; 1664 alfa = 1.; 1665 if (p3 > nmaxCont2) 1666 { 1667 dp3 = G4float(p3); 1668 rfac = dp3/(G4float(nmaxCont2)+dp 1669 namean = G4float(p3)*rfac; 1670 sa = G4float(nmaxCont1)*rfac; 1671 na = G4RandGauss::shoot(namean, 1672 if (na > 0.) 1673 { 1674 alfa = w1*G4float(nmaxCont2+p3)/ 1675 (w1*G4float(nmaxCont2)+G4float(p3)); 1676 alfa1 = alfa*std::log(alfa)/(alfa- 1677 ea = na*ipotFluct*alfa1; 1678 sea = ipotFluct*std::sqrt(na*(al 1679 lossc += G4RandGauss::shoot(ea,sea) 1680 } 1681 } 1682 1683 nb = G4int(G4float(p3)-na); 1684 if (nb > 0) 1685 { 1686 w2 = alfa*ipotFluct; 1687 w = (tMax-w2)/tMax; 1688 for (G4int k=0; k<nb; k++) lossc += w2/ 1689 } 1690 } 1691 loss += lossc; 1692 } 1693 } 1694 1695 return loss ; 1696 } 1697 1698 1699 1700 void G4hImpactIonisation::SetCutForSecondaryP 1701 { 1702 minGammaEnergy = cut; 1703 } 1704 1705 1706 1707 void G4hImpactIonisation::SetCutForAugerElect 1708 { 1709 minElectronEnergy = cut; 1710 } 1711 1712 1713 1714 void G4hImpactIonisation::ActivateAugerElectr 1715 { 1716 atomicDeexcitation.ActivateAugerElectronPro 1717 } 1718 1719 1720 1721 void G4hImpactIonisation::PrintInfoDefinition 1722 { 1723 G4String comments = " Knock-on electron cr 1724 comments += "\n Good description abo 1725 comments += " Delta ray energy sampl 1726 1727 G4cout << G4endl << GetProcessName() << ": 1728 << "\n PhysicsTables from " < 1729 << " to " << HighestKineticEnergy / 1730 << " in " << TotBin << " bins." 1731 << "\n Electronic stopping power mo 1732 << protonTable 1733 << "\n from " << protonLowEne 1734 << " to " << protonHighEnergy / MeV 1735 G4cout << "\n Parametrisation model 1736 << antiprotonTable 1737 << "\n from " << antiprotonLo 1738 << " to " << antiprotonHighEnergy / 1739 if(theBarkas){ 1740 G4cout << " Parametrization of the 1741 << G4endl ; 1742 } 1743 if(nStopping) { 1744 G4cout << " Nuclear stopping power 1745 << G4endl ; 1746 } 1747 1748 G4bool printHead = true; 1749 1750 const G4ProductionCutsTable* theCoupleTable 1751 G4ProductionCutsTable::GetProductionCutsT 1752 G4int numOfCouples = (G4int)theCoupleTable- 1753 1754 // loop for materials 1755 1756 for (G4int j=0 ; j < numOfCouples; ++j) { 1757 1758 const G4MaterialCutsCouple* couple = theC 1759 const G4Material* material= couple->GetMa 1760 G4double deltaCutNow = cutForDelta[(coupl 1761 G4double excitationEnergy = material->Get 1762 1763 if(excitationEnergy > deltaCutNow) { 1764 if(printHead) { 1765 printHead = false ; 1766 1767 G4cout << " material 1768 G4cout << G4endl; 1769 } 1770 1771 G4cout << std::setw(20) << material->Ge 1772 << std::setw(15) << excitationEnergy/k 1773 } 1774 } 1775 } 1776