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