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 #include "G4AdjointBremsstrahlungModel.hh" 28 29 #include "G4AdjointCSManager.hh" 30 #include "G4AdjointElectron.hh" 31 #include "G4AdjointGamma.hh" 32 #include "G4Electron.hh" 33 #include "G4EmModelManager.hh" 34 #include "G4Gamma.hh" 35 #include "G4ParticleChange.hh" 36 #include "G4PhysicalConstants.hh" 37 #include "G4SeltzerBergerModel.hh" 38 #include "G4SystemOfUnits.hh" 39 #include "G4TrackStatus.hh" 40 41 ////////////////////////////////////////////// 42 G4AdjointBremsstrahlungModel::G4AdjointBremsst 43 : G4VEmAdjointModel("AdjointeBremModel") 44 { 45 fDirectModel = aModel; 46 Initialize(); 47 } 48 49 ////////////////////////////////////////////// 50 G4AdjointBremsstrahlungModel::G4AdjointBremsst 51 : G4VEmAdjointModel("AdjointeBremModel") 52 { 53 fDirectModel = new G4SeltzerBergerModel(); 54 Initialize(); 55 } 56 57 ////////////////////////////////////////////// 58 void G4AdjointBremsstrahlungModel::Initialize( 59 { 60 SetUseMatrix(false); 61 SetUseMatrixPerElement(false); 62 63 fEmModelManagerForFwdModels = new G4EmModelM 64 fEmModelManagerForFwdModels->AddEmModel(1, f 65 SetApplyCutInRange(true); 66 67 fElectron = G4Electron::Electron(); 68 fGamma = G4Gamma::Gamma(); 69 70 fAdjEquivDirectPrimPart = G4AdjointElectro 71 fAdjEquivDirectSecondPart = G4AdjointGamma:: 72 fDirectPrimaryPart = fElectron; 73 fSecondPartSameType = false; 74 75 fCSManager = G4AdjointCSManager::GetAdjointC 76 } 77 78 ////////////////////////////////////////////// 79 G4AdjointBremsstrahlungModel::~G4AdjointBremss 80 { 81 if(fEmModelManagerForFwdModels) 82 delete fEmModelManagerForFwdModels; 83 } 84 85 ////////////////////////////////////////////// 86 void G4AdjointBremsstrahlungModel::SampleSecon 87 const G4Track& aTrack, G4bool isScatProjToPr 88 G4ParticleChange* fParticleChange) 89 { 90 if(!fUseMatrix) 91 return RapidSampleSecondaries(aTrack, isSc 92 93 const G4DynamicParticle* theAdjointPrimary = 94 DefineCurrentMaterial(aTrack.GetMaterialCuts 95 96 G4double adjointPrimKinEnergy = theAdjoint 97 G4double adjointPrimTotalEnergy = theAdjoint 98 99 if(adjointPrimKinEnergy > GetHighEnergyLimit 100 { 101 return; 102 } 103 104 G4double projectileKinEnergy = 105 SampleAdjSecEnergyFromCSMatrix(adjointPrim 106 107 // Weight correction 108 CorrectPostStepWeight(fParticleChange, aTrac 109 adjointPrimKinEnergy, 110 isScatProjToProj); 111 112 // Kinematic 113 G4double projectileM0 = fAdjEquivDi 114 G4double projectileTotalEnergy = projectileM 115 G4double projectileP2 = 116 projectileTotalEnergy * projectileTotalEne 117 G4double projectileP = std::sqrt(projectileP 118 119 // Angle of the gamma direction with the pro 120 // G4eBremsstrahlungModel 121 G4double u; 122 if(0.25 > G4UniformRand()) 123 u = -std::log(G4UniformRand() * G4UniformR 124 else 125 u = -std::log(G4UniformRand() * G4UniformR 126 127 G4double theta = u * electron_mass_c2 / proj 128 G4double sint = std::sin(theta); 129 G4double cost = std::cos(theta); 130 131 G4double phi = twopi * G4UniformRand(); 132 133 G4ThreeVector projectileMomentum = 134 G4ThreeVector(std::cos(phi) * sint, std::s 135 projectileP; // gamma frame 136 if(isScatProjToProj) 137 { // the adjoint primary is the scattered e 138 G4ThreeVector gammaMomentum = 139 (projectileTotalEnergy - adjointPrimTota 140 G4ThreeVector(0., 0., 1.); 141 G4ThreeVector dirProd = projectileMomentum 142 G4double cost1 = std::cos(dirProd.a 143 G4double sint1 = std::sqrt(1. - cos 144 projectileMomentum = 145 G4ThreeVector(std::cos(phi) * sint1, std 146 projectileP; 147 } 148 149 projectileMomentum.rotateUz(theAdjointPrimar 150 151 if(!isScatProjToProj) 152 { // kill the primary and add a secondary 153 fParticleChange->ProposeTrackStatus(fStopA 154 fParticleChange->AddSecondary( 155 new G4DynamicParticle(fAdjEquivDirectPri 156 } 157 else 158 { 159 fParticleChange->ProposeEnergy(projectileK 160 fParticleChange->ProposeMomentumDirection( 161 } 162 } 163 164 ////////////////////////////////////////////// 165 void G4AdjointBremsstrahlungModel::RapidSample 166 const G4Track& aTrack, G4bool isScatProjToPr 167 G4ParticleChange* fParticleChange) 168 { 169 const G4DynamicParticle* theAdjointPrimary = 170 DefineCurrentMaterial(aTrack.GetMaterialCuts 171 172 G4double adjointPrimKinEnergy = theAdjoint 173 G4double adjointPrimTotalEnergy = theAdjoint 174 175 if(adjointPrimKinEnergy > GetHighEnergyLimit 176 { 177 return; 178 } 179 180 G4double projectileKinEnergy = 0.; 181 G4double gammaEnergy = 0.; 182 G4double diffCSUsed = 0.; 183 if(!isScatProjToProj) 184 { 185 gammaEnergy = adjointPrimKinEnergy; 186 G4double Emax = GetSecondAdjEnergyMaxForPr 187 G4double Emin = GetSecondAdjEnergyMinForPr 188 if(Emin >= Emax) 189 return; 190 projectileKinEnergy = Emin * std::pow(Emax 191 diffCSUsed = fCsBiasingFactor * f 192 } 193 else 194 { 195 G4double Emax = 196 GetSecondAdjEnergyMaxForScatProjToProj(a 197 G4double Emin = 198 GetSecondAdjEnergyMinForScatProjToProj(a 199 if(Emin >= Emax) 200 return; 201 G4double f1 = (Emin - adjointPrimKinEnergy 202 G4double f2 = (Emax - adjointPrimKinEnergy 203 projectileKinEnergy = 204 adjointPrimKinEnergy / (1. - f1 * std::p 205 gammaEnergy = projectileKinEnergy - adjoin 206 diffCSUsed = 207 fLastCZ * adjointPrimKinEnergy / project 208 } 209 210 // Weight correction: 211 // First w_corr is set to the ratio between 212 // if this has to be done in the model. 213 // For the case of forced interaction this w 214 // the forced interaction. It is important 215 // creation of the secondary 216 G4double w_corr = fOutsideWeightFactor; 217 if(fInModelWeightCorr) 218 { 219 w_corr = fCSManager->GetPostStepWeightCorr 220 } 221 222 // Then another correction is needed due to 223 // differential CS has been used rather than 224 // direct model Here we consider the true di 225 // numerical differentiation over Tcut of th 226 // Migdal term. Basically any other differen 227 // (example Penelope). 228 G4double diffCS = DiffCrossSectionPerVolumeP 229 fCurrentMaterial, projectileKinEnergy, gam 230 w_corr *= diffCS / diffCSUsed; 231 232 G4double new_weight = aTrack.GetWeight() * w 233 fParticleChange->SetParentWeightByProcess(fa 234 fParticleChange->SetSecondaryWeightByProcess 235 fParticleChange->ProposeParentWeight(new_wei 236 237 // Kinematic 238 G4double projectileM0 = fAdjEquivDi 239 G4double projectileTotalEnergy = projectileM 240 G4double projectileP2 = 241 projectileTotalEnergy * projectileTotalEne 242 G4double projectileP = std::sqrt(projectileP 243 244 // Use the angular model of the forward mode 245 // Dummy dynamic particle to use the model 246 G4DynamicParticle* aDynPart = 247 new G4DynamicParticle(fElectron, G4ThreeVe 248 249 // Get the element from the direct model 250 const G4Element* elm = fDirectModel->SelectR 251 fCurrentCouple, fElectron, projectileKinEn 252 G4int Z = elm->GetZasInt(); 253 G4double energy = aDynPart->GetTotalEnergy() 254 G4ThreeVector projectileMomentum = 255 fDirectModel->GetAngularDistribution()->Sa 256 fC 257 G4double phi = projectileMomentum.getPhi(); 258 259 if(isScatProjToProj) 260 { // the adjoint primary is the scattered e 261 G4ThreeVector gammaMomentum = 262 (projectileTotalEnergy - adjointPrimTota 263 G4ThreeVector(0., 0., 1.); 264 G4ThreeVector dirProd = projectileMomentum 265 G4double cost1 = std::cos(dirProd.a 266 G4double sint1 = std::sqrt(1. - cos 267 projectileMomentum = 268 G4ThreeVector(std::cos(phi) * sint1, std 269 projectileP; 270 } 271 272 projectileMomentum.rotateUz(theAdjointPrimar 273 274 if(!isScatProjToProj) 275 { // kill the primary and add a secondary 276 fParticleChange->ProposeTrackStatus(fStopA 277 fParticleChange->AddSecondary( 278 new G4DynamicParticle(fAdjEquivDirectPri 279 } 280 else 281 { 282 fParticleChange->ProposeEnergy(projectileK 283 fParticleChange->ProposeMomentumDirection( 284 } 285 } 286 287 ////////////////////////////////////////////// 288 G4double G4AdjointBremsstrahlungModel::DiffCro 289 const G4Material* aMaterial, 290 G4double kinEnergyProj, // kin energy of pr 291 G4double kinEnergyProd // kinetic energy o 292 ) 293 { 294 if(!fIsDirectModelInitialised) 295 { 296 fEmModelManagerForFwdModels->Initialise(fE 297 fIsDirectModelInitialised = true; 298 } 299 return G4VEmAdjointModel::DiffCrossSectionPe 300 aMaterial, kinEnergyProj, kinEnergyProd); 301 } 302 303 ////////////////////////////////////////////// 304 G4double G4AdjointBremsstrahlungModel::Adjoint 305 const G4MaterialCutsCouple* aCouple, G4doubl 306 G4bool isScatProjToProj) 307 { 308 static constexpr G4double maxEnergy = 100. * 309 // 2.78.. == std::exp(1.) 310 if(!fIsDirectModelInitialised) 311 { 312 fEmModelManagerForFwdModels->Initialise(fE 313 fIsDirectModelInitialised = true; 314 } 315 if(fUseMatrix) 316 return G4VEmAdjointModel::AdjointCrossSect 317 318 DefineCurrentMaterial(aCouple); 319 G4double Cross = 0.; 320 // this gives the constant above 321 fLastCZ = fDirectModel->CrossSectionPerVolum 322 aCouple->GetMaterial(), fDirectPrimaryPart 323 324 if(!isScatProjToProj) 325 { 326 G4double Emax_proj = GetSecondAdjEnergyMax 327 G4double Emin_proj = GetSecondAdjEnergyMin 328 if(Emax_proj > Emin_proj && primEnergy > f 329 Cross = fCsBiasingFactor * fLastCZ * std 330 } 331 else 332 { 333 G4double Emax_proj = GetSecondAdjEnergyMax 334 G4double Emin_proj = 335 GetSecondAdjEnergyMinForScatProjToProj(p 336 if(Emax_proj > Emin_proj) 337 Cross = fLastCZ * std::log((Emax_proj - 338 Emax_proj / ( 339 } 340 return Cross; 341 } 342