00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include "G4AdjointComptonModel.hh"
00029 #include "G4AdjointCSManager.hh"
00030
00031 #include "G4PhysicalConstants.hh"
00032 #include "G4SystemOfUnits.hh"
00033 #include "G4Integrator.hh"
00034 #include "G4TrackStatus.hh"
00035 #include "G4ParticleChange.hh"
00036 #include "G4AdjointElectron.hh"
00037 #include "G4AdjointGamma.hh"
00038 #include "G4Gamma.hh"
00039 #include "G4KleinNishinaCompton.hh"
00040
00041
00043
00044 G4AdjointComptonModel::G4AdjointComptonModel():
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 }
00059
00060 G4AdjointComptonModel::~G4AdjointComptonModel()
00061 {;}
00063
00064 void G4AdjointComptonModel::SampleSecondaries(const G4Track& aTrack,
00065 G4bool IsScatProjToProjCase,
00066 G4ParticleChange* fParticleChange)
00067 {
00068 if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange);
00069
00070
00071
00072
00073
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
00084
00085 G4double gammaE1;
00086 gammaE1 = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy,
00087 IsScatProjToProjCase);
00088
00089
00090
00091
00092
00093 G4double gammaE2 = adjointPrimKinEnergy;
00094 if (!IsScatProjToProjCase) gammaE2 = gammaE1 - adjointPrimKinEnergy;
00095
00096
00097
00098
00099
00100
00101
00102
00103
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
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
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
00132
00133
00134
00135
00136 CorrectPostStepWeight(fParticleChange,
00137 aTrack.GetWeight(),
00138 adjointPrimKinEnergy,
00139 gammaE1,
00140 IsScatProjToProjCase);
00141
00142 if (!IsScatProjToProjCase){
00143 fParticleChange->ProposeTrackStatus(fStopAndKill);
00144 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
00145
00146 }
00147 else {
00148 fParticleChange->ProposeEnergy(gammaE1);
00149 fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
00150 }
00151
00152
00153 }
00155
00156 void G4AdjointComptonModel::RapidSampleSecondaries(const G4Track& aTrack,
00157 G4bool IsScatProjToProjCase,
00158 G4ParticleChange* fParticleChange)
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
00203
00204
00205 G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
00206
00207
00208
00209
00210
00211 G4double diffCS = DiffCrossSectionPerAtomPrimToScatPrim(gammaE1, gammaE2,1,0.);
00212 if (diffCS >0) diffCS /=G4direct_CS;
00213 diffCS*=theDirectEMProcess->GetLambda(gammaE1,currentCouple);
00214
00215
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
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
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
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){
00261 fParticleChange->ProposeTrackStatus(fStopAndKill);
00262 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,gammaMomentum1));
00263
00264 }
00265 else {
00266 fParticleChange->ProposeEnergy(gammaE1);
00267 fParticleChange->ProposeMomentumDirection(gammaMomentum1.unit());
00268 }
00269
00270
00271
00272 }
00273
00274
00276
00277
00278 G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToSecond(
00279 G4double gamEnergy0,
00280 G4double kinEnergyElec,
00281 G4double Z,
00282 G4double A)
00283 {
00284 G4double gamEnergy1 = gamEnergy0 - kinEnergyElec;
00285 G4double dSigmadEprod=0.;
00286 if (gamEnergy1>0.) dSigmadEprod=DiffCrossSectionPerAtomPrimToScatPrim(gamEnergy0,gamEnergy1,Z,A);
00287 return dSigmadEprod;
00288 }
00289
00290
00292
00293 G4double G4AdjointComptonModel::DiffCrossSectionPerAtomPrimToScatPrim(
00294 G4double gamEnergy0,
00295 G4double gamEnergy1,
00296 G4double Z,
00297 G4double )
00298 {
00299
00300
00301
00302
00303
00304
00305
00306
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
00314
00315
00316
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
00329
00330
00331
00332
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
00342
00343
00344 G4direct_CS = theDirectEMModel->ComputeCrossSectionPerAtom(G4Gamma::Gamma(),
00345 gamEnergy0,
00346 Z, 0., 0.,0.);
00347
00348 dCS_dE1 *= G4direct_CS/CS;
00349
00350
00351
00352
00353 return dCS_dE1;
00354
00355
00356 }
00357
00359
00360 G4double G4AdjointComptonModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
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 }
00367
00368 G4double G4AdjointComptonModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
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 }
00375
00376 G4double G4AdjointComptonModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
00377 G4double primEnergy,
00378 G4bool IsScatProjToProjCase)
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
00401 }
00402
00403
00404 }
00405
00406 Cross*=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2;
00407 lastCS=Cross;
00408 return double(Cross);
00409 }
00411
00412 G4double G4AdjointComptonModel::GetAdjointCrossSection(const G4MaterialCutsCouple* aCouple,
00413 G4double primEnergy,
00414 G4bool IsScatProjToProjCase)
00415 { return AdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase);
00416
00417 }