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 "G4AdjointPhotoElectricModel.hh" 28 29 #include "G4AdjointCSManager.hh" 30 #include "G4AdjointElectron.hh" 31 #include "G4AdjointGamma.hh" 32 #include "G4Gamma.hh" 33 #include "G4ParticleChange.hh" 34 #include "G4PEEffectFluoModel.hh" 35 #include "G4PhysicalConstants.hh" 36 #include "G4TrackStatus.hh" 37 38 ////////////////////////////////////////////// 39 G4AdjointPhotoElectricModel::G4AdjointPhotoEle 40 : G4VEmAdjointModel("AdjointPEEffect") 41 42 { 43 SetUseMatrix(false); 44 SetApplyCutInRange(false); 45 46 fAdjEquivDirectPrimPart = G4AdjointGamma:: 47 fAdjEquivDirectSecondPart = G4AdjointElectro 48 fDirectPrimaryPart = G4Gamma::Gamma() 49 fSecondPartSameType = false; 50 fDirectModel = new G4PEEffectFl 51 } 52 53 ////////////////////////////////////////////// 54 G4AdjointPhotoElectricModel::~G4AdjointPhotoEl 55 56 ////////////////////////////////////////////// 57 void G4AdjointPhotoElectricModel::SampleSecond 58 const G4Track& aTrack, G4bool isScatProjToPr 59 G4ParticleChange* fParticleChange) 60 { 61 if(isScatProjToProj) 62 return; 63 64 // Compute the fTotAdjointCS vectors if not 65 // couple and electron energy 66 const G4DynamicParticle* aDynPart = aTrack.G 67 G4double electronEnergy = aDynPart 68 G4ThreeVector electronDirection = aDynPart 69 fPreStepAdjointCS = 70 fTotAdjointCS; // The last computed CS wa 71 AdjointCrossSection(aTrack.GetMaterialCutsCo 72 isScatProjToProj); 73 fPostStepAdjointCS = fTotAdjointCS; 74 75 // Sample element 76 const G4ElementVector* theElementVector = 77 fCurrentMaterial->GetElementVector(); 78 size_t nelm = fCurrentMaterial->GetNumb 79 G4double rand_CS = G4UniformRand() * fXsec[n 80 for(fIndexElement = 0; fIndexElement < nelm 81 { 82 if(rand_CS < fXsec[fIndexElement]) 83 break; 84 } 85 86 // Sample shell and binding energy 87 G4int nShells = (*theElementVector)[fIndexEl 88 rand_CS = fShellProb[fIndexElement][nS 89 G4int i; 90 for(i = 0; i < nShells - 1; ++i) 91 { 92 if(rand_CS < fShellProb[fIndexElement][i]) 93 break; 94 } 95 G4double gammaEnergy = 96 electronEnergy + (*theElementVector)[fInde 97 98 // Sample cos theta 99 // Copy of the G4PEEfectFluoModel cos theta 100 // ElecCosThetaDistribution. This method can 101 // G4PEEffectFluoModel because it is a frien 102 G4double cos_theta = 1.; 103 G4double gamma = 1. + electronEnergy / e 104 if(gamma <= 5.) 105 { 106 G4double beta = std::sqrt(gamma * gamma - 107 G4double b = 0.5 * gamma * (gamma - 1.) 108 109 G4double rndm, term, greject, grejsup; 110 if(gamma < 2.) 111 grejsup = gamma * gamma * (1. + b - beta 112 else 113 grejsup = gamma * gamma * (1. + b + beta 114 115 do 116 { 117 rndm = 1. - 2. * G4UniformRand(); 118 cos_theta = (rndm + beta) / (rndm * beta 119 term = 1. - beta * cos_theta; 120 greject = (1. - cos_theta * cos_theta) * 121 // Loop checking, 07-Aug-2015, Vladimir 122 } while(greject < G4UniformRand() * grejsu 123 } 124 125 // direction of the adjoint gamma electron 126 G4double sin_theta = std::sqrt(1. - cos_thet 127 G4double phi = twopi * G4UniformRand() 128 G4double dirx = sin_theta * std::cos(ph 129 G4double diry = sin_theta * std::sin(ph 130 G4double dirz = cos_theta; 131 G4ThreeVector adjoint_gammaDirection(dirx, d 132 adjoint_gammaDirection.rotateUz(electronDire 133 134 // Weight correction 135 CorrectPostStepWeight(fParticleChange, aTrac 136 gammaEnergy, isScatPro 137 138 // Create secondary and modify fParticleChan 139 G4DynamicParticle* anAdjointGamma = new G4Dy 140 G4AdjointGamma::AdjointGamma(), adjoint_ga 141 142 fParticleChange->ProposeTrackStatus(fStopAnd 143 fParticleChange->AddSecondary(anAdjointGamma 144 } 145 146 ////////////////////////////////////////////// 147 void G4AdjointPhotoElectricModel::CorrectPostS 148 G4ParticleChange* fParticleChange, G4double 149 G4double adjointPrimKinEnergy, G4double proj 150 { 151 G4double new_weight = old_weight; 152 153 G4double w_corr = 154 G4AdjointCSManager::GetAdjointCSManager()- 155 fFactorCSBiasing; 156 w_corr *= fPostStepAdjointCS / fPreStepAdjoi 157 158 new_weight *= w_corr * projectileKinEnergy / 159 fParticleChange->SetParentWeightByProcess(fa 160 fParticleChange->SetSecondaryWeightByProcess 161 fParticleChange->ProposeParentWeight(new_wei 162 } 163 164 ////////////////////////////////////////////// 165 G4double G4AdjointPhotoElectricModel::AdjointC 166 const G4MaterialCutsCouple* aCouple, G4doubl 167 G4bool isScatProjToProj) 168 { 169 if(isScatProjToProj) 170 return 0.; 171 172 G4double totBiasedAdjointCS = 0.; 173 if(aCouple != fCurrentCouple || fCurrenteEne 174 { 175 fTotAdjointCS = 0.; 176 DefineCurrentMaterialAndElectronEnergy(aCo 177 const G4ElementVector* theElementVector = 178 fCurrentMaterial->GetElementVector(); 179 const G4double* theAtomNumDensityVector = 180 fCurrentMaterial->GetVecNbOfAtomsPerVolu 181 size_t nelm = fCurrentMaterial->GetNumberO 182 for(fIndexElement = 0; fIndexElement < nel 183 { 184 fTotAdjointCS += AdjointCrossSectionPerA 185 (*theElementVector)[f 186 theAtomNumDensityVector 187 fXsec[fIndexElement] = fTotAdjointCS; 188 } 189 190 totBiasedAdjointCS = std::min(fTotAdjointC 191 fFactorCSBiasing = totBiasedAdjointCS / 192 } 193 return totBiasedAdjointCS; 194 } 195 196 ////////////////////////////////////////////// 197 G4double G4AdjointPhotoElectricModel::AdjointC 198 const G4Element* anElement, G4double electro 199 { 200 G4int nShells = anElement->GetNbOfAto 201 G4double Z = anElement->GetZ(); 202 G4double gammaEnergy = electronEnergy + anEl 203 G4double CS = fDirectModel->Compute 204 G4Gamma::Gamma(), gammaEnergy, Z, 0., 0., 205 G4double adjointCS = 0.; 206 if(CS > 0.) 207 adjointCS += CS / gammaEnergy; 208 fShellProb[fIndexElement][0] = adjointCS; 209 for(G4int i = 1; i < nShells; ++i) 210 { 211 G4double Bi1 = anElement->GetAtomicShell(i 212 G4double Bi = anElement->GetAtomicShell(i 213 if(electronEnergy < Bi1 - Bi) 214 { 215 gammaEnergy = electronEnergy + Bi; 216 CS = fDirectModel->ComputeCross 217 218 if(CS > 0.) 219 adjointCS += CS / gammaEnergy; 220 } 221 fShellProb[fIndexElement][i] = adjointCS; 222 } 223 adjointCS *= electronEnergy; 224 return adjointCS; 225 } 226 227 ////////////////////////////////////////////// 228 void G4AdjointPhotoElectricModel::DefineCurren 229 const G4MaterialCutsCouple* couple, G4double 230 { 231 fCurrentCouple = const_cast<G4MaterialCuts 232 fCurrentMaterial = const_cast<G4Material*>(c 233 fCurrenteEnergy = anEnergy; 234 fDirectModel->SetCurrentCouple(couple); 235 } 236