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 // $Id: G4AdjointPhotoElectricModel.cc,v 1.6 2010/11/11 11:51:56 ldesorgh Exp $ >> 27 // GEANT4 tag $Name: geant4-09-04 $ >> 28 // 27 #include "G4AdjointPhotoElectricModel.hh" 29 #include "G4AdjointPhotoElectricModel.hh" 28 << 29 #include "G4AdjointCSManager.hh" 30 #include "G4AdjointCSManager.hh" >> 31 >> 32 >> 33 #include "G4Integrator.hh" >> 34 #include "G4TrackStatus.hh" >> 35 #include "G4ParticleChange.hh" 30 #include "G4AdjointElectron.hh" 36 #include "G4AdjointElectron.hh" >> 37 #include "G4Gamma.hh" 31 #include "G4AdjointGamma.hh" 38 #include "G4AdjointGamma.hh" 32 #include "G4Gamma.hh" << 39 33 #include "G4ParticleChange.hh" << 34 #include "G4PEEffectFluoModel.hh" << 35 #include "G4PhysicalConstants.hh" << 36 #include "G4TrackStatus.hh" << 37 40 38 ////////////////////////////////////////////// 41 //////////////////////////////////////////////////////////////////////////////// 39 G4AdjointPhotoElectricModel::G4AdjointPhotoEle << 42 // 40 : G4VEmAdjointModel("AdjointPEEffect") << 43 G4AdjointPhotoElectricModel::G4AdjointPhotoElectricModel(): >> 44 G4VEmAdjointModel("AdjointPEEffect") 41 45 42 { << 46 { SetUseMatrix(false); 43 SetUseMatrix(false); << 44 SetApplyCutInRange(false); 47 SetApplyCutInRange(false); 45 << 48 current_eEnergy =0.; 46 fAdjEquivDirectPrimPart = G4AdjointGamma:: << 49 totAdjointCS=0.; 47 fAdjEquivDirectSecondPart = G4AdjointElectro << 50 theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma(); 48 fDirectPrimaryPart = G4Gamma::Gamma() << 51 theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron(); 49 fSecondPartSameType = false; << 52 theDirectPrimaryPartDef=G4Gamma::Gamma(); 50 fDirectModel = new G4PEEffectFl << 53 second_part_of_same_type=false; >> 54 theDirectPEEffectModel = new G4PEEffectModel(); >> 55 51 } 56 } 52 << 53 ////////////////////////////////////////////// 57 //////////////////////////////////////////////////////////////////////////////// 54 G4AdjointPhotoElectricModel::~G4AdjointPhotoEl << 58 // 55 << 59 G4AdjointPhotoElectricModel::~G4AdjointPhotoElectricModel() 56 ////////////////////////////////////////////// << 60 {;} 57 void G4AdjointPhotoElectricModel::SampleSecond << 61 58 const G4Track& aTrack, G4bool isScatProjToPr << 62 //////////////////////////////////////////////////////////////////////////////// 59 G4ParticleChange* fParticleChange) << 63 // 60 { << 64 void G4AdjointPhotoElectricModel::SampleSecondaries(const G4Track& aTrack, 61 if(isScatProjToProj) << 65 G4bool IsScatProjToProjCase, 62 return; << 66 G4ParticleChange* fParticleChange) 63 << 67 { if (IsScatProjToProjCase) return ; 64 // Compute the fTotAdjointCS vectors if not << 68 65 // couple and electron energy << 69 //Compute the totAdjointCS vectors if not already done for the current couple and electron energy 66 const G4DynamicParticle* aDynPart = aTrack.G << 70 //----------------------------------------------------------------------------------------------- 67 G4double electronEnergy = aDynPart << 71 const G4MaterialCutsCouple* aCouple = aTrack.GetMaterialCutsCouple(); 68 G4ThreeVector electronDirection = aDynPart << 72 const G4DynamicParticle* aDynPart = aTrack.GetDynamicParticle() ; 69 fPreStepAdjointCS = << 73 G4double electronEnergy = aDynPart->GetKineticEnergy(); 70 fTotAdjointCS; // The last computed CS wa << 74 G4ThreeVector electronDirection= aDynPart->GetMomentumDirection() ; 71 AdjointCrossSection(aTrack.GetMaterialCutsCo << 75 pre_step_AdjointCS = totAdjointCS; //The last computed CS was at pre step point 72 isScatProjToProj); << 76 G4double adjCS; 73 fPostStepAdjointCS = fTotAdjointCS; << 77 adjCS = AdjointCrossSection(aCouple, electronEnergy,IsScatProjToProjCase); 74 << 78 post_step_AdjointCS = totAdjointCS; 75 // Sample element << 79 76 const G4ElementVector* theElementVector = << 80 77 fCurrentMaterial->GetElementVector(); << 81 78 size_t nelm = fCurrentMaterial->GetNumb << 82 79 G4double rand_CS = G4UniformRand() * fXsec[n << 83 //Sample element 80 for(fIndexElement = 0; fIndexElement < nelm << 84 //------------- 81 { << 85 const G4ElementVector* theElementVector = currentMaterial->GetElementVector(); 82 if(rand_CS < fXsec[fIndexElement]) << 86 size_t nelm = currentMaterial->GetNumberOfElements(); 83 break; << 87 G4double rand_CS= G4UniformRand()*xsec[nelm-1]; 84 } << 88 for (index_element=0; index_element<nelm-1; index_element++){ 85 << 89 if (rand_CS<xsec[index_element]) break; 86 // Sample shell and binding energy << 90 } 87 G4int nShells = (*theElementVector)[fIndexEl << 91 88 rand_CS = fShellProb[fIndexElement][nS << 92 //Sample shell and binding energy 89 G4int i; << 93 //------------- 90 for(i = 0; i < nShells - 1; ++i) << 94 G4int nShells = (*theElementVector)[index_element]->GetNbOfAtomicShells(); 91 { << 95 rand_CS= shell_prob[index_element][nShells-1]*G4UniformRand(); 92 if(rand_CS < fShellProb[fIndexElement][i]) << 96 G4int i = 0; 93 break; << 97 for (i=0; i<nShells-1; i++){ 94 } << 98 if (rand_CS<shell_prob[index_element][i]) break; 95 G4double gammaEnergy = << 99 } 96 electronEnergy + (*theElementVector)[fInde << 100 G4double gammaEnergy= electronEnergy+(*theElementVector)[index_element]->GetAtomicShell(i); 97 << 101 98 // Sample cos theta << 102 //Sample cos theta 99 // Copy of the G4PEEfectFluoModel cos theta << 103 //Copy of the G4PEEffectModel cos theta sampling method ElecCosThetaDistribution. 100 // ElecCosThetaDistribution. This method can << 104 //This method cannot be used directly from G4PEEffectModel because it is a friend method. I should ask Vladimir to change that 101 // G4PEEffectFluoModel because it is a frien << 105 //------------------------------------------------------------------------------------------------ 102 G4double cos_theta = 1.; << 106 //G4double cos_theta = theDirectPEEffectModel->ElecCosThetaDistribution(electronEnergy); 103 G4double gamma = 1. + electronEnergy / e << 107 104 if(gamma <= 5.) << 108 G4double cos_theta = 1.; 105 { << 109 G4double gamma = 1. + electronEnergy/electron_mass_c2; 106 G4double beta = std::sqrt(gamma * gamma - << 110 if (gamma <= 5.) { 107 G4double b = 0.5 * gamma * (gamma - 1.) << 111 G4double beta = std::sqrt(gamma*gamma-1.)/gamma; 108 << 112 G4double b = 0.5*gamma*(gamma-1.)*(gamma-2); 109 G4double rndm, term, greject, grejsup; << 113 110 if(gamma < 2.) << 114 G4double rndm,term,greject,grejsup; 111 grejsup = gamma * gamma * (1. + b - beta << 115 if (gamma < 2.) grejsup = gamma*gamma*(1.+b-beta*b); 112 else << 116 else grejsup = gamma*gamma*(1.+b+beta*b); 113 grejsup = gamma * gamma * (1. + b + beta << 117 114 << 118 do { rndm = 1.-2*G4UniformRand(); 115 do << 119 cos_theta = (rndm+beta)/(rndm*beta+1.); 116 { << 120 term = 1.-beta*cos_theta; 117 rndm = 1. - 2. * G4UniformRand(); << 121 greject = (1.-cos_theta*cos_theta)*(1.+b*term)/(term*term); 118 cos_theta = (rndm + beta) / (rndm * beta << 122 } while(greject < G4UniformRand()*grejsup); 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 } 123 } 124 << 124 125 // direction of the adjoint gamma electron 125 // direction of the adjoint gamma electron 126 G4double sin_theta = std::sqrt(1. - cos_thet << 126 //--------------------------------------- 127 G4double phi = twopi * G4UniformRand() << 127 128 G4double dirx = sin_theta * std::cos(ph << 128 129 G4double diry = sin_theta * std::sin(ph << 129 G4double sin_theta = std::sqrt(1.-cos_theta*cos_theta); 130 G4double dirz = cos_theta; << 130 G4double Phi = twopi * G4UniformRand(); 131 G4ThreeVector adjoint_gammaDirection(dirx, d << 131 G4double dirx = sin_theta*std::cos(Phi),diry = sin_theta*std::sin(Phi),dirz = cos_theta; >> 132 G4ThreeVector adjoint_gammaDirection(dirx,diry,dirz); 132 adjoint_gammaDirection.rotateUz(electronDire 133 adjoint_gammaDirection.rotateUz(electronDirection); 133 << 134 134 // Weight correction << 135 135 CorrectPostStepWeight(fParticleChange, aTrac << 136 136 gammaEnergy, isScatPro << 137 //Weight correction 137 << 138 //----------------------- 138 // Create secondary and modify fParticleChan << 139 CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(), electronEnergy,gammaEnergy,IsScatProjToProjCase); 139 G4DynamicParticle* anAdjointGamma = new G4Dy << 140 140 G4AdjointGamma::AdjointGamma(), adjoint_ga << 141 141 << 142 >> 143 //Create secondary and modify fParticleChange >> 144 //-------------------------------------------- >> 145 G4DynamicParticle* anAdjointGamma = new G4DynamicParticle ( >> 146 G4AdjointGamma::AdjointGamma(),adjoint_gammaDirection, gammaEnergy); >> 147 >> 148 >> 149 >> 150 >> 151 142 fParticleChange->ProposeTrackStatus(fStopAnd 152 fParticleChange->ProposeTrackStatus(fStopAndKill); 143 fParticleChange->AddSecondary(anAdjointGamma << 153 fParticleChange->AddSecondary(anAdjointGamma); 144 } << 154 145 155 146 ////////////////////////////////////////////// << 156 147 void G4AdjointPhotoElectricModel::CorrectPostS << 148 G4ParticleChange* fParticleChange, G4double << 149 G4double adjointPrimKinEnergy, G4double proj << 150 { << 151 G4double new_weight = old_weight; << 152 157 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 } 158 } 163 159 164 ////////////////////////////////////////////// 160 //////////////////////////////////////////////////////////////////////////////// 165 G4double G4AdjointPhotoElectricModel::AdjointC << 161 // 166 const G4MaterialCutsCouple* aCouple, G4doubl << 162 void G4AdjointPhotoElectricModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange, 167 G4bool isScatProjToProj) << 163 G4double old_weight, 168 { << 164 G4double adjointPrimKinEnergy, 169 if(isScatProjToProj) << 165 G4double projectileKinEnergy , 170 return 0.; << 166 G4bool ) 171 << 167 { 172 G4double totBiasedAdjointCS = 0.; << 168 G4double new_weight=old_weight; 173 if(aCouple != fCurrentCouple || fCurrenteEne << 169 174 { << 170 G4double w_corr =G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection()/factorCSBiasing; 175 fTotAdjointCS = 0.; << 171 w_corr*=post_step_AdjointCS/pre_step_AdjointCS; 176 DefineCurrentMaterialAndElectronEnergy(aCo << 172 new_weight*=w_corr; 177 const G4ElementVector* theElementVector = << 173 new_weight*=projectileKinEnergy/adjointPrimKinEnergy; 178 fCurrentMaterial->GetElementVector(); << 174 fParticleChange->SetParentWeightByProcess(false); 179 const G4double* theAtomNumDensityVector = << 175 fParticleChange->SetSecondaryWeightByProcess(false); 180 fCurrentMaterial->GetVecNbOfAtomsPerVolu << 176 fParticleChange->ProposeParentWeight(new_weight); 181 size_t nelm = fCurrentMaterial->GetNumberO << 177 } 182 for(fIndexElement = 0; fIndexElement < nel << 178 183 { << 179 //////////////////////////////////////////////////////////////////////////////// 184 fTotAdjointCS += AdjointCrossSectionPerA << 180 // 185 (*theElementVector)[f << 181 186 theAtomNumDensityVector << 182 G4double G4AdjointPhotoElectricModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple, 187 fXsec[fIndexElement] = fTotAdjointCS; << 183 G4double electronEnergy, 188 } << 184 G4bool IsScatProjToProjCase) 189 << 185 { 190 totBiasedAdjointCS = std::min(fTotAdjointC << 186 191 fFactorCSBiasing = totBiasedAdjointCS / << 187 >> 188 if (IsScatProjToProjCase) return 0.; >> 189 >> 190 >> 191 if (aCouple !=currentCouple || current_eEnergy !=electronEnergy) { >> 192 totAdjointCS = 0.; >> 193 DefineCurrentMaterialAndElectronEnergy(aCouple, electronEnergy); >> 194 const G4ElementVector* theElementVector = currentMaterial->GetElementVector(); >> 195 const double* theAtomNumDensityVector = currentMaterial->GetVecNbOfAtomsPerVolume(); >> 196 size_t nelm = currentMaterial->GetNumberOfElements(); >> 197 for (index_element=0;index_element<nelm;index_element++){ >> 198 >> 199 totAdjointCS +=AdjointCrossSectionPerAtom((*theElementVector)[index_element],electronEnergy)*theAtomNumDensityVector[index_element]; >> 200 xsec[index_element] = totAdjointCS; >> 201 } >> 202 >> 203 totBiasedAdjointCS=std::min(totAdjointCS,0.01); >> 204 // totBiasedAdjointCS=totAdjointCS; >> 205 factorCSBiasing = totBiasedAdjointCS/totAdjointCS; >> 206 lastCS=totBiasedAdjointCS; >> 207 >> 208 192 } 209 } 193 return totBiasedAdjointCS; 210 return totBiasedAdjointCS; >> 211 >> 212 194 } 213 } >> 214 //////////////////////////////////////////////////////////////////////////////// >> 215 // 195 216 >> 217 G4double G4AdjointPhotoElectricModel::GetAdjointCrossSection(const G4MaterialCutsCouple* aCouple, >> 218 G4double electronEnergy, >> 219 G4bool IsScatProjToProjCase) >> 220 { return AdjointCrossSection(aCouple,electronEnergy,IsScatProjToProjCase); >> 221 } 196 ////////////////////////////////////////////// 222 //////////////////////////////////////////////////////////////////////////////// 197 G4double G4AdjointPhotoElectricModel::AdjointC << 223 // 198 const G4Element* anElement, G4double electro << 224 199 { << 225 G4double G4AdjointPhotoElectricModel::AdjointCrossSectionPerAtom(const G4Element* anElement,G4double electronEnergy) 200 G4int nShells = anElement->GetNbOfAto << 226 { 201 G4double Z = anElement->GetZ(); << 227 G4int nShells = anElement->GetNbOfAtomicShells(); 202 G4double gammaEnergy = electronEnergy + anEl << 228 G4double Z= anElement->GetZ(); 203 G4double CS = fDirectModel->Compute << 229 G4int i = 0; 204 G4Gamma::Gamma(), gammaEnergy, Z, 0., 0., << 230 G4double B0=anElement->GetAtomicShell(0); 205 G4double adjointCS = 0.; << 231 G4double gammaEnergy = electronEnergy+B0; 206 if(CS > 0.) << 232 G4double CS= theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,Z,0.,0.,0.); 207 adjointCS += CS / gammaEnergy; << 233 G4double adjointCS =0.; 208 fShellProb[fIndexElement][0] = adjointCS; << 234 if (CS >0) adjointCS += CS/gammaEnergy; 209 for(G4int i = 1; i < nShells; ++i) << 235 shell_prob[index_element][0] = adjointCS; 210 { << 236 for (i=1;i<nShells;i++){ 211 G4double Bi1 = anElement->GetAtomicShell(i << 237 //G4cout<<i<<G4endl; 212 G4double Bi = anElement->GetAtomicShell(i << 238 G4double Bi_= anElement->GetAtomicShell(i-1); 213 if(electronEnergy < Bi1 - Bi) << 239 G4double Bi = anElement->GetAtomicShell(i); 214 { << 240 //G4cout<<Bi_<<'\t'<<Bi<<G4endl; 215 gammaEnergy = electronEnergy + Bi; << 241 if (electronEnergy <Bi_-Bi) { 216 CS = fDirectModel->ComputeCross << 242 gammaEnergy = electronEnergy+Bi; 217 << 243 218 if(CS > 0.) << 244 CS=theDirectPEEffectModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),gammaEnergy,Z,0.,0.,0.); 219 adjointCS += CS / gammaEnergy; << 245 if (CS>0) adjointCS +=CS/gammaEnergy; 220 } << 246 } 221 fShellProb[fIndexElement][i] = adjointCS; << 247 shell_prob[index_element][i] = adjointCS; >> 248 222 } 249 } 223 adjointCS *= electronEnergy; << 250 adjointCS*=electronEnergy; 224 return adjointCS; 251 return adjointCS; 225 } << 252 226 << 253 } 227 ////////////////////////////////////////////// 254 //////////////////////////////////////////////////////////////////////////////// 228 void G4AdjointPhotoElectricModel::DefineCurren << 255 // 229 const G4MaterialCutsCouple* couple, G4double << 256 230 { << 257 void G4AdjointPhotoElectricModel::DefineCurrentMaterialAndElectronEnergy(const G4MaterialCutsCouple* couple, G4double anEnergy) 231 fCurrentCouple = const_cast<G4MaterialCuts << 258 { currentCouple = const_cast<G4MaterialCutsCouple*> (couple); 232 fCurrentMaterial = const_cast<G4Material*>(c << 259 currentMaterial = const_cast<G4Material*> (couple->GetMaterial()); 233 fCurrenteEnergy = anEnergy; << 260 currentCoupleIndex = couple->GetIndex(); 234 fDirectModel->SetCurrentCouple(couple); << 261 currentMaterialIndex = currentMaterial->GetIndex(); >> 262 current_eEnergy = anEnergy; 235 } 263 } 236 264