Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 26 27 #include "G4AdjointIonIonisationModel.hh" 27 #include "G4AdjointIonIonisationModel.hh" 28 28 29 #include "G4AdjointCSManager.hh" 29 #include "G4AdjointCSManager.hh" 30 #include "G4AdjointElectron.hh" 30 #include "G4AdjointElectron.hh" 31 #include "G4BetheBlochModel.hh" 31 #include "G4BetheBlochModel.hh" 32 #include "G4BraggIonModel.hh" 32 #include "G4BraggIonModel.hh" 33 #include "G4GenericIon.hh" 33 #include "G4GenericIon.hh" 34 #include "G4NistManager.hh" 34 #include "G4NistManager.hh" 35 #include "G4ParticleChange.hh" 35 #include "G4ParticleChange.hh" 36 #include "G4PhysicalConstants.hh" 36 #include "G4PhysicalConstants.hh" 37 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh" 38 #include "G4TrackStatus.hh" 38 #include "G4TrackStatus.hh" 39 39 40 ////////////////////////////////////////////// 40 //////////////////////////////////////////////////////////////////////////////// 41 G4AdjointIonIonisationModel::G4AdjointIonIonis 41 G4AdjointIonIonisationModel::G4AdjointIonIonisationModel() 42 : G4VEmAdjointModel("Adjoint_IonIonisation") 42 : G4VEmAdjointModel("Adjoint_IonIonisation") 43 { 43 { 44 fUseMatrix = true; 44 fUseMatrix = true; 45 fUseMatrixPerElement = true; 45 fUseMatrixPerElement = true; 46 fApplyCutInRange = true; 46 fApplyCutInRange = true; 47 fOneMatrixForAllElements = true; 47 fOneMatrixForAllElements = true; 48 fSecondPartSameType = false; 48 fSecondPartSameType = false; 49 49 50 // The direct EM Model is taken as BetheBloc 50 // The direct EM Model is taken as BetheBloch. It is only used for the 51 // computation of the differential cross sec 51 // computation of the differential cross section. 52 // The Bragg model could be used as an alter 52 // The Bragg model could be used as an alternative as it offers the same 53 // differential cross section 53 // differential cross section 54 54 55 fBetheBlochDirectEMModel = new G4BetheBloch 55 fBetheBlochDirectEMModel = new G4BetheBlochModel(G4GenericIon::GenericIon()); 56 fBraggIonDirectEMModel = new G4BraggIonMo 56 fBraggIonDirectEMModel = new G4BraggIonModel(G4GenericIon::GenericIon()); 57 fAdjEquivDirectSecondPart = G4AdjointElectro 57 fAdjEquivDirectSecondPart = G4AdjointElectron::AdjointElectron(); 58 fDirectPrimaryPart = nullptr; 58 fDirectPrimaryPart = nullptr; 59 } 59 } 60 60 61 ////////////////////////////////////////////// 61 //////////////////////////////////////////////////////////////////////////////// 62 G4AdjointIonIonisationModel::~G4AdjointIonIoni 62 G4AdjointIonIonisationModel::~G4AdjointIonIonisationModel() {} 63 63 64 ////////////////////////////////////////////// 64 //////////////////////////////////////////////////////////////////////////////// 65 void G4AdjointIonIonisationModel::SampleSecond 65 void G4AdjointIonIonisationModel::SampleSecondaries( 66 const G4Track& aTrack, G4bool isScatProjToPr 66 const G4Track& aTrack, G4bool isScatProjToProj, 67 G4ParticleChange* fParticleChange) 67 G4ParticleChange* fParticleChange) 68 { 68 { 69 const G4DynamicParticle* theAdjointPrimary = 69 const G4DynamicParticle* theAdjointPrimary = aTrack.GetDynamicParticle(); 70 70 71 // Elastic inverse scattering 71 // Elastic inverse scattering 72 G4double adjointPrimKinEnergy = theAdjointPr 72 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); 73 G4double adjointPrimP = theAdjointPr 73 G4double adjointPrimP = theAdjointPrimary->GetTotalMomentum(); 74 74 75 if(adjointPrimKinEnergy > GetHighEnergyLimit 75 if(adjointPrimKinEnergy > GetHighEnergyLimit() * 0.999) 76 { 76 { 77 return; 77 return; 78 } 78 } 79 79 80 // Sample secondary energy 80 // Sample secondary energy 81 G4double projectileKinEnergy = 81 G4double projectileKinEnergy = 82 SampleAdjSecEnergyFromCSMatrix(adjointPrim 82 SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, isScatProjToProj); 83 // Caution !!!this weight correction should 83 // Caution !!!this weight correction should be always applied 84 CorrectPostStepWeight(fParticleChange, aTrac 84 CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), 85 adjointPrimKinEnergy, 85 adjointPrimKinEnergy, projectileKinEnergy, 86 isScatProjToProj); 86 isScatProjToProj); 87 87 88 // Kinematics: 88 // Kinematics: 89 // we consider a two body elastic scattering 89 // we consider a two body elastic scattering for the forward processes where 90 // the projectile knock on an e- at rest and 90 // the projectile knock on an e- at rest and gives it part of its energy 91 G4double projectileM0 = fAdjEquivDi 91 G4double projectileM0 = fAdjEquivDirectPrimPart->GetPDGMass(); 92 G4double projectileTotalEnergy = projectileM 92 G4double projectileTotalEnergy = projectileM0 + projectileKinEnergy; 93 G4double projectileP2 = 93 G4double projectileP2 = 94 projectileTotalEnergy * projectileTotalEne 94 projectileTotalEnergy * projectileTotalEnergy - projectileM0 * projectileM0; 95 95 96 // Companion 96 // Companion 97 G4double companionM0 = fAdjEquivDirectPrimPa 97 G4double companionM0 = fAdjEquivDirectPrimPart->GetPDGMass(); 98 if(isScatProjToProj) 98 if(isScatProjToProj) 99 { 99 { 100 companionM0 = fAdjEquivDirectSecondPart->G 100 companionM0 = fAdjEquivDirectSecondPart->GetPDGMass(); 101 } 101 } 102 G4double companionTotalEnergy = 102 G4double companionTotalEnergy = 103 companionM0 + projectileKinEnergy - adjoin 103 companionM0 + projectileKinEnergy - adjointPrimKinEnergy; 104 G4double companionP2 = 104 G4double companionP2 = 105 companionTotalEnergy * companionTotalEnerg 105 companionTotalEnergy * companionTotalEnergy - companionM0 * companionM0; 106 106 107 // Projectile momentum 107 // Projectile momentum 108 G4double P_parallel = 108 G4double P_parallel = 109 (adjointPrimP * adjointPrimP + projectileP 109 (adjointPrimP * adjointPrimP + projectileP2 - companionP2) / 110 (2. * adjointPrimP); 110 (2. * adjointPrimP); 111 G4double P_perp = std::sqrt(projectileP2 - P 111 G4double P_perp = std::sqrt(projectileP2 - P_parallel * P_parallel); 112 G4ThreeVector dir_parallel = theAdjointPrima 112 G4ThreeVector dir_parallel = theAdjointPrimary->GetMomentumDirection(); 113 G4double phi = G4UniformRand() 113 G4double phi = G4UniformRand() * twopi; 114 G4ThreeVector projectileMomentum = 114 G4ThreeVector projectileMomentum = 115 G4ThreeVector(P_perp * std::cos(phi), P_pe 115 G4ThreeVector(P_perp * std::cos(phi), P_perp * std::sin(phi), P_parallel); 116 projectileMomentum.rotateUz(dir_parallel); 116 projectileMomentum.rotateUz(dir_parallel); 117 117 118 if(!isScatProjToProj) 118 if(!isScatProjToProj) 119 { // kill the primary and add a secondary 119 { // kill the primary and add a secondary 120 fParticleChange->ProposeTrackStatus(fStopA 120 fParticleChange->ProposeTrackStatus(fStopAndKill); 121 fParticleChange->AddSecondary( 121 fParticleChange->AddSecondary( 122 new G4DynamicParticle(fAdjEquivDirectPri 122 new G4DynamicParticle(fAdjEquivDirectPrimPart, projectileMomentum)); 123 } 123 } 124 else 124 else 125 { 125 { 126 fParticleChange->ProposeEnergy(projectileK 126 fParticleChange->ProposeEnergy(projectileKinEnergy); 127 fParticleChange->ProposeMomentumDirection( 127 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit()); 128 } 128 } 129 } 129 } 130 130 131 ////////////////////////////////////////////// 131 //////////////////////////////////////////////////////////////////////////////// 132 G4double G4AdjointIonIonisationModel::DiffCros 132 G4double G4AdjointIonIonisationModel::DiffCrossSectionPerAtomPrimToSecond( 133 G4double kinEnergyProj, G4double kinEnergyPr 133 G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A) 134 { 134 { 135 // Probably that here the Bragg Model should 135 // Probably that here the Bragg Model should be also used for 136 // kinEnergyProj/nuc<2MeV 136 // kinEnergyProj/nuc<2MeV 137 G4double dSigmadEprod = 0.; 137 G4double dSigmadEprod = 0.; 138 G4double Emax_proj = GetSecondAdjEnergyMa 138 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProj(kinEnergyProd); 139 G4double Emin_proj = GetSecondAdjEnergyMi 139 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProj(kinEnergyProd); 140 140 141 G4double kinEnergyProjScaled = fMassRatio * 141 G4double kinEnergyProjScaled = fMassRatio * kinEnergyProj; 142 142 143 // the produced particle should have a kinet 143 // the produced particle should have a kinetic energy smaller than the 144 // projectile 144 // projectile 145 if(kinEnergyProj > Emin_proj && kinEnergyPro 145 if(kinEnergyProj > Emin_proj && kinEnergyProj <= Emax_proj) 146 { 146 { 147 G4double Tmax = kinEnergyProj; 147 G4double Tmax = kinEnergyProj; 148 148 149 G4double E1 = kinEnergyProd; 149 G4double E1 = kinEnergyProd; 150 G4double E2 = kinEnergyProd * 1.000001; 150 G4double E2 = kinEnergyProd * 1.000001; 151 G4double dE = (E2 - E1); 151 G4double dE = (E2 - E1); 152 G4double sigma1, sigma2; 152 G4double sigma1, sigma2; 153 fDirectModel = fBraggIonDirectEMModel; 153 fDirectModel = fBraggIonDirectEMModel; 154 if(kinEnergyProjScaled > 2. * MeV && !fUse 154 if(kinEnergyProjScaled > 2. * MeV && !fUseOnlyBragg) 155 fDirectModel = fBetheBlochDirectEMModel; 155 fDirectModel = fBetheBlochDirectEMModel; 156 sigma1 = fDirectModel->ComputeCrossSection 156 sigma1 = fDirectModel->ComputeCrossSectionPerAtom( 157 fDirectPrimaryPart, kinEnergyProj, Z, A, 157 fDirectPrimaryPart, kinEnergyProj, Z, A, E1, 1.e20); 158 sigma2 = fDirectModel->ComputeCrossSection 158 sigma2 = fDirectModel->ComputeCrossSectionPerAtom( 159 fDirectPrimaryPart, kinEnergyProj, Z, A, 159 fDirectPrimaryPart, kinEnergyProj, Z, A, E2, 1.e20); 160 160 161 dSigmadEprod = (sigma1 - sigma2) / dE; 161 dSigmadEprod = (sigma1 - sigma2) / dE; 162 162 163 if(dSigmadEprod > 1.) 163 if(dSigmadEprod > 1.) 164 { 164 { 165 G4cout << "sigma1 " << kinEnergyProj / M 165 G4cout << "sigma1 " << kinEnergyProj / MeV << '\t' << kinEnergyProd / MeV 166 << '\t' << sigma1 << G4endl; 166 << '\t' << sigma1 << G4endl; 167 G4cout << "sigma2 " << kinEnergyProj / M 167 G4cout << "sigma2 " << kinEnergyProj / MeV << '\t' << kinEnergyProd / MeV 168 << '\t' << sigma2 << G4endl; 168 << '\t' << sigma2 << G4endl; 169 G4cout << "dsigma " << kinEnergyProj / M 169 G4cout << "dsigma " << kinEnergyProj / MeV << '\t' << kinEnergyProd / MeV 170 << '\t' << dSigmadEprod << G4endl 170 << '\t' << dSigmadEprod << G4endl; 171 } 171 } 172 172 173 if(fDirectModel == fBetheBlochDirectEMMode 173 if(fDirectModel == fBetheBlochDirectEMModel) 174 { 174 { 175 // correction of differential cross sect 175 // correction of differential cross section at high energy to correct for 176 // the suppression of particle at second 176 // the suppression of particle at secondary at high energy used in the 177 // Bethe Bloch Model. This correction co 177 // Bethe Bloch Model. This correction consist to multiply by g the 178 // probability function used to test the 178 // probability function used to test the rejection of a secondary Source 179 // code taken from G4BetheBlochModel::Sa 179 // code taken from G4BetheBlochModel::SampleSecondaries 180 180 181 G4double deltaKinEnergy = kinEnergyProd; 181 G4double deltaKinEnergy = kinEnergyProd; 182 182 183 G4double x = fFormFact * deltaKinEnergy; 183 G4double x = fFormFact * deltaKinEnergy; 184 if(x > 1.e-6) 184 if(x > 1.e-6) 185 { 185 { 186 G4double totEnergy = kinEnergyProj + f 186 G4double totEnergy = kinEnergyProj + fMass; 187 G4double etot2 = totEnergy * totEn 187 G4double etot2 = totEnergy * totEnergy; 188 G4double beta2 = kinEnergyProj * (kinE 188 G4double beta2 = kinEnergyProj * (kinEnergyProj + 2.0 * fMass) / etot2; 189 G4double f1 = 0.0; 189 G4double f1 = 0.0; 190 G4double f = 1.0 - beta2 * deltaKi 190 G4double f = 1.0 - beta2 * deltaKinEnergy / Tmax; 191 if(0.5 == fSpin) 191 if(0.5 == fSpin) 192 { 192 { 193 f1 = 0.5 * deltaKinEnergy * deltaKin 193 f1 = 0.5 * deltaKinEnergy * deltaKinEnergy / etot2; 194 f += f1; 194 f += f1; 195 } 195 } 196 G4double x1 = 1.0 + x; 196 G4double x1 = 1.0 + x; 197 G4double gg = 1.0 / (x1 * x1); 197 G4double gg = 1.0 / (x1 * x1); 198 if(0.5 == fSpin) 198 if(0.5 == fSpin) 199 { 199 { 200 G4double x2 = 200 G4double x2 = 201 0.5 * electron_mass_c2 * deltaKinE 201 0.5 * electron_mass_c2 * deltaKinEnergy / (fMass * fMass); 202 gg *= (1.0 + fMagMoment2 * (x2 - f1 202 gg *= (1.0 + fMagMoment2 * (x2 - f1 / f) / (1.0 + x2)); 203 } 203 } 204 if(gg > 1.0) 204 if(gg > 1.0) 205 { 205 { 206 G4cout << "### G4BetheBlochModel in 206 G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: gg= " << gg 207 << G4endl; 207 << G4endl; 208 gg = 1.; 208 gg = 1.; 209 } 209 } 210 dSigmadEprod *= gg; 210 dSigmadEprod *= gg; 211 } 211 } 212 } 212 } 213 } 213 } 214 return dSigmadEprod; 214 return dSigmadEprod; 215 } 215 } 216 216 217 ////////////////////////////////////////////// 217 //////////////////////////////////////////////////////////////////////////////// 218 void G4AdjointIonIonisationModel::SetIon(G4Par 218 void G4AdjointIonIonisationModel::SetIon(G4ParticleDefinition* adj_ion, 219 G4Par 219 G4ParticleDefinition* fwd_ion) 220 { 220 { 221 fDirectPrimaryPart = fwd_ion; 221 fDirectPrimaryPart = fwd_ion; 222 fAdjEquivDirectPrimPart = adj_ion; 222 fAdjEquivDirectPrimPart = adj_ion; 223 223 224 DefineProjectileProperty(); 224 DefineProjectileProperty(); 225 } 225 } 226 226 227 ////////////////////////////////////////////// 227 //////////////////////////////////////////////////////////////////////////////// 228 void G4AdjointIonIonisationModel::CorrectPostS 228 void G4AdjointIonIonisationModel::CorrectPostStepWeight( 229 G4ParticleChange* fParticleChange, G4double 229 G4ParticleChange* fParticleChange, G4double old_weight, 230 G4double adjointPrimKinEnergy, G4double proj 230 G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool) 231 { 231 { 232 // It is needed because the direct cross sec 232 // It is needed because the direct cross section used to compute the 233 // differential cross section is not the one 233 // differential cross section is not the one used in 234 // the direct model where the GenericIon stu 234 // the direct model where the GenericIon stuff is considered with correction 235 // of effective charge. In the direct model 235 // of effective charge. In the direct model the samnepl of secondaries does 236 // not reflect the integral cross section. T 236 // not reflect the integral cross section. The integral fwd cross section that 237 // we used to compute the differential CS ma 237 // we used to compute the differential CS match the sample of secondaries in 238 // the forward case despite the fact that it 238 // the forward case despite the fact that its is not the same total CS than in 239 // the FWD case. For this reason an extra we 239 // the FWD case. For this reason an extra weight correction is needed at the 240 // end. 240 // end. 241 241 242 G4double new_weight = old_weight; 242 G4double new_weight = old_weight; 243 243 244 // the correction of CS due to the problem e 244 // the correction of CS due to the problem explained above 245 G4double kinEnergyProjScaled = fMassRatio * 245 G4double kinEnergyProjScaled = fMassRatio * projectileKinEnergy; 246 fDirectModel = fBraggIonDire 246 fDirectModel = fBraggIonDirectEMModel; 247 if(kinEnergyProjScaled > 2. * MeV && !fUseOn 247 if(kinEnergyProjScaled > 2. * MeV && !fUseOnlyBragg) 248 fDirectModel = fBetheBlochDirectEMModel; 248 fDirectModel = fBetheBlochDirectEMModel; 249 G4double UsedFwdCS = fDirectModel->ComputeCr 249 G4double UsedFwdCS = fDirectModel->ComputeCrossSectionPerAtom( 250 fDirectPrimaryPart, projectileKinEnergy, 1 250 fDirectPrimaryPart, projectileKinEnergy, 1, 1, fTcutSecond, 1.e20); 251 G4double chargeSqRatio = 1.; 251 G4double chargeSqRatio = 1.; 252 if(fChargeSquare > 1.) 252 if(fChargeSquare > 1.) 253 chargeSqRatio = fDirectModel->GetChargeSqu 253 chargeSqRatio = fDirectModel->GetChargeSquareRatio( 254 fDirectPrimaryPart, fCurrentMaterial, pr 254 fDirectPrimaryPart, fCurrentMaterial, projectileKinEnergy); 255 G4double CorrectFwdCS = 255 G4double CorrectFwdCS = 256 chargeSqRatio * fDirectModel->ComputeCross 256 chargeSqRatio * fDirectModel->ComputeCrossSectionPerAtom( 257 G4GenericIon::GenericIon 257 G4GenericIon::GenericIon(), kinEnergyProjScaled, 1, 1, 258 fTcutSecond, 1.e20); 258 fTcutSecond, 1.e20); 259 // May be some check is needed if UsedFwdCS 259 // May be some check is needed if UsedFwdCS ==0 probably that then we should 260 // avoid a secondary to be produced, 260 // avoid a secondary to be produced, 261 if(UsedFwdCS > 0.) 261 if(UsedFwdCS > 0.) 262 new_weight *= CorrectFwdCS / UsedFwdCS; 262 new_weight *= CorrectFwdCS / UsedFwdCS; 263 263 264 // additional CS correction needed for cross 264 // additional CS correction needed for cross section biasing in general. 265 // May be wrong for ions. Most of the time n 265 // May be wrong for ions. Most of the time not used. 266 new_weight *= 266 new_weight *= 267 G4AdjointCSManager::GetAdjointCSManager()- 267 G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection() / 268 fCsBiasingFactor; 268 fCsBiasingFactor; 269 269 270 new_weight *= projectileKinEnergy / adjointP 270 new_weight *= projectileKinEnergy / adjointPrimKinEnergy; 271 271 272 fParticleChange->SetParentWeightByProcess(fa 272 fParticleChange->SetParentWeightByProcess(false); 273 fParticleChange->SetSecondaryWeightByProcess 273 fParticleChange->SetSecondaryWeightByProcess(false); 274 fParticleChange->ProposeParentWeight(new_wei 274 fParticleChange->ProposeParentWeight(new_weight); 275 } 275 } 276 276 277 ////////////////////////////////////////////// 277 //////////////////////////////////////////////////////////////////////////////// 278 void G4AdjointIonIonisationModel::DefineProjec 278 void G4AdjointIonIonisationModel::DefineProjectileProperty() 279 { 279 { 280 // Slightly modified code taken from G4Bethe 280 // Slightly modified code taken from G4BetheBlochModel::SetParticle 281 G4String pname = fDirectPrimaryPart->GetPart 281 G4String pname = fDirectPrimaryPart->GetParticleName(); 282 282 283 fMass = fDirectPrimaryPart->GetPDG 283 fMass = fDirectPrimaryPart->GetPDGMass(); 284 fMassRatio = G4GenericIon::GenericIon() 284 fMassRatio = G4GenericIon::GenericIon()->GetPDGMass() / fMass; 285 fSpin = fDirectPrimaryPart->GetPDG 285 fSpin = fDirectPrimaryPart->GetPDGSpin(); 286 G4double q = fDirectPrimaryPart->GetPDG 286 G4double q = fDirectPrimaryPart->GetPDGCharge() / eplus; 287 fChargeSquare = q * q; 287 fChargeSquare = q * q; 288 fRatio = electron_mass_c2 / fMass; 288 fRatio = electron_mass_c2 / fMass; 289 fOnePlusRatio2 = (1. + fRatio) * (1. + fRat 289 fOnePlusRatio2 = (1. + fRatio) * (1. + fRatio); 290 fOneMinusRatio2 = (1. - fRatio) * (1. - fRat 290 fOneMinusRatio2 = (1. - fRatio) * (1. - fRatio); 291 G4double magmom = fDirectPrimaryPart->GetPDG 291 G4double magmom = fDirectPrimaryPart->GetPDGMagneticMoment() * fMass / 292 (0.5 * eplus * hbar_Planck 292 (0.5 * eplus * hbar_Planck * c_squared); 293 fMagMoment2 = magmom * magmom - 1.0; 293 fMagMoment2 = magmom * magmom - 1.0; 294 if(fDirectPrimaryPart->GetLeptonNumber() == 294 if(fDirectPrimaryPart->GetLeptonNumber() == 0) 295 { 295 { 296 G4double x = 0.8426 * GeV; 296 G4double x = 0.8426 * GeV; 297 if(fSpin == 0.0 && fMass < GeV) 297 if(fSpin == 0.0 && fMass < GeV) 298 { 298 { 299 x = 0.736 * GeV; 299 x = 0.736 * GeV; 300 } 300 } 301 else if(fMass > GeV) 301 else if(fMass > GeV) 302 { 302 { 303 x /= G4NistManager::Instance()->GetZ13(f 303 x /= G4NistManager::Instance()->GetZ13(fMass / proton_mass_c2); 304 } 304 } 305 fFormFact = 2.0 * electron_mass_c2 / (x * 305 fFormFact = 2.0 * electron_mass_c2 / (x * x); 306 } 306 } 307 } 307 } 308 308 309 ////////////////////////////////////////////// 309 ////////////////////////////////////////////////////////////////////////////// 310 G4double G4AdjointIonIonisationModel::GetSecon 310 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMaxForScatProjToProj( 311 G4double primAdjEnergy) 311 G4double primAdjEnergy) 312 { 312 { 313 return primAdjEnergy * fOnePlusRatio2 / 313 return primAdjEnergy * fOnePlusRatio2 / 314 (fOneMinusRatio2 - 2. * fRatio * prim 314 (fOneMinusRatio2 - 2. * fRatio * primAdjEnergy / fMass); 315 } 315 } 316 316 317 ////////////////////////////////////////////// 317 ////////////////////////////////////////////////////////////////////////////// 318 G4double G4AdjointIonIonisationModel::GetSecon 318 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMinForScatProjToProj( 319 G4double primAdjEnergy, G4double tcut) 319 G4double primAdjEnergy, G4double tcut) 320 { 320 { 321 return primAdjEnergy + tcut; 321 return primAdjEnergy + tcut; 322 } 322 } 323 323 324 ////////////////////////////////////////////// 324 ////////////////////////////////////////////////////////////////////////////// 325 G4double G4AdjointIonIonisationModel::GetSecon 325 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMaxForProdToProj( 326 G4double) 326 G4double) 327 { 327 { 328 return GetHighEnergyLimit(); 328 return GetHighEnergyLimit(); 329 } 329 } 330 330 331 ////////////////////////////////////////////// 331 ////////////////////////////////////////////////////////////////////////////// 332 G4double G4AdjointIonIonisationModel::GetSecon 332 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMinForProdToProj( 333 G4double primAdjEnergy) 333 G4double primAdjEnergy) 334 { 334 { 335 return (2. * primAdjEnergy - 4. * fMass + 335 return (2. * primAdjEnergy - 4. * fMass + 336 std::sqrt(4. * primAdjEnergy * primA 336 std::sqrt(4. * primAdjEnergy * primAdjEnergy + 16. * fMass * fMass + 337 8. * primAdjEnergy * fMass 337 8. * primAdjEnergy * fMass * (1. / fRatio + fRatio))) / 338 4.; 338 4.; 339 } 339 } 340 340