Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 27 #include "G4AdjointeIonisationModel.hh" 28 29 #include "G4AdjointElectron.hh" 30 #include "G4Electron.hh" 31 #include "G4ParticleChange.hh" 32 #include "G4PhysicalConstants.hh" 33 #include "G4TrackStatus.hh" 34 35 //////////////////////////////////////////////////////////////////////////////// 36 G4AdjointeIonisationModel::G4AdjointeIonisationModel() 37 : G4VEmAdjointModel("Inv_eIon_model") 38 39 { 40 fUseMatrix = true; 41 fUseMatrixPerElement = true; 42 fApplyCutInRange = true; 43 fOneMatrixForAllElements = true; 44 45 fAdjEquivDirectPrimPart = G4AdjointElectron::AdjointElectron(); 46 fAdjEquivDirectSecondPart = G4AdjointElectron::AdjointElectron(); 47 fDirectPrimaryPart = G4Electron::Electron(); 48 fSecondPartSameType = true; 49 } 50 51 //////////////////////////////////////////////////////////////////////////////// 52 G4AdjointeIonisationModel::~G4AdjointeIonisationModel() {} 53 54 //////////////////////////////////////////////////////////////////////////////// 55 void G4AdjointeIonisationModel::SampleSecondaries( 56 const G4Track& aTrack, G4bool IsScatProjToProj, 57 G4ParticleChange* fParticleChange) 58 { 59 const G4DynamicParticle* theAdjointPrimary = aTrack.GetDynamicParticle(); 60 61 // Elastic inverse scattering 62 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); 63 G4double adjointPrimP = theAdjointPrimary->GetTotalMomentum(); 64 65 if(adjointPrimKinEnergy > GetHighEnergyLimit() * 0.999) 66 { 67 return; 68 } 69 70 // Sample secondary energy 71 G4double projectileKinEnergy; 72 if(!fWithRapidSampling) 73 { // used by default 74 projectileKinEnergy = 75 SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProj); 76 77 // Caution!!! this weight correction should be always applied 78 CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), 79 adjointPrimKinEnergy, projectileKinEnergy, 80 IsScatProjToProj); 81 } 82 else 83 { // only for testing 84 G4double Emin, Emax; 85 if(IsScatProjToProj) 86 { 87 Emin = GetSecondAdjEnergyMinForScatProjToProj(adjointPrimKinEnergy, 88 fTcutSecond); 89 Emax = GetSecondAdjEnergyMaxForScatProjToProj(adjointPrimKinEnergy); 90 } 91 else 92 { 93 Emin = GetSecondAdjEnergyMinForProdToProj(adjointPrimKinEnergy); 94 Emax = GetSecondAdjEnergyMaxForProdToProj(adjointPrimKinEnergy); 95 } 96 projectileKinEnergy = Emin * std::pow(Emax / Emin, G4UniformRand()); 97 98 fLastCS = fLastAdjointCSForScatProjToProj; 99 if(!IsScatProjToProj) 100 fLastCS = fLastAdjointCSForProdToProj; 101 102 G4double new_weight = aTrack.GetWeight(); 103 G4double used_diffCS = 104 fLastCS * std::log(Emax / Emin) / projectileKinEnergy; 105 G4double needed_diffCS = adjointPrimKinEnergy / projectileKinEnergy; 106 if(!IsScatProjToProj) 107 needed_diffCS *= DiffCrossSectionPerVolumePrimToSecond( 108 fCurrentMaterial, projectileKinEnergy, adjointPrimKinEnergy); 109 else 110 needed_diffCS *= DiffCrossSectionPerVolumePrimToScatPrim( 111 fCurrentMaterial, projectileKinEnergy, adjointPrimKinEnergy); 112 new_weight *= needed_diffCS / used_diffCS; 113 fParticleChange->SetParentWeightByProcess(false); 114 fParticleChange->SetSecondaryWeightByProcess(false); 115 fParticleChange->ProposeParentWeight(new_weight); 116 } 117 118 // Kinematic: 119 // we consider a two body elastic scattering for the forward processes where 120 // the projectile knock on an e- at rest and gives it part of its energy 121 G4double projectileM0 = fAdjEquivDirectPrimPart->GetPDGMass(); 122 G4double projectileTotalEnergy = projectileM0 + projectileKinEnergy; 123 G4double projectileP2 = 124 projectileTotalEnergy * projectileTotalEnergy - projectileM0 * projectileM0; 125 126 // Companion 127 G4double companionM0 = fAdjEquivDirectPrimPart->GetPDGMass(); 128 if(IsScatProjToProj) 129 { 130 companionM0 = fAdjEquivDirectSecondPart->GetPDGMass(); 131 } 132 G4double companionTotalEnergy = 133 companionM0 + projectileKinEnergy - adjointPrimKinEnergy; 134 G4double companionP2 = 135 companionTotalEnergy * companionTotalEnergy - companionM0 * companionM0; 136 137 // Projectile momentum 138 G4double P_parallel = 139 (adjointPrimP * adjointPrimP + projectileP2 - companionP2) / 140 (2. * adjointPrimP); 141 G4double P_perp = std::sqrt(projectileP2 - P_parallel * P_parallel); 142 G4ThreeVector dir_parallel = theAdjointPrimary->GetMomentumDirection(); 143 G4double phi = G4UniformRand() * twopi; 144 G4ThreeVector projectileMomentum = 145 G4ThreeVector(P_perp * std::cos(phi), P_perp * std::sin(phi), P_parallel); 146 projectileMomentum.rotateUz(dir_parallel); 147 148 if(!IsScatProjToProj) 149 { // kill the primary and add a secondary 150 fParticleChange->ProposeTrackStatus(fStopAndKill); 151 fParticleChange->AddSecondary( 152 new G4DynamicParticle(fAdjEquivDirectPrimPart, projectileMomentum)); 153 } 154 else 155 { 156 fParticleChange->ProposeEnergy(projectileKinEnergy); 157 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit()); 158 } 159 } 160 161 //////////////////////////////////////////////////////////////////////////////// 162 // The implementation here is correct for energy loss process, for the 163 // photoelectric and compton scattering the method should be redefined 164 G4double G4AdjointeIonisationModel::DiffCrossSectionPerAtomPrimToSecond( 165 G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double) 166 { 167 G4double dSigmadEprod = 0.; 168 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProj(kinEnergyProd); 169 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProj(kinEnergyProd); 170 171 // the produced particle should have a kinetic energy smaller than the 172 // projectile 173 if(kinEnergyProj > Emin_proj && kinEnergyProj <= Emax_proj) 174 { 175 dSigmadEprod = Z * DiffCrossSectionMoller(kinEnergyProj, kinEnergyProd); 176 } 177 return dSigmadEprod; 178 } 179 180 ////////////////////////////////////////////////////////////////////////////// 181 G4double G4AdjointeIonisationModel::DiffCrossSectionMoller( 182 G4double kinEnergyProj, G4double kinEnergyProd) 183 { 184 // G4double energy = kinEnergyProj + electron_mass_c2; 185 G4double x = kinEnergyProd / kinEnergyProj; 186 G4double gam = (kinEnergyProj + electron_mass_c2) / electron_mass_c2; 187 G4double gamma2 = gam * gam; 188 G4double beta2 = 1.0 - 1.0 / gamma2; 189 190 G4double gg = (2.0 * gam - 1.0) / gamma2; 191 G4double y = 1.0 - x; 192 G4double fac = twopi_mc2_rcl2 / electron_mass_c2; 193 G4double dCS = 194 fac * (1. - gg + ((1.0 - gg * x) / (x * x)) + ((1.0 - gg * y) / (y * y))) / 195 (beta2 * (gam - 1.)); 196 return dCS / kinEnergyProj; 197 } 198