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 #include "G4EmSaturation.hh"
00047 #include "G4PhysicalConstants.hh"
00048 #include "G4SystemOfUnits.hh"
00049 #include "G4LossTableManager.hh"
00050 #include "G4NistManager.hh"
00051 #include "G4Material.hh"
00052 #include "G4MaterialCutsCouple.hh"
00053 #include "G4Electron.hh"
00054 #include "G4Proton.hh"
00055
00056
00057
00058 G4EmSaturation::G4EmSaturation()
00059 {
00060 verbose = 1;
00061 manager = 0;
00062
00063 curMaterial = 0;
00064 curBirks = 0.0;
00065 curRatio = 1.0;
00066 curChargeSq = 1.0;
00067 nMaterials = 0;
00068
00069 electron = 0;
00070 proton = 0;
00071 nist = G4NistManager::Instance();
00072
00073 Initialise();
00074 }
00075
00076
00077
00078 G4EmSaturation::~G4EmSaturation()
00079 {}
00080
00081
00082
00083 G4double G4EmSaturation::VisibleEnergyDeposition(
00084 const G4ParticleDefinition* p,
00085 const G4MaterialCutsCouple* couple,
00086 G4double length,
00087 G4double edep,
00088 G4double niel)
00089 {
00090 if(edep <= 0.0) { return 0.0; }
00091
00092 G4double evis = edep;
00093 G4double bfactor = FindBirksCoefficient(couple->GetMaterial());
00094
00095 if(bfactor > 0.0) {
00096
00097 G4int pdgCode = p->GetPDGEncoding();
00098
00099 if(22 == pdgCode) {
00100 evis /= (1.0 + bfactor*edep/manager->GetRange(electron,edep,couple));
00101
00102
00103 } else {
00104
00105
00106 G4double nloss = niel;
00107 if(nloss < 0.0) nloss = 0.0;
00108 G4double eloss = edep - nloss;
00109
00110
00111 if(2112 == pdgCode || eloss < 0.0 || length <= 0.0) {
00112 nloss = edep;
00113 eloss = 0.0;
00114 }
00115
00116
00117 if(eloss > 0.0) { eloss /= (1.0 + bfactor*eloss/length); }
00118
00119
00120 if(nloss > 0.0) {
00121 if(!proton) { proton = G4Proton::Proton(); }
00122 G4double escaled = nloss*curRatio;
00123 G4double range = manager->GetRange(proton,escaled,couple)/curChargeSq;
00124 nloss /= (1.0 + bfactor*nloss/range);
00125 }
00126
00127 evis = eloss + nloss;
00128 }
00129 }
00130
00131 return evis;
00132 }
00133
00134
00135
00136 G4double G4EmSaturation::FindG4BirksCoefficient(const G4Material* mat)
00137 {
00138 G4String name = mat->GetName();
00139
00140
00141 for(G4int j=0; j<nG4Birks; ++j) {
00142 if(name == g4MatNames[j]) {
00143 if(verbose > 0)
00144 G4cout << "### G4EmSaturation::FindG4BirksCoefficient for "
00145 << name << " is " << g4MatData[j]*MeV/mm << " mm/MeV "
00146 << G4endl;
00147 return g4MatData[j];
00148 }
00149 }
00150 return FindBirksCoefficient(mat);
00151 }
00152
00153
00154
00155 G4double G4EmSaturation::FindBirksCoefficient(const G4Material* mat)
00156 {
00157
00158 if(!manager) {
00159 manager = G4LossTableManager::Instance();
00160 electron= G4Electron::Electron();
00161 }
00162
00163 if(mat == curMaterial) { return curBirks; }
00164
00165 curMaterial = mat;
00166 curBirks = 0.0;
00167 curRatio = 1.0;
00168 curChargeSq = 1.0;
00169
00170
00171 for(G4int i=0; i<nMaterials; ++i) {
00172 if(mat == matPointers[i]) {
00173 curBirks = mat->GetIonisation()->GetBirksConstant();
00174 curRatio = massFactors[i];
00175 curChargeSq = effCharges[i];
00176 return curBirks;
00177 }
00178 }
00179
00180 G4String name = mat->GetName();
00181 curBirks = mat->GetIonisation()->GetBirksConstant();
00182
00183
00184
00185 if(curBirks == 0.0) {
00186 for(G4int j=0; j<nG4Birks; ++j) {
00187 if(name == g4MatNames[j]) {
00188 mat->GetIonisation()->SetBirksConstant(g4MatData[j]);
00189 curBirks = g4MatData[j];
00190 break;
00191 }
00192 }
00193 }
00194
00195 if(curBirks == 0.0 && verbose > 0) {
00196 G4cout << "### G4EmSaturation::FindBirksCoefficient fails "
00197 " for material " << name << G4endl;
00198 }
00199
00200
00201 curRatio = 0.0;
00202 curChargeSq = 0.0;
00203 G4double norm = 0.0;
00204 const G4ElementVector* theElementVector = mat->GetElementVector();
00205 const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
00206 size_t nelm = mat->GetNumberOfElements();
00207 for (size_t i=0; i<nelm; ++i) {
00208 const G4Element* elm = (*theElementVector)[i];
00209 G4double Z = elm->GetZ();
00210 G4double w = Z*Z*theAtomNumDensityVector[i];
00211 curRatio += w/nist->GetAtomicMassAmu(G4int(Z));
00212 curChargeSq = Z*Z*w;
00213 norm += w;
00214 }
00215 curRatio *= proton_mass_c2/norm;
00216 curChargeSq /= norm;
00217
00218
00219 matPointers.push_back(mat);
00220 matNames.push_back(name);
00221 massFactors.push_back(curRatio);
00222 effCharges.push_back(curChargeSq);
00223 nMaterials++;
00224 if(curBirks > 0.0 && verbose > 0) {
00225 G4cout << "### G4EmSaturation::FindBirksCoefficient Birks coefficient for "
00226 << name << " " << curBirks*MeV/mm << " mm/MeV" << G4endl;
00227 }
00228 return curBirks;
00229 }
00230
00231
00232
00233 void G4EmSaturation::DumpBirksCoefficients()
00234 {
00235 if(nMaterials > 0) {
00236 G4cout << "### Birks coeffitients used in run time" << G4endl;
00237 for(G4int i=0; i<nMaterials; ++i) {
00238 G4double br = matPointers[i]->GetIonisation()->GetBirksConstant();
00239 G4cout << " " << matNames[i] << " "
00240 << br*MeV/mm << " mm/MeV" << " "
00241 << br*matPointers[i]->GetDensity()*MeV*cm2/g
00242 << " g/cm^2/MeV"
00243 << G4endl;
00244 }
00245 }
00246 }
00247
00248
00249
00250 void G4EmSaturation::DumpG4BirksCoefficients()
00251 {
00252 if(nG4Birks > 0) {
00253 G4cout << "### Birks coeffitients for Geant4 materials" << G4endl;
00254 for(G4int i=0; i<nG4Birks; ++i) {
00255 G4cout << " " << g4MatNames[i] << " "
00256 << g4MatData[i]*MeV/mm << " mm/MeV" << G4endl;
00257 }
00258 }
00259 }
00260
00261
00262
00263 void G4EmSaturation::Initialise()
00264 {
00265
00266
00267 g4MatNames.push_back("G4_POLYSTYRENE");
00268 g4MatData.push_back(0.07943*mm/MeV);
00269
00270
00271
00272 g4MatNames.push_back("G4_BGO");
00273 g4MatData.push_back(0.008415*mm/MeV);
00274
00275
00276
00277
00278
00279
00280 g4MatNames.push_back("G4_lAr");
00281 g4MatData.push_back(0.1576*mm/MeV);
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294 nG4Birks = g4MatData.size();
00295 }
00296
00297