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