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 "G4AdjointeIonisationModel.hh"
00029 #include "G4AdjointCSManager.hh"
00030
00031 #include "G4PhysicalConstants.hh"
00032 #include "G4Integrator.hh"
00033 #include "G4TrackStatus.hh"
00034 #include "G4ParticleChange.hh"
00035 #include "G4AdjointElectron.hh"
00036 #include "G4Gamma.hh"
00037 #include "G4AdjointGamma.hh"
00038
00039
00041
00042 G4AdjointeIonisationModel::G4AdjointeIonisationModel():
00043 G4VEmAdjointModel("Inv_eIon_model")
00044
00045 {
00046
00047 UseMatrix =true;
00048 UseMatrixPerElement = true;
00049 ApplyCutInRange = true;
00050 UseOnlyOneMatrixForAllElements = true;
00051 CS_biasing_factor =1.;
00052 WithRapidSampling = false;
00053
00054 theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron();
00055 theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron();
00056 theDirectPrimaryPartDef=G4Electron::Electron();
00057 second_part_of_same_type=true;
00058 }
00060
00061 G4AdjointeIonisationModel::~G4AdjointeIonisationModel()
00062 {;}
00064
00065 void G4AdjointeIonisationModel::SampleSecondaries(const G4Track& aTrack,
00066 G4bool IsScatProjToProjCase,
00067 G4ParticleChange* fParticleChange)
00068 {
00069
00070
00071 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle();
00072
00073
00074
00075 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
00076 G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum();
00077
00078 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){
00079 return;
00080 }
00081
00082
00083
00084 G4double projectileKinEnergy;
00085 if (!WithRapidSampling ) {
00086 projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase);
00087
00088 CorrectPostStepWeight(fParticleChange,
00089 aTrack.GetWeight(),
00090 adjointPrimKinEnergy,
00091 projectileKinEnergy,
00092 IsScatProjToProjCase);
00093 }
00094 else {
00095
00096 G4double Emin,Emax;
00097 if (IsScatProjToProjCase) {
00098 Emin=GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond);
00099 Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy);
00100 }
00101 else {
00102 Emin=GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);
00103 Emax=GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy);
00104 }
00105 projectileKinEnergy = Emin*std::pow(Emax/Emin,G4UniformRand());
00106
00107
00108
00109 lastCS=lastAdjointCSForScatProjToProjCase;
00110 if ( !IsScatProjToProjCase) lastCS=lastAdjointCSForProdToProjCase;
00111
00112 G4double new_weight=aTrack.GetWeight();
00113 G4double used_diffCS=lastCS*std::log(Emax/Emin)/projectileKinEnergy;
00114 G4double needed_diffCS=adjointPrimKinEnergy/projectileKinEnergy;
00115 if (!IsScatProjToProjCase) needed_diffCS *=DiffCrossSectionPerVolumePrimToSecond(currentMaterial,projectileKinEnergy,adjointPrimKinEnergy);
00116 else needed_diffCS *=DiffCrossSectionPerVolumePrimToScatPrim(currentMaterial,projectileKinEnergy,adjointPrimKinEnergy);
00117 new_weight*=needed_diffCS/used_diffCS;
00118 fParticleChange->SetParentWeightByProcess(false);
00119 fParticleChange->SetSecondaryWeightByProcess(false);
00120 fParticleChange->ProposeParentWeight(new_weight);
00121
00122
00123 }
00124
00125
00126
00127
00128
00129
00130
00131
00132 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00133 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
00134 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
00135
00136
00137
00138
00139
00140 G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass();
00141 if (IsScatProjToProjCase) {
00142 companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass();
00143 }
00144 G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
00145 G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
00146
00147
00148
00149
00150 G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
00151 G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
00152 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection();
00153 G4double phi =G4UniformRand()*2.*3.1415926;
00154 G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel);
00155 projectileMomentum.rotateUz(dir_parallel);
00156
00157
00158
00159 if (!IsScatProjToProjCase ){
00160 fParticleChange->ProposeTrackStatus(fStopAndKill);
00161 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum));
00162
00163 }
00164 else {
00165 fParticleChange->ProposeEnergy(projectileKinEnergy);
00166 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
00167 }
00168
00169
00170
00171
00172 }
00174
00175
00176 G4double G4AdjointeIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
00177 G4double kinEnergyProj,
00178 G4double kinEnergyProd,
00179 G4double Z,
00180 G4double )
00181 {
00182 G4double dSigmadEprod=0;
00183 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
00184 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
00185
00186
00187 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
00188 dSigmadEprod=Z*DiffCrossSectionMoller(kinEnergyProj,kinEnergyProd);
00189 }
00190 return dSigmadEprod;
00191
00192
00193
00194 }
00195
00197
00198 G4double G4AdjointeIonisationModel::DiffCrossSectionMoller(G4double kinEnergyProj,G4double kinEnergyProd){
00199
00200 G4double energy = kinEnergyProj + electron_mass_c2;
00201 G4double x = kinEnergyProd/kinEnergyProj;
00202 G4double gam = energy/electron_mass_c2;
00203 G4double gamma2 = gam*gam;
00204 G4double beta2 = 1.0 - 1.0/gamma2;
00205
00206 G4double gg = (2.0*gam - 1.0)/gamma2;
00207 G4double y = 1.0 - x;
00208 G4double fac=twopi_mc2_rcl2/electron_mass_c2;
00209 G4double dCS = fac*( 1.-gg + ((1.0 - gg*x)/(x*x))
00210 + ((1.0 - gg*y)/(y*y)))/(beta2*(gam-1));
00211 return dCS/kinEnergyProj;
00212
00213
00214
00215 }
00216