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
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 #include "G4PenelopeAnnihilationModel.hh"
00039 #include "G4PhysicalConstants.hh"
00040 #include "G4SystemOfUnits.hh"
00041 #include "G4ParticleDefinition.hh"
00042 #include "G4MaterialCutsCouple.hh"
00043 #include "G4ProductionCutsTable.hh"
00044 #include "G4DynamicParticle.hh"
00045 #include "G4Gamma.hh"
00046
00047
00048
00049
00050 G4PenelopeAnnihilationModel::G4PenelopeAnnihilationModel(const G4ParticleDefinition*,
00051 const G4String& nam)
00052 :G4VEmModel(nam),fParticleChange(0),isInitialised(false)
00053 {
00054 fIntrinsicLowEnergyLimit = 0.0;
00055 fIntrinsicHighEnergyLimit = 100.0*GeV;
00056
00057 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
00058
00059
00060 fPielr2 = pi*classic_electr_radius*classic_electr_radius;
00061
00062 verboseLevel= 0;
00063
00064
00065
00066
00067
00068
00069
00070 }
00071
00072
00073
00074 G4PenelopeAnnihilationModel::~G4PenelopeAnnihilationModel()
00075 {;}
00076
00077
00078
00079 void G4PenelopeAnnihilationModel::Initialise(const G4ParticleDefinition*,
00080 const G4DataVector&)
00081 {
00082 if (verboseLevel > 3)
00083 G4cout << "Calling G4PenelopeAnnihilationModel::Initialise()" << G4endl;
00084
00085 if(verboseLevel > 0) {
00086 G4cout << "Penelope Annihilation model is initialized " << G4endl
00087 << "Energy range: "
00088 << LowEnergyLimit() / keV << " keV - "
00089 << HighEnergyLimit() / GeV << " GeV"
00090 << G4endl;
00091 }
00092
00093 if(isInitialised) return;
00094 fParticleChange = GetParticleChangeForGamma();
00095 isInitialised = true;
00096 }
00097
00098
00099
00100 G4double G4PenelopeAnnihilationModel::ComputeCrossSectionPerAtom(
00101 const G4ParticleDefinition*,
00102 G4double energy,
00103 G4double Z, G4double,
00104 G4double, G4double)
00105 {
00106 if (verboseLevel > 3)
00107 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeAnnihilationModel" <<
00108 G4endl;
00109
00110 G4double cs = Z*ComputeCrossSectionPerElectron(energy);
00111
00112 if (verboseLevel > 2)
00113 G4cout << "Annihilation cross Section at " << energy/keV << " keV for Z=" << Z <<
00114 " = " << cs/barn << " barn" << G4endl;
00115 return cs;
00116 }
00117
00118
00119
00120 void G4PenelopeAnnihilationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00121 const G4MaterialCutsCouple*,
00122 const G4DynamicParticle* aDynamicPositron,
00123 G4double,
00124 G4double)
00125 {
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141 if (verboseLevel > 3)
00142 G4cout << "Calling SamplingSecondaries() of G4PenelopeAnnihilationModel" << G4endl;
00143
00144 G4double kineticEnergy = aDynamicPositron->GetKineticEnergy();
00145
00146
00147 fParticleChange->SetProposedKineticEnergy(0.);
00148 fParticleChange->ProposeTrackStatus(fStopAndKill);
00149
00150 if (kineticEnergy == 0.0)
00151 {
00152
00153 G4double cosTheta = -1.0+2.0*G4UniformRand();
00154 G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta);
00155 G4double phi = twopi*G4UniformRand();
00156 G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta);
00157 G4DynamicParticle* firstGamma = new G4DynamicParticle (G4Gamma::Gamma(),
00158 direction, electron_mass_c2);
00159 G4DynamicParticle* secondGamma = new G4DynamicParticle (G4Gamma::Gamma(),
00160 -direction, electron_mass_c2);
00161
00162 fvect->push_back(firstGamma);
00163 fvect->push_back(secondGamma);
00164 return;
00165 }
00166
00167
00168 G4ParticleMomentum positronDirection =
00169 aDynamicPositron->GetMomentumDirection();
00170 G4double gamma = 1.0 + std::max(kineticEnergy,1.0*eV)/electron_mass_c2;
00171 G4double gamma21 = std::sqrt(gamma*gamma-1);
00172 G4double ani = 1.0+gamma;
00173 G4double chimin = 1.0/(ani+gamma21);
00174 G4double rchi = (1.0-chimin)/chimin;
00175 G4double gt0 = ani*ani-2.0;
00176 G4double test=0.0;
00177 G4double epsilon = 0;
00178 do{
00179 epsilon = chimin*std::pow(rchi,G4UniformRand());
00180 G4double reject = ani*ani*(1.0-epsilon)+2.0*gamma-(1.0/epsilon);
00181 test = G4UniformRand()*gt0-reject;
00182 }while(test>0);
00183
00184 G4double totalAvailableEnergy = kineticEnergy + 2.0*electron_mass_c2;
00185 G4double photon1Energy = epsilon*totalAvailableEnergy;
00186 G4double photon2Energy = (1.0-epsilon)*totalAvailableEnergy;
00187 G4double cosTheta1 = (ani-1.0/epsilon)/gamma21;
00188 G4double cosTheta2 = (ani-1.0/(1.0-epsilon))/gamma21;
00189
00190
00191
00192 G4double sinTheta1 = std::sqrt(1.-cosTheta1*cosTheta1);
00193 G4double phi1 = twopi * G4UniformRand();
00194 G4double dirx1 = sinTheta1 * std::cos(phi1);
00195 G4double diry1 = sinTheta1 * std::sin(phi1);
00196 G4double dirz1 = cosTheta1;
00197
00198 G4double sinTheta2 = std::sqrt(1.-cosTheta2*cosTheta2);
00199 G4double phi2 = phi1+pi;
00200 G4double dirx2 = sinTheta2 * std::cos(phi2);
00201 G4double diry2 = sinTheta2 * std::sin(phi2);
00202 G4double dirz2 = cosTheta2;
00203
00204 G4ThreeVector photon1Direction (dirx1,diry1,dirz1);
00205 photon1Direction.rotateUz(positronDirection);
00206
00207 G4DynamicParticle* aParticle1= new G4DynamicParticle (G4Gamma::Gamma(),
00208 photon1Direction,
00209 photon1Energy);
00210 fvect->push_back(aParticle1);
00211
00212 G4ThreeVector photon2Direction(dirx2,diry2,dirz2);
00213 photon2Direction.rotateUz(positronDirection);
00214
00215 G4DynamicParticle* aParticle2= new G4DynamicParticle (G4Gamma::Gamma(),
00216 photon2Direction,
00217 photon2Energy);
00218 fvect->push_back(aParticle2);
00219
00220 if (verboseLevel > 1)
00221 {
00222 G4cout << "-----------------------------------------------------------" << G4endl;
00223 G4cout << "Energy balance from G4PenelopeAnnihilation" << G4endl;
00224 G4cout << "Kinetic positron energy: " << kineticEnergy/keV << " keV" << G4endl;
00225 G4cout << "Total available energy: " << totalAvailableEnergy/keV << " keV " << G4endl;
00226 G4cout << "-----------------------------------------------------------" << G4endl;
00227 G4cout << "Photon energy 1: " << photon1Energy/keV << " keV" << G4endl;
00228 G4cout << "Photon energy 2: " << photon2Energy/keV << " keV" << G4endl;
00229 G4cout << "Total final state: " << (photon1Energy+photon2Energy)/keV <<
00230 " keV" << G4endl;
00231 G4cout << "-----------------------------------------------------------" << G4endl;
00232 }
00233 if (verboseLevel > 0)
00234 {
00235 G4double energyDiff = std::fabs(totalAvailableEnergy-photon1Energy-photon2Energy);
00236 if (energyDiff > 0.05*keV)
00237 G4cout << "Warning from G4PenelopeAnnihilation: problem with energy conservation: " <<
00238 (photon1Energy+photon2Energy)/keV <<
00239 " keV (final) vs. " <<
00240 totalAvailableEnergy/keV << " keV (initial)" << G4endl;
00241 }
00242 return;
00243 }
00244
00245
00246
00247 G4double G4PenelopeAnnihilationModel:: ComputeCrossSectionPerElectron(G4double energy)
00248 {
00249
00250
00251
00252
00253
00254
00255
00256 G4double gamma = 1.0+std::max(energy,1.0*eV)/electron_mass_c2;
00257 G4double gamma2 = gamma*gamma;
00258 G4double f2 = gamma2-1.0;
00259 G4double f1 = std::sqrt(f2);
00260 G4double crossSection = fPielr2*((gamma2+4.0*gamma+1.0)*std::log(gamma+f1)/f2
00261 - (gamma+3.0)/f1)/(gamma+1.0);
00262 return crossSection;
00263 }