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 // $Id: G4AdjointComptonModel.cc,v 1.7 2010-11-11 11:51:56 ldesorgh Exp $ >> 27 // GEANT4 tag $Name: not supported by cvs2svn $ 26 // 28 // 27 << 28 #include "G4AdjointComptonModel.hh" 29 #include "G4AdjointComptonModel.hh" 29 << 30 #include "G4AdjointCSManager.hh" 30 #include "G4AdjointCSManager.hh" >> 31 >> 32 >> 33 #include "G4Integrator.hh" >> 34 #include "G4TrackStatus.hh" >> 35 #include "G4ParticleChange.hh" 31 #include "G4AdjointElectron.hh" 36 #include "G4AdjointElectron.hh" 32 #include "G4AdjointGamma.hh" 37 #include "G4AdjointGamma.hh" 33 #include "G4Gamma.hh" 38 #include "G4Gamma.hh" 34 #include "G4KleinNishinaCompton.hh" 39 #include "G4KleinNishinaCompton.hh" 35 #include "G4ParticleChange.hh" << 40 36 #include "G4PhysicalConstants.hh" << 37 #include "G4SystemOfUnits.hh" << 38 #include "G4TrackStatus.hh" << 39 #include "G4VEmProcess.hh" << 40 41 41 ////////////////////////////////////////////// 42 //////////////////////////////////////////////////////////////////////////////// 42 G4AdjointComptonModel::G4AdjointComptonModel() << 43 // 43 : G4VEmAdjointModel("AdjointCompton") << 44 G4AdjointComptonModel::G4AdjointComptonModel(): 44 { << 45 G4VEmAdjointModel("AdjointCompton") 45 SetApplyCutInRange(false); << 46 >> 47 { SetApplyCutInRange(false); 46 SetUseMatrix(false); 48 SetUseMatrix(false); 47 SetUseMatrixPerElement(true); 49 SetUseMatrixPerElement(true); 48 SetUseOnlyOneMatrixForAllElements(true); 50 SetUseOnlyOneMatrixForAllElements(true); 49 fAdjEquivDirectPrimPart = G4AdjointGamma:: << 51 theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma(); 50 fAdjEquivDirectSecondPart = G4AdjointElectro << 52 theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron(); 51 fDirectPrimaryPart = G4Gamma::Gamma() << 53 theDirectPrimaryPartDef=G4Gamma::Gamma(); 52 fSecondPartSameType = false; << 54 second_part_of_same_type=false; 53 fDirectModel = << 55 theDirectEMModel=new G4KleinNishinaCompton(G4Gamma::Gamma(),"ComptonDirectModel"); 54 new G4KleinNishinaCompton(G4Gamma::Gamma() << 56 55 } 57 } 56 << 57 ////////////////////////////////////////////// 58 //////////////////////////////////////////////////////////////////////////////// 58 G4AdjointComptonModel::~G4AdjointComptonModel( << 59 // 59 << 60 G4AdjointComptonModel::~G4AdjointComptonModel() >> 61 {;} 60 ////////////////////////////////////////////// 62 //////////////////////////////////////////////////////////////////////////////// >> 63 // 61 void G4AdjointComptonModel::SampleSecondaries( 64 void G4AdjointComptonModel::SampleSecondaries(const G4Track& aTrack, 62 << 65 G4bool IsScatProjToProjCase, 63 << 66 G4ParticleChange* fParticleChange) 64 { << 67 { 65 if(!fUseMatrix) << 68 if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange); 66 return RapidSampleSecondaries(aTrack, isSc << 69 67 << 70 //A recall of the compton scattering law is 68 // A recall of the compton scattering law: << 71 //Egamma2=Egamma1/(1+(Egamma1/E0_electron)(1.-cos_th)) 69 // Egamma2=Egamma1/(1+(Egamma1/E0_electron)( << 72 //Therefore Egamma2_max= Egamma2(cos_th=1) = Egamma1 70 // Therefore Egamma2_max= Egamma2(cos_th=1) << 73 //Therefore Egamma2_min= Egamma2(cos_th=-1) = Egamma1/(1+2.(Egamma1/E0_electron)) 71 // and Egamma2_min= Egamma2(cos_th=-1) << 74 72 // Egamma1/(1+2.(Egamma1/E0_elec << 75 73 const G4DynamicParticle* theAdjointPrimary = << 76 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); 74 G4double adjointPrimKinEnergy = theAdjointPr 77 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); 75 if(adjointPrimKinEnergy > GetHighEnergyLimit << 78 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ 76 { << 79 return; 77 return; << 78 } 80 } 79 << 81 80 // Sample secondary energy << 82 81 G4double gammaE1; << 83 //Sample secondary energy 82 gammaE1 = << 84 //----------------------- 83 SampleAdjSecEnergyFromCSMatrix(adjointPrim << 85 G4double gammaE1; 84 << 86 gammaE1 = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, 85 // gammaE2 << 87 IsScatProjToProjCase); 86 G4double gammaE2 = adjointPrimKinEnergy; << 88 87 if(!isScatProjToProj) << 89 88 gammaE2 = gammaE1 - adjointPrimKinEnergy; << 90 //gammaE2 89 << 91 //----------- 90 // Cos th << 92 91 G4double cos_th = 1. + electron_mass_c2 * (1 << 93 G4double gammaE2 = adjointPrimKinEnergy; 92 if(!isScatProjToProj) << 94 if (!IsScatProjToProjCase) gammaE2 = gammaE1 - adjointPrimKinEnergy; 93 { << 95 94 cos_th = << 96 95 (gammaE1 - gammaE2 * cos_th) / theAdjoin << 97 96 } << 98 97 G4double sin_th = 0.; << 99 98 if(std::abs(cos_th) > 1.) << 100 99 { << 101 //Cos th 100 if(cos_th > 0.) << 102 //------- 101 { << 103 // G4cout<<"Compton scattering "<<gammaE1<<'\t'<<gammaE2<<G4endl; 102 cos_th = 1.; << 104 G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2); 103 } << 105 if (!IsScatProjToProjCase) { 104 else << 106 G4double p_elec=theAdjointPrimary->GetTotalMomentum(); 105 cos_th = -1.; << 107 cos_th = (gammaE1 - gammaE2*cos_th)/p_elec; 106 sin_th = 0.; << 108 } 107 } << 109 G4double sin_th = 0.; 108 else << 110 if (std::abs(cos_th)>1){ 109 sin_th = std::sqrt(1. - cos_th * cos_th); << 111 //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl; 110 << 112 if (cos_th>0) { 111 // gamma0 momentum << 113 cos_th=1.; 112 G4ThreeVector dir_parallel = theAdjointPrima << 114 } 113 G4double phi = G4UniformRand() << 115 else cos_th=-1.; 114 G4ThreeVector gammaMomentum1 = << 116 sin_th=0.; 115 gammaE1 * << 117 } 116 G4ThreeVector(std::cos(phi) * sin_th, std: << 118 else sin_th = std::sqrt(1.-cos_th*cos_th); 117 gammaMomentum1.rotateUz(dir_parallel); << 119 118 << 120 119 // correct the weight of particles before ad << 121 120 CorrectPostStepWeight(fParticleChange, aTrac << 122 121 adjointPrimKinEnergy, << 123 //gamma0 momentum 122 << 124 //-------------------- 123 if(!isScatProjToProj) << 125 124 { // kill the primary and add a secondary << 126 125 fParticleChange->ProposeTrackStatus(fStopA << 127 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection(); 126 fParticleChange->AddSecondary( << 128 G4double phi =G4UniformRand()*2.*3.1415926; 127 new G4DynamicParticle(fAdjEquivDirectPri << 129 G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th); 128 } << 130 gammaMomentum1.rotateUz(dir_parallel); 129 else << 131 // G4cout<<gamma0Energy<<'\t'<<gamma0Momentum<<G4endl; 130 { << 132 131 fParticleChange->ProposeEnergy(gammaE1); << 133 132 fParticleChange->ProposeMomentumDirection( << 134 //It is important to correct the weight of particles before adding the secondary 133 } << 135 //------------------------------------------------------------------------------ 134 } << 136 CorrectPostStepWeight(fParticleChange, 135 << 137 aTrack.GetWeight(), >> 138 adjointPrimKinEnergy, >> 139 gammaE1, >> 140 IsScatProjToProjCase); >> 141 >> 142 if (!IsScatProjToProjCase){ //kill the primary and add a secondary >> 143 fParticleChange->ProposeTrackStatus(fStopAndKill); >> 144 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1)); >> 145 //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl; >> 146 } >> 147 else { >> 148 fParticleChange->ProposeEnergy(gammaE1); >> 149 fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit()); >> 150 } >> 151 >> 152 >> 153 } 136 ////////////////////////////////////////////// 154 //////////////////////////////////////////////////////////////////////////////// 137 void G4AdjointComptonModel::RapidSampleSeconda << 155 // 138 const G4Track& aTrack, G4bool isScatProjToPr << 156 void G4AdjointComptonModel::RapidSampleSecondaries(const G4Track& aTrack, 139 G4ParticleChange* fParticleChange) << 157 G4bool IsScatProjToProjCase, 140 { << 158 G4ParticleChange* fParticleChange) 141 const G4DynamicParticle* theAdjointPrimary = << 159 { 142 DefineCurrentMaterial(aTrack.GetMaterialCuts << 160 143 << 161 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); 144 G4double adjointPrimKinEnergy = theAdjointPr << 162 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple()); 145 << 163 146 if(adjointPrimKinEnergy > GetHighEnergyLimit << 164 147 { << 165 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); 148 return; << 166 149 } << 167 150 << 168 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ 151 G4double diffCSUsed = << 169 return; 152 0.1 * fCurrentMaterial->GetElectronDensity << 170 } 153 G4double gammaE1 = 0.; << 171 154 G4double gammaE2 = 0.; << 172 155 if(!isScatProjToProj) << 173 156 { << 174 G4double diffCSUsed=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2; 157 G4double Emax = GetSecondAdjEnergyMaxForPr << 175 G4double gammaE1=0.; 158 G4double Emin = GetSecondAdjEnergyMinForPr << 176 G4double gammaE2=0.; 159 if(Emin >= Emax) << 177 if (!IsScatProjToProjCase){ 160 return; << 178 161 G4double f1 = (Emin - adjointPrimKinEnergy << 179 G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy); 162 G4double f2 = (Emax - adjointPrimKinEnergy << 180 G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);; 163 gammaE1 = adjointPrimKinEnergy / (1. - f1 << 181 if (Emin>=Emax) return; 164 gammaE2 = gammaE1 - adjointPrimKinEnergy; << 182 G4double f1=(Emin-adjointPrimKinEnergy)/Emin; 165 diffCSUsed = << 183 G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1; 166 diffCSUsed * << 184 gammaE1=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand()));; 167 (1. + 2. * std::log(1. + electron_mass_c << 185 gammaE2=gammaE1-adjointPrimKinEnergy; 168 adjointPrimKinEnergy / gammaE1 / gammaE2 << 186 diffCSUsed= diffCSUsed*(1.+2.*std::log(1.+electron_mass_c2/adjointPrimKinEnergy))*adjointPrimKinEnergy/gammaE1/gammaE2; 169 } << 187 170 else << 188 171 { << 189 } 172 G4double Emax = << 190 else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy); 173 GetSecondAdjEnergyMaxForScatProjToProj(a << 191 G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond); 174 G4double Emin = << 192 if (Emin>=Emax) return; 175 GetSecondAdjEnergyMinForScatProjToProj(a << 193 gammaE2 =adjointPrimKinEnergy; 176 if(Emin >= Emax) << 194 gammaE1=Emin*std::pow(Emax/Emin,G4UniformRand()); 177 return; << 195 diffCSUsed= diffCSUsed/gammaE1; 178 gammaE2 = adjointPrimKinEnergy; << 196 179 gammaE1 = Emin * std::pow(Emax / Emin, << 197 } 180 diffCSUsed = diffCSUsed / gammaE1; << 198 181 } << 199 182 << 200 183 // Weight correction << 201 184 // First w_corr is set to the ratio between << 202 //Weight correction 185 G4double w_corr = fOutsideWeightFactor; << 203 //----------------------- 186 if(fInModelWeightCorr) << 204 //First w_corr is set to the ratio between adjoint total CS and fwd total CS 187 { << 205 G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection(); 188 w_corr = << 206 189 G4AdjointCSManager::GetAdjointCSManager( << 207 //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the 190 } << 208 //one consistent with the direct model 191 // Then another correction is needed because << 209 192 // been used rather than the one consistent << 210 193 << 211 G4double diffCS = DiffCrossSectionPerAtomPrimToScatPrim(gammaE1, gammaE2,1,0.); 194 G4double diffCS = << 212 if (diffCS >0) diffCS /=G4direct_CS; // here we have the normalised diffCS 195 DiffCrossSectionPerAtomPrimToScatPrim(gamm << 213 diffCS*=theDirectEMProcess->GetLambda(gammaE1,currentCouple); 196 if(diffCS > 0.) << 214 //diffCS*=theDirectEMModel->CrossSectionPerVolume(currentMaterial,G4Gamma::Gamma(),gammaE1,0.,2.*gammaE1); 197 diffCS /= fDirectCS; // here we have the << 215 //G4cout<<"diffCS/diffCSUsed "<<diffCS/diffCSUsed<<'\t'<<gammaE1<<'\t'<<gammaE2<<G4endl; 198 // And we remultiply by the lambda of the fo << 216 199 diffCS *= fDirectProcess->GetCrossSection(ga << 217 w_corr*=diffCS/diffCSUsed; 200 << 218 201 w_corr *= diffCS / diffCSUsed; << 219 G4double new_weight = aTrack.GetWeight()*w_corr; 202 << 220 fParticleChange->SetParentWeightByProcess(false); 203 G4double new_weight = aTrack.GetWeight() * w << 221 fParticleChange->SetSecondaryWeightByProcess(false); 204 fParticleChange->SetParentWeightByProcess(fa << 222 fParticleChange->ProposeParentWeight(new_weight); 205 fParticleChange->SetSecondaryWeightByProcess << 223 206 fParticleChange->ProposeParentWeight(new_wei << 224 207 << 225 208 G4double cos_th = 1. + electron_mass_c2 * (1 << 226 //Cos th 209 if(!isScatProjToProj) << 227 //------- 210 { << 228 211 G4double p_elec = theAdjointPrimary->GetTo << 229 G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2); 212 cos_th = (gammaE1 - gammaE2 * cos << 230 if (!IsScatProjToProjCase) { 213 } << 231 G4double p_elec=theAdjointPrimary->GetTotalMomentum(); 214 G4double sin_th = 0.; << 232 cos_th = (gammaE1 - gammaE2*cos_th)/p_elec; 215 if(std::abs(cos_th) > 1.) << 233 } 216 { << 234 G4double sin_th = 0.; 217 if(cos_th > 0.) << 235 if (std::abs(cos_th)>1){ 218 { << 236 //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl; 219 cos_th = 1.; << 237 if (cos_th>0) { 220 } << 238 cos_th=1.; 221 else << 239 } 222 cos_th = -1.; << 240 else cos_th=-1.; 223 } << 241 sin_th=0.; 224 else << 242 } 225 sin_th = std::sqrt(1. - cos_th * cos_th); << 243 else sin_th = std::sqrt(1.-cos_th*cos_th); 226 << 244 227 // gamma0 momentum << 245 228 G4ThreeVector dir_parallel = theAdjointPrima << 246 229 G4double phi = G4UniformRand() << 247 230 G4ThreeVector gammaMomentum1 = << 248 //gamma0 momentum 231 gammaE1 * << 249 //-------------------- 232 G4ThreeVector(std::cos(phi) * sin_th, std: << 250 233 gammaMomentum1.rotateUz(dir_parallel); << 251 234 << 252 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection(); 235 if(!isScatProjToProj) << 253 G4double phi =G4UniformRand()*2.*3.1415926; 236 { // kill the primary and add a secondary << 254 G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th); 237 fParticleChange->ProposeTrackStatus(fStopA << 255 gammaMomentum1.rotateUz(dir_parallel); 238 fParticleChange->AddSecondary( << 256 239 new G4DynamicParticle(fAdjEquivDirectPri << 257 240 } << 258 241 else << 259 242 { << 260 if (!IsScatProjToProjCase){ //kill the primary and add a secondary 243 fParticleChange->ProposeEnergy(gammaE1); << 261 fParticleChange->ProposeTrackStatus(fStopAndKill); 244 fParticleChange->ProposeMomentumDirection( << 262 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1)); 245 } << 263 //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl; 246 } << 264 } >> 265 else { >> 266 fParticleChange->ProposeEnergy(gammaE1); >> 267 fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit()); >> 268 } >> 269 >> 270 >> 271 >> 272 } 247 273 >> 274 248 ////////////////////////////////////////////// 275 //////////////////////////////////////////////////////////////////////////////// 249 // The implementation here is correct for ener << 276 // 250 // photoelectric and compton scattering the me << 277 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine 251 G4double G4AdjointComptonModel::DiffCrossSecti 278 G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToSecond( 252 G4double gamEnergy0, G4double kinEnergyElec, << 279 G4double gamEnergy0, 253 { << 280 G4double kinEnergyElec, 254 G4double gamEnergy1 = gamEnergy0 - kinEner << 281 G4double Z, 255 G4double dSigmadEprod = 0.; << 282 G4double A) 256 if(gamEnergy1 > 0.) << 283 { 257 dSigmadEprod = << 284 G4double gamEnergy1 = gamEnergy0 - kinEnergyElec; 258 DiffCrossSectionPerAtomPrimToScatPrim(ga << 285 G4double dSigmadEprod=0.; 259 return dSigmadEprod; << 286 if (gamEnergy1>0.) dSigmadEprod=DiffCrossSectionPerAtomPrimToScatPrim(gamEnergy0,gamEnergy1,Z,A); >> 287 return dSigmadEprod; 260 } 288 } 261 289 >> 290 262 ////////////////////////////////////////////// 291 //////////////////////////////////////////////////////////////////////////////// >> 292 // 263 G4double G4AdjointComptonModel::DiffCrossSecti 293 G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim( 264 G4double gamEnergy0, G4double gamEnergy1, G4 << 294 G4double gamEnergy0, 265 { << 295 G4double gamEnergy1, 266 // Based on Klein Nishina formula << 296 G4double Z, 267 // In the forward case (see G4KleinNishinaCo << 297 G4double ) 268 // parametrised but the secondaries are samp << 298 { //Based on Klein Nishina formula 269 // differential cross section. The different << 299 // In the forward case (see G4KleinNishinaModel) the cross section is parametrised while 270 // is therefore the cross section multiplied << 300 // the secondaries are sampled from the 271 // differential Klein Nishina cross section << 301 // Klein Nishida differential cross section 272 << 302 // The used diffrential cross section here is therefore the cross section multiplied by the normalised 273 // Klein Nishina Cross Section << 303 //differential Klein Nishida cross section 274 G4double epsilon = gamEnergy0 / el << 304 275 G4double one_plus_two_epsi = 1. + 2. * epsil << 305 276 if(gamEnergy1 > gamEnergy0 || gamEnergy1 < g << 306 //Klein Nishida Cross Section 277 { << 307 //----------------------------- 278 return 0.; << 308 G4double epsilon = gamEnergy0 / electron_mass_c2 ; 279 } << 309 G4double one_plus_two_epsi =1.+2.*epsilon; >> 310 G4double gamEnergy1_max = gamEnergy0; >> 311 G4double gamEnergy1_min = gamEnergy0/one_plus_two_epsi; >> 312 if (gamEnergy1 >gamEnergy1_max || gamEnergy1<gamEnergy1_min) { >> 313 /*G4cout<<"the differential CS is null"<<G4endl; >> 314 G4cout<<gamEnergy0<<G4endl; >> 315 G4cout<<gamEnergy1<<G4endl; >> 316 G4cout<<gamEnergy1_min<<G4endl;*/ >> 317 return 0.; >> 318 } >> 319 >> 320 >> 321 G4double epsi2 = epsilon *epsilon ; >> 322 G4double one_plus_two_epsi_2=one_plus_two_epsi*one_plus_two_epsi; >> 323 >> 324 >> 325 G4double CS=std::log(one_plus_two_epsi)*(1.- 2.*(1.+epsilon)/epsi2); >> 326 CS+=4./epsilon +0.5*(1.-1./one_plus_two_epsi_2); >> 327 CS/=epsilon; >> 328 //Note that the pi*re2*Z factor is neglected because it is suppresed when computing dCS_dE1/CS; >> 329 // in the differential cross section >> 330 >> 331 >> 332 //Klein Nishida Differential Cross Section >> 333 //----------------------------------------- >> 334 G4double epsilon1 = gamEnergy1 / electron_mass_c2 ; >> 335 G4double v= epsilon1/epsilon; >> 336 G4double term1 =1.+ 1./epsilon -1/epsilon1; >> 337 G4double dCS_dE1= 1./v +v + term1*term1 -1.; >> 338 dCS_dE1 *=1./epsilon/gamEnergy0; >> 339 >> 340 >> 341 //Normalised to the CS used in G4 >> 342 //------------------------------- >> 343 >> 344 G4direct_CS = theDirectEMModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(), >> 345 gamEnergy0, >> 346 Z, 0., 0.,0.); >> 347 >> 348 dCS_dE1 *= G4direct_CS/CS; >> 349 /* G4cout<<"the differential CS is not null"<<G4endl; >> 350 G4cout<<gamEnergy0<<G4endl; >> 351 G4cout<<gamEnergy1<<G4endl;*/ >> 352 >> 353 return dCS_dE1; 280 354 281 G4double CS = std::log(one_plus_two_epsi) * << 282 (1. - 2. * (1. + epsilon) / (e << 283 CS += << 284 4. / epsilon + 0.5 * (1. - 1. / (one_plus_ << 285 CS /= epsilon; << 286 // Note that the pi*re2*Z factor is neglecte << 287 // computing dCS_dE1/CS in the differential << 288 << 289 // Klein Nishina Differential Cross Section << 290 G4double epsilon1 = gamEnergy1 / electron_ma << 291 G4double v = epsilon1 / epsilon; << 292 G4double term1 = 1. + 1. / epsilon - 1. / << 293 G4double dCS_dE1 = 1. / v + v + term1 * ter << 294 dCS_dE1 *= 1. / epsilon / gamEnergy0; << 295 << 296 // Normalised to the CS used in G4 << 297 fDirectCS = fDirectModel->ComputeCrossSectio << 298 G4Gamma::Gamma(), gamEnergy0, Z, 0., 0., 0 << 299 355 300 dCS_dE1 *= fDirectCS / CS; << 301 << 302 return dCS_dE1; << 303 } 356 } 304 << 357 305 ////////////////////////////////////////////// 358 //////////////////////////////////////////////////////////////////////////////// 306 G4double G4AdjointComptonModel::GetSecondAdjEn << 359 // 307 G4double primAdjEnergy) << 360 G4double G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy) 308 { << 361 { G4double inv_e_max = 1./PrimAdjEnergy - 2./electron_mass_c2; 309 G4double inv_e_max = 1. / primAdjEnergy - 2. << 362 G4double e_max = HighEnergyLimit; 310 G4double e_max = GetHighEnergyLimit(); << 363 if (inv_e_max > 0. ) e_max=std::min(1./inv_e_max,HighEnergyLimit); 311 if(inv_e_max > 0.) << 364 return e_max; 312 e_max = std::min(1. / inv_e_max, e_max); << 313 return e_max; << 314 } 365 } 315 << 316 ////////////////////////////////////////////// 366 //////////////////////////////////////////////////////////////////////////////// 317 G4double G4AdjointComptonModel::GetSecondAdjEn << 367 // 318 G4double primAdjEnergy) << 368 G4double G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy) 319 { << 369 { G4double half_e=PrimAdjEnergy/2.; 320 G4double half_e = primAdjEnergy / 2.; << 370 G4double term=std::sqrt(half_e*(electron_mass_c2+half_e)); 321 return half_e + std::sqrt(half_e * (electron << 371 G4double emin=half_e+term; >> 372 return emin; 322 } 373 } 323 << 324 ////////////////////////////////////////////// 374 //////////////////////////////////////////////////////////////////////////////// 325 G4double G4AdjointComptonModel::AdjointCrossSe << 375 // 326 const G4MaterialCutsCouple* aCouple, G4doubl << 376 G4double G4AdjointComptonModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple, 327 G4bool isScatProjToProj) << 377 G4double primEnergy, 328 { << 378 G4bool IsScatProjToProjCase) 329 if(fUseMatrix) << 379 { 330 return G4VEmAdjointModel::AdjointCrossSect << 380 if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase); 331 << 332 DefineCurrentMaterial(aCouple); 381 DefineCurrentMaterial(aCouple); 333 << 382 334 G4float Cross = 0.; << 383 335 G4float Emax_proj = 0.; << 384 float Cross=0.; 336 G4float Emin_proj = 0.; << 385 float Emax_proj =0.; 337 if(!isScatProjToProj) << 386 float Emin_proj =0.; 338 { << 387 if (!IsScatProjToProjCase ){ 339 Emax_proj = GetSecondAdjEnergyMaxForProdTo << 388 Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy); 340 Emin_proj = GetSecondAdjEnergyMinForProdTo << 389 Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy); 341 if(Emax_proj > Emin_proj) << 390 if (Emax_proj>Emin_proj ){ 342 { << 391 Cross= std::log((Emax_proj-float (primEnergy))*Emin_proj/Emax_proj/(Emin_proj-primEnergy)) 343 Cross = 0.1 * << 392 *(1.+2.*std::log(float(1.+electron_mass_c2/primEnergy))); 344 std::log((Emax_proj - G4float(pr << 393 } 345 Emax_proj / (Emin_proj << 394 } 346 (1. + 2. * std::log(G4float(1. + << 395 else { 347 } << 396 Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy); 348 } << 397 Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,0.); 349 else << 398 if (Emax_proj>Emin_proj) { 350 { << 399 Cross = std::log(Emax_proj/Emin_proj); 351 Emax_proj = GetSecondAdjEnergyMaxForScatPr << 400 //+0.5*primEnergy*primEnergy(1./(Emin_proj*Emin_proj) - 1./(Emax_proj*Emax_proj)); neglected at the moment 352 Emin_proj = GetSecondAdjEnergyMinForScatPr << 401 } 353 if(Emax_proj > Emin_proj) << 402 354 { << 403 355 Cross = 0.1 * std::log(Emax_proj / Emin_ << 404 } 356 } << 405 357 } << 406 Cross*=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2; 358 << 407 lastCS=Cross; 359 Cross *= fCurrentMaterial->GetElectronDensit << 408 return double(Cross); 360 fLastCS = Cross; << 409 } 361 return double(Cross); << 410 //////////////////////////////////////////////////////////////////////////////// >> 411 // >> 412 G4double G4AdjointComptonModel::GetAdjointCrossSection(const G4MaterialCutsCouple* aCouple, >> 413 G4double primEnergy, >> 414 G4bool IsScatProjToProjCase) >> 415 { return AdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase); >> 416 //return G4VEmAdjointModel::GetAdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase); 362 } 417 } 363 418