#include <G4AdjointComptonModel.hh>
Inheritance diagram for G4AdjointComptonModel:
Definition at line 54 of file G4AdjointComptonModel.hh.
| G4AdjointComptonModel::G4AdjointComptonModel | ( | ) |
Definition at line 44 of file G4AdjointComptonModel.cc.
References G4AdjointElectron::AdjointElectron(), G4AdjointGamma::AdjointGamma(), G4Gamma::Gamma(), G4VEmAdjointModel::second_part_of_same_type, G4VEmAdjointModel::SetApplyCutInRange(), G4VEmAdjointModel::SetUseMatrix(), G4VEmAdjointModel::SetUseMatrixPerElement(), G4VEmAdjointModel::SetUseOnlyOneMatrixForAllElements(), G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef, G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.
00044 : 00045 G4VEmAdjointModel("AdjointCompton") 00046 00047 { SetApplyCutInRange(false); 00048 SetUseMatrix(false); 00049 SetUseMatrixPerElement(true); 00050 SetUseOnlyOneMatrixForAllElements(true); 00051 theAdjEquivOfDirectPrimPartDef =G4AdjointGamma::AdjointGamma(); 00052 theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron(); 00053 theDirectPrimaryPartDef=G4Gamma::Gamma(); 00054 second_part_of_same_type=false; 00055 theDirectEMModel=new G4KleinNishinaCompton(G4Gamma::Gamma(),"ComptonDirectModel"); 00056 G4direct_CS = 0.; 00057 }
| G4AdjointComptonModel::~G4AdjointComptonModel | ( | ) |
| G4double G4AdjointComptonModel::AdjointCrossSection | ( | const G4MaterialCutsCouple * | aCouple, | |
| G4double | primEnergy, | |||
| G4bool | IsScatProjToProjCase | |||
| ) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 376 of file G4AdjointComptonModel.cc.
References G4VEmAdjointModel::AdjointCrossSection(), G4VEmAdjointModel::currentMaterial, G4VEmAdjointModel::DefineCurrentMaterial(), G4Material::GetElectronDensity(), G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(), GetSecondAdjEnergyMaxForScatProjToProjCase(), GetSecondAdjEnergyMinForProdToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(), G4VEmAdjointModel::lastCS, and G4VEmAdjointModel::UseMatrix.
Referenced by GetAdjointCrossSection().
00379 { 00380 if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase); 00381 DefineCurrentMaterial(aCouple); 00382 00383 00384 float Cross=0.; 00385 float Emax_proj =0.; 00386 float Emin_proj =0.; 00387 if (!IsScatProjToProjCase ){ 00388 Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy); 00389 Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy); 00390 if (Emax_proj>Emin_proj ){ 00391 Cross= std::log((Emax_proj-float (primEnergy))*Emin_proj/Emax_proj/(Emin_proj-primEnergy)) 00392 *(1.+2.*std::log(float(1.+electron_mass_c2/primEnergy))); 00393 } 00394 } 00395 else { 00396 Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy); 00397 Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,0.); 00398 if (Emax_proj>Emin_proj) { 00399 Cross = std::log(Emax_proj/Emin_proj); 00400 //+0.5*primEnergy*primEnergy(1./(Emin_proj*Emin_proj) - 1./(Emax_proj*Emax_proj)); neglected at the moment 00401 } 00402 00403 00404 } 00405 00406 Cross*=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2; 00407 lastCS=Cross; 00408 return double(Cross); 00409 }
| G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim | ( | G4double | kinEnergyProj, | |
| G4double | kinEnergyScatProj, | |||
| G4double | Z, | |||
| G4double | A = 0. | |||
| ) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 293 of file G4AdjointComptonModel.cc.
References G4VEmModel::ComputeCrossSectionPerAtom(), G4Gamma::Gamma(), and G4VEmAdjointModel::theDirectEMModel.
Referenced by DiffCrossSectionPerAtomPrimToSecond(), and RapidSampleSecondaries().
00298 { //Based on Klein Nishina formula 00299 // In the forward case (see G4KleinNishinaModel) the cross section is parametrised while 00300 // the secondaries are sampled from the 00301 // Klein Nishida differential cross section 00302 // The used diffrential cross section here is therefore the cross section multiplied by the normalised 00303 //differential Klein Nishida cross section 00304 00305 00306 //Klein Nishida Cross Section 00307 //----------------------------- 00308 G4double epsilon = gamEnergy0 / electron_mass_c2 ; 00309 G4double one_plus_two_epsi =1.+2.*epsilon; 00310 G4double gamEnergy1_max = gamEnergy0; 00311 G4double gamEnergy1_min = gamEnergy0/one_plus_two_epsi; 00312 if (gamEnergy1 >gamEnergy1_max || gamEnergy1<gamEnergy1_min) { 00313 /*G4cout<<"the differential CS is null"<<G4endl; 00314 G4cout<<gamEnergy0<<G4endl; 00315 G4cout<<gamEnergy1<<G4endl; 00316 G4cout<<gamEnergy1_min<<G4endl;*/ 00317 return 0.; 00318 } 00319 00320 00321 G4double epsi2 = epsilon *epsilon ; 00322 G4double one_plus_two_epsi_2=one_plus_two_epsi*one_plus_two_epsi; 00323 00324 00325 G4double CS=std::log(one_plus_two_epsi)*(1.- 2.*(1.+epsilon)/epsi2); 00326 CS+=4./epsilon +0.5*(1.-1./one_plus_two_epsi_2); 00327 CS/=epsilon; 00328 //Note that the pi*re2*Z factor is neglected because it is suppresed when computing dCS_dE1/CS; 00329 // in the differential cross section 00330 00331 00332 //Klein Nishida Differential Cross Section 00333 //----------------------------------------- 00334 G4double epsilon1 = gamEnergy1 / electron_mass_c2 ; 00335 G4double v= epsilon1/epsilon; 00336 G4double term1 =1.+ 1./epsilon -1/epsilon1; 00337 G4double dCS_dE1= 1./v +v + term1*term1 -1.; 00338 dCS_dE1 *=1./epsilon/gamEnergy0; 00339 00340 00341 //Normalised to the CS used in G4 00342 //------------------------------- 00343 00344 G4direct_CS = theDirectEMModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(), 00345 gamEnergy0, 00346 Z, 0., 0.,0.); 00347 00348 dCS_dE1 *= G4direct_CS/CS; 00349 /* G4cout<<"the differential CS is not null"<<G4endl; 00350 G4cout<<gamEnergy0<<G4endl; 00351 G4cout<<gamEnergy1<<G4endl;*/ 00352 00353 return dCS_dE1; 00354 00355 00356 }
| G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToSecond | ( | G4double | kinEnergyProj, | |
| G4double | kinEnergyProd, | |||
| G4double | Z, | |||
| G4double | A = 0. | |||
| ) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 278 of file G4AdjointComptonModel.cc.
References DiffCrossSectionPerAtomPrimToScatPrim().
00283 { 00284 G4double gamEnergy1 = gamEnergy0 - kinEnergyElec; 00285 G4double dSigmadEprod=0.; 00286 if (gamEnergy1>0.) dSigmadEprod=DiffCrossSectionPerAtomPrimToScatPrim(gamEnergy0,gamEnergy1,Z,A); 00287 return dSigmadEprod; 00288 }
| G4double G4AdjointComptonModel::GetAdjointCrossSection | ( | const G4MaterialCutsCouple * | aCouple, | |
| G4double | primEnergy, | |||
| G4bool | IsScatProjToProjCase | |||
| ) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 412 of file G4AdjointComptonModel.cc.
References AdjointCrossSection().
00415 { return AdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase); 00416 //return G4VEmAdjointModel::GetAdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase); 00417 }
| G4double G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase | ( | G4double | PrimAdjEnergy | ) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 360 of file G4AdjointComptonModel.cc.
References G4VEmAdjointModel::HighEnergyLimit.
Referenced by AdjointCrossSection(), and RapidSampleSecondaries().
00361 { G4double inv_e_max = 1./PrimAdjEnergy - 2./electron_mass_c2; 00362 G4double e_max = HighEnergyLimit; 00363 if (inv_e_max > 0. ) e_max=std::min(1./inv_e_max,HighEnergyLimit); 00364 return e_max; 00365 }
| G4double G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase | ( | G4double | PrimAdjEnergy | ) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 368 of file G4AdjointComptonModel.cc.
Referenced by AdjointCrossSection(), and RapidSampleSecondaries().
00369 { G4double half_e=PrimAdjEnergy/2.; 00370 G4double term=std::sqrt(half_e*(electron_mass_c2+half_e)); 00371 G4double emin=half_e+term; 00372 return emin; 00373 }
| void G4AdjointComptonModel::RapidSampleSecondaries | ( | const G4Track & | aTrack, | |
| G4bool | IsScatProjToProjCase, | |||
| G4ParticleChange * | fParticleChange | |||
| ) |
Definition at line 156 of file G4AdjointComptonModel.cc.
References G4ParticleChange::AddSecondary(), G4VEmAdjointModel::currentCouple, G4VEmAdjointModel::currentMaterial, G4VEmAdjointModel::currentTcutForDirectSecond, G4VEmAdjointModel::DefineCurrentMaterial(), DiffCrossSectionPerAtomPrimToScatPrim(), fStopAndKill, G4UniformRand, G4AdjointCSManager::GetAdjointCSManager(), G4Track::GetDynamicParticle(), G4Material::GetElectronDensity(), G4DynamicParticle::GetKineticEnergy(), G4VEmProcess::GetLambda(), G4Track::GetMaterialCutsCouple(), G4DynamicParticle::GetMomentumDirection(), G4AdjointCSManager::GetPostStepWeightCorrection(), G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(), GetSecondAdjEnergyMaxForScatProjToProjCase(), GetSecondAdjEnergyMinForProdToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(), G4DynamicParticle::GetTotalMomentum(), G4Track::GetWeight(), G4VEmAdjointModel::HighEnergyLimit, G4ParticleChange::ProposeEnergy(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeParentWeight(), G4VParticleChange::ProposeTrackStatus(), G4VParticleChange::SetParentWeightByProcess(), G4VParticleChange::SetSecondaryWeightByProcess(), and G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef.
Referenced by SampleSecondaries().
00159 { 00160 00161 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); 00162 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple()); 00163 00164 00165 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); 00166 00167 00168 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ 00169 return; 00170 } 00171 00172 00173 00174 G4double diffCSUsed=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2; 00175 G4double gammaE1=0.; 00176 G4double gammaE2=0.; 00177 if (!IsScatProjToProjCase){ 00178 00179 G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy); 00180 G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);; 00181 if (Emin>=Emax) return; 00182 G4double f1=(Emin-adjointPrimKinEnergy)/Emin; 00183 G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1; 00184 gammaE1=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand()));; 00185 gammaE2=gammaE1-adjointPrimKinEnergy; 00186 diffCSUsed= diffCSUsed*(1.+2.*std::log(1.+electron_mass_c2/adjointPrimKinEnergy))*adjointPrimKinEnergy/gammaE1/gammaE2; 00187 00188 00189 } 00190 else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy); 00191 G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond); 00192 if (Emin>=Emax) return; 00193 gammaE2 =adjointPrimKinEnergy; 00194 gammaE1=Emin*std::pow(Emax/Emin,G4UniformRand()); 00195 diffCSUsed= diffCSUsed/gammaE1; 00196 00197 } 00198 00199 00200 00201 00202 //Weight correction 00203 //----------------------- 00204 //First w_corr is set to the ratio between adjoint total CS and fwd total CS 00205 G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection(); 00206 00207 //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the 00208 //one consistent with the direct model 00209 00210 00211 G4double diffCS = DiffCrossSectionPerAtomPrimToScatPrim(gammaE1, gammaE2,1,0.); 00212 if (diffCS >0) diffCS /=G4direct_CS; // here we have the normalised diffCS 00213 diffCS*=theDirectEMProcess->GetLambda(gammaE1,currentCouple); 00214 //diffCS*=theDirectEMModel->CrossSectionPerVolume(currentMaterial,G4Gamma::Gamma(),gammaE1,0.,2.*gammaE1); 00215 //G4cout<<"diffCS/diffCSUsed "<<diffCS/diffCSUsed<<'\t'<<gammaE1<<'\t'<<gammaE2<<G4endl; 00216 00217 w_corr*=diffCS/diffCSUsed; 00218 00219 G4double new_weight = aTrack.GetWeight()*w_corr; 00220 fParticleChange->SetParentWeightByProcess(false); 00221 fParticleChange->SetSecondaryWeightByProcess(false); 00222 fParticleChange->ProposeParentWeight(new_weight); 00223 00224 00225 00226 //Cos th 00227 //------- 00228 00229 G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2); 00230 if (!IsScatProjToProjCase) { 00231 G4double p_elec=theAdjointPrimary->GetTotalMomentum(); 00232 cos_th = (gammaE1 - gammaE2*cos_th)/p_elec; 00233 } 00234 G4double sin_th = 0.; 00235 if (std::abs(cos_th)>1){ 00236 //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl; 00237 if (cos_th>0) { 00238 cos_th=1.; 00239 } 00240 else cos_th=-1.; 00241 sin_th=0.; 00242 } 00243 else sin_th = std::sqrt(1.-cos_th*cos_th); 00244 00245 00246 00247 00248 //gamma0 momentum 00249 //-------------------- 00250 00251 00252 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection(); 00253 G4double phi =G4UniformRand()*2.*3.1415926; 00254 G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th); 00255 gammaMomentum1.rotateUz(dir_parallel); 00256 00257 00258 00259 00260 if (!IsScatProjToProjCase){ //kill the primary and add a secondary 00261 fParticleChange->ProposeTrackStatus(fStopAndKill); 00262 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1)); 00263 //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl; 00264 } 00265 else { 00266 fParticleChange->ProposeEnergy(gammaE1); 00267 fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit()); 00268 } 00269 00270 00271 00272 }
| void G4AdjointComptonModel::SampleSecondaries | ( | const G4Track & | aTrack, | |
| G4bool | IsScatProjToProjCase, | |||
| G4ParticleChange * | fParticleChange | |||
| ) | [virtual] |
Implements G4VEmAdjointModel.
Definition at line 64 of file G4AdjointComptonModel.cc.
References G4ParticleChange::AddSecondary(), G4VEmAdjointModel::CorrectPostStepWeight(), fStopAndKill, G4UniformRand, G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), G4DynamicParticle::GetTotalMomentum(), G4VParticleChange::GetWeight(), G4VEmAdjointModel::HighEnergyLimit, G4ParticleChange::ProposeEnergy(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), RapidSampleSecondaries(), G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(), G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, and G4VEmAdjointModel::UseMatrix.
00067 { 00068 if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange); 00069 00070 //A recall of the compton scattering law is 00071 //Egamma2=Egamma1/(1+(Egamma1/E0_electron)(1.-cos_th)) 00072 //Therefore Egamma2_max= Egamma2(cos_th=1) = Egamma1 00073 //Therefore Egamma2_min= Egamma2(cos_th=-1) = Egamma1/(1+2.(Egamma1/E0_electron)) 00074 00075 00076 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); 00077 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); 00078 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ 00079 return; 00080 } 00081 00082 00083 //Sample secondary energy 00084 //----------------------- 00085 G4double gammaE1; 00086 gammaE1 = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, 00087 IsScatProjToProjCase); 00088 00089 00090 //gammaE2 00091 //----------- 00092 00093 G4double gammaE2 = adjointPrimKinEnergy; 00094 if (!IsScatProjToProjCase) gammaE2 = gammaE1 - adjointPrimKinEnergy; 00095 00096 00097 00098 00099 00100 00101 //Cos th 00102 //------- 00103 // G4cout<<"Compton scattering "<<gammaE1<<'\t'<<gammaE2<<G4endl; 00104 G4double cos_th = 1.+ electron_mass_c2*(1./gammaE1 -1./gammaE2); 00105 if (!IsScatProjToProjCase) { 00106 G4double p_elec=theAdjointPrimary->GetTotalMomentum(); 00107 cos_th = (gammaE1 - gammaE2*cos_th)/p_elec; 00108 } 00109 G4double sin_th = 0.; 00110 if (std::abs(cos_th)>1){ 00111 //G4cout<<"Problem in compton scattering with cos_th "<<cos_th<<G4endl; 00112 if (cos_th>0) { 00113 cos_th=1.; 00114 } 00115 else cos_th=-1.; 00116 sin_th=0.; 00117 } 00118 else sin_th = std::sqrt(1.-cos_th*cos_th); 00119 00120 00121 00122 00123 //gamma0 momentum 00124 //-------------------- 00125 00126 00127 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection(); 00128 G4double phi =G4UniformRand()*2.*3.1415926; 00129 G4ThreeVector gammaMomentum1 = gammaE1*G4ThreeVector(std::cos(phi)*sin_th,std::sin(phi)*sin_th,cos_th); 00130 gammaMomentum1.rotateUz(dir_parallel); 00131 // G4cout<<gamma0Energy<<'\t'<<gamma0Momentum<<G4endl; 00132 00133 00134 //It is important to correct the weight of particles before adding the secondary 00135 //------------------------------------------------------------------------------ 00136 CorrectPostStepWeight(fParticleChange, 00137 aTrack.GetWeight(), 00138 adjointPrimKinEnergy, 00139 gammaE1, 00140 IsScatProjToProjCase); 00141 00142 if (!IsScatProjToProjCase){ //kill the primary and add a secondary 00143 fParticleChange->ProposeTrackStatus(fStopAndKill); 00144 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1)); 00145 //G4cout<<"gamma0Momentum "<<gamma0Momentum<<G4endl; 00146 } 00147 else { 00148 fParticleChange->ProposeEnergy(gammaE1); 00149 fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit()); 00150 } 00151 00152 00153 }
| void G4AdjointComptonModel::SetDirectProcess | ( | G4VEmProcess * | aProcess | ) | [inline] |
1.4.7