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
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 #include "G4eCoulombScatteringModel.hh"
00060 #include "G4PhysicalConstants.hh"
00061 #include "G4SystemOfUnits.hh"
00062 #include "Randomize.hh"
00063 #include "G4DataVector.hh"
00064 #include "G4ElementTable.hh"
00065 #include "G4ParticleChangeForGamma.hh"
00066 #include "G4Proton.hh"
00067 #include "G4ParticleTable.hh"
00068 #include "G4ProductionCutsTable.hh"
00069 #include "G4NucleiProperties.hh"
00070 #include "G4Pow.hh"
00071 #include "G4LossTableManager.hh"
00072 #include "G4LossTableBuilder.hh"
00073 #include "G4NistManager.hh"
00074
00075
00076
00077 using namespace std;
00078
00079 G4eCoulombScatteringModel::G4eCoulombScatteringModel(const G4String& nam)
00080 : G4VEmModel(nam),
00081 cosThetaMin(1.0),
00082 cosThetaMax(-1.0),
00083 isInitialised(false)
00084 {
00085 fParticleChange = 0;
00086 fNistManager = G4NistManager::Instance();
00087 theParticleTable = G4ParticleTable::GetParticleTable();
00088 theProton = G4Proton::Proton();
00089 currentMaterial = 0;
00090
00091 pCuts = 0;
00092
00093 lowEnergyThreshold = 1*keV;
00094 recoilThreshold = 0.*keV;
00095
00096 particle = 0;
00097 currentCouple = 0;
00098 wokvi = new G4WentzelOKandVIxSection();
00099
00100 currentMaterialIndex = 0;
00101
00102 cosTetMinNuc = 1.0;
00103 cosTetMaxNuc = -1.0;
00104 elecRatio = 0.0;
00105 mass = proton_mass_c2;
00106 }
00107
00108
00109
00110 G4eCoulombScatteringModel::~G4eCoulombScatteringModel()
00111 {
00112 delete wokvi;
00113 }
00114
00115
00116
00117 void G4eCoulombScatteringModel::Initialise(const G4ParticleDefinition* p,
00118 const G4DataVector& cuts)
00119 {
00120 SetupParticle(p);
00121 currentCouple = 0;
00122 cosThetaMin = cos(PolarAngleLimit());
00123 wokvi->Initialise(p, cosThetaMin);
00124
00125
00126
00127
00128
00129
00130 pCuts = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3);
00131
00132
00133
00134
00135
00136
00137 if(!isInitialised) {
00138 isInitialised = true;
00139 fParticleChange = GetParticleChangeForGamma();
00140 }
00141 if(mass < GeV) {
00142 InitialiseElementSelectors(p,cuts);
00143 }
00144 }
00145
00146
00147
00148 G4double G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(
00149 const G4ParticleDefinition* p,
00150 G4double kinEnergy,
00151 G4double Z, G4double,
00152 G4double cutEnergy, G4double)
00153 {
00154
00155
00156 G4double cross = 0.0;
00157 if(p != particle) { SetupParticle(p); }
00158
00159
00160 if(kinEnergy <= 0.0) { return cross; }
00161 DefineMaterial(CurrentCouple());
00162 cosTetMinNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
00163 if(cosThetaMax < cosTetMinNuc) {
00164 G4int iz = G4int(Z);
00165 cosTetMinNuc = wokvi->SetupTarget(iz, cutEnergy);
00166 cosTetMaxNuc = cosThetaMax;
00167 if(iz == 1 && cosTetMaxNuc < 0.0 && particle == theProton) {
00168 cosTetMaxNuc = 0.0;
00169 }
00170 cross = wokvi->ComputeNuclearCrossSection(cosTetMinNuc, cosTetMaxNuc);
00171 elecRatio = wokvi->ComputeElectronCrossSection(cosTetMinNuc, cosThetaMax);
00172 cross += elecRatio;
00173 if(cross > 0.0) { elecRatio /= cross; }
00174 }
00175
00176
00177
00178
00179
00180
00181
00182 return cross;
00183 }
00184
00185
00186
00187 void G4eCoulombScatteringModel::SampleSecondaries(
00188 std::vector<G4DynamicParticle*>* fvect,
00189 const G4MaterialCutsCouple* couple,
00190 const G4DynamicParticle* dp,
00191 G4double cutEnergy,
00192 G4double)
00193 {
00194 G4double kinEnergy = dp->GetKineticEnergy();
00195
00196
00197
00198 if(kinEnergy < lowEnergyThreshold) {
00199 fParticleChange->SetProposedKineticEnergy(0.0);
00200 fParticleChange->ProposeLocalEnergyDeposit(kinEnergy);
00201 fParticleChange->ProposeNonIonizingEnergyDeposit(kinEnergy);
00202 return;
00203 }
00204 SetupParticle(dp->GetDefinition());
00205 DefineMaterial(couple);
00206
00207
00208
00209
00210
00211
00212 const G4Element* currentElement =
00213 SelectRandomAtom(couple,particle,kinEnergy,cutEnergy,kinEnergy);
00214
00215 G4double Z = currentElement->GetZ();
00216
00217 if(ComputeCrossSectionPerAtom(particle,kinEnergy, Z,
00218 kinEnergy, cutEnergy, kinEnergy) == 0.0)
00219 { return; }
00220
00221 G4int iz = G4int(Z);
00222 G4int ia = SelectIsotopeNumber(currentElement);
00223 G4double targetMass = G4NucleiProperties::GetNuclearMass(ia, iz);
00224 wokvi->SetTargetMass(targetMass);
00225
00226 G4ThreeVector newDirection =
00227 wokvi->SampleSingleScattering(cosTetMinNuc, cosThetaMax, elecRatio);
00228 G4double cost = newDirection.z();
00229
00230 G4ThreeVector direction = dp->GetMomentumDirection();
00231 newDirection.rotateUz(direction);
00232
00233 fParticleChange->ProposeMomentumDirection(newDirection);
00234
00235
00236
00237 G4double mom2 = wokvi->GetMomentumSquare();
00238 G4double trec = mom2*(1.0 - cost)/(targetMass + (mass + kinEnergy)*(1.0 - cost));
00239 G4double finalT = kinEnergy - trec;
00240
00241 if(finalT <= lowEnergyThreshold) {
00242 trec = kinEnergy;
00243 finalT = 0.0;
00244 }
00245
00246 fParticleChange->SetProposedKineticEnergy(finalT);
00247 G4double tcut = recoilThreshold;
00248 if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); }
00249
00250 if(trec > tcut) {
00251 G4ParticleDefinition* ion = theParticleTable->GetIon(iz, ia, 0.0);
00252 G4ThreeVector dir = (direction*sqrt(mom2) -
00253 newDirection*sqrt(finalT*(2*mass + finalT))).unit();
00254 G4DynamicParticle* newdp = new G4DynamicParticle(ion, dir, trec);
00255 fvect->push_back(newdp);
00256 } else {
00257 fParticleChange->ProposeLocalEnergyDeposit(trec);
00258 fParticleChange->ProposeNonIonizingEnergyDeposit(trec);
00259 }
00260
00261 return;
00262 }
00263
00264
00265
00266