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 #include <fstream>
00046
00047 #include "G4BremsstrahlungParameters.hh"
00048 #include "G4PhysicalConstants.hh"
00049 #include "G4VEMDataSet.hh"
00050 #include "G4EMDataSet.hh"
00051 #include "G4LogLogInterpolation.hh"
00052 #include "G4Material.hh"
00053
00054 G4BremsstrahlungParameters:: G4BremsstrahlungParameters(const G4String& name,
00055 size_t num, G4int minZ, G4int maxZ)
00056 : zMin(minZ),
00057 zMax(maxZ),
00058 length(num)
00059 {
00060 LoadData(name);
00061 }
00062
00063
00064 G4BremsstrahlungParameters::~G4BremsstrahlungParameters()
00065 {
00066
00067 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::iterator pos;
00068
00069 for (pos = param.begin(); pos != param.end(); ++pos)
00070 {
00071 G4VEMDataSet* dataSet = (*pos).second;
00072 delete dataSet;
00073 }
00074
00075 activeZ.clear();
00076 paramC.clear();
00077 }
00078
00079
00080 G4double G4BremsstrahlungParameters::Parameter(G4int parameterIndex,
00081 G4int Z,
00082 G4double energy) const
00083 {
00084 G4double value = 0.;
00085 G4int id = Z*length + parameterIndex;
00086 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
00087
00088 pos = param.find(id);
00089 if (pos!= param.end()) {
00090
00091 G4VEMDataSet* dataSet = (*pos).second;
00092 const G4DataVector ener = dataSet->GetEnergies(0);
00093 G4double ee = std::max(ener.front(),std::min(ener.back(),energy));
00094 value = dataSet->FindValue(ee);
00095
00096 } else {
00097 G4cout << "WARNING: G4BremsstrahlungParameters::FindValue "
00098 << "did not find ID = "
00099 << id << G4endl;
00100 }
00101
00102 return value;
00103 }
00104
00105 void G4BremsstrahlungParameters::LoadData(const G4String& name)
00106 {
00107 const G4double mConst =
00108 classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi;
00109
00110
00111
00112
00113 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
00114 if (materialTable == 0)
00115 G4Exception("G4BremsstrahlungParameters::LoadData",
00116 "em1001",FatalException,"Unable to find MaterialTable");
00117
00118 G4int nMaterials = G4Material::GetNumberOfMaterials();
00119
00120 G4double x = 1.e-9;
00121 for (G4int mmLocal=0; mmLocal<100; mmLocal++) {
00122 paramC.push_back(x);
00123 }
00124
00125 for (G4int mLocal=0; mLocal<nMaterials; mLocal++) {
00126
00127 const G4Material* material= (*materialTable)[mLocal];
00128 const G4ElementVector* elementVector = material->GetElementVector();
00129 const G4int nElements = material->GetNumberOfElements();
00130
00131 for (G4int iEl=0; iEl<nElements; iEl++) {
00132 G4Element* element = (*elementVector)[iEl];
00133 G4double Z = element->GetZ();
00134 G4int iz = (G4int)Z;
00135 if(iz < 100) {
00136 paramC[iz] = mConst*material->GetTotNbOfElectPerVolume();
00137
00138 }
00139 if (!(activeZ.contains(Z))) {
00140 activeZ.push_back(Z);
00141 }
00142 }
00143 }
00144
00145
00146
00147 char* path = getenv("G4LEDATA");
00148 if (path == 0)
00149 {
00150 G4Exception("G4BremsstrahlungParameters::LoadData",
00151 "em0006",FatalException,"G4LEDATA environment variable not set");
00152 return;
00153 }
00154
00155 G4String pathString_a(path);
00156 G4String name_a = pathString_a + name;
00157 std::ifstream file_a(name_a);
00158 std::filebuf* lsdp_a = file_a.rdbuf();
00159
00160 if (! (lsdp_a->is_open()) )
00161 {
00162 G4String stringConversion2("G4BremsstrahlungParameters::LoadData");
00163 G4String excep = stringConversion2 + name_a;
00164 G4Exception("G4BremsstrahlungParameters::LoadData",
00165 "em0003",FatalException,excep);
00166 }
00167
00168
00169
00170
00171
00172
00173
00174 G4double ener = 0.0;
00175 G4double sum = 0.0;
00176 G4int z = 0;
00177
00178 std::vector<G4DataVector*> a;
00179 a.resize(length);
00180
00181 G4DataVector e;
00182 e.clear();
00183
00184 G4bool isReady = false;
00185
00186 do {
00187 file_a >> ener >> sum;
00188
00189
00190 if (ener == (G4double)(-2)) {
00191 break;
00192
00193
00194 } else if (ener == (G4double)(-1)) {
00195
00196 ++z;
00197 G4double Z = (G4double)z;
00198
00199
00200 if (activeZ.contains(Z)) {
00201
00202 for (size_t k=0; k<length; ++k) {
00203
00204 G4int id = z*length + k;
00205 G4VDataSetAlgorithm* inter = new G4LogLogInterpolation();
00206 G4DataVector* eVector = new G4DataVector;
00207 size_t eSize = e.size();
00208 for (size_t sLocal=0; sLocal<eSize; sLocal++) {
00209 eVector->push_back(e[sLocal]);
00210 }
00211 G4VEMDataSet* set = new G4EMDataSet(id,eVector,a[k],inter,1.,1.);
00212 param[id] = set;
00213 }
00214 } else {
00215 for (size_t j=0; j<length; j++) {
00216 delete a[j];
00217 }
00218 }
00219 isReady = false;
00220
00221 } else {
00222
00223 if(!isReady) {
00224 isReady = true;
00225 e.clear();
00226 for (size_t j=0; j<length; ++j) {
00227 a[j] = new G4DataVector();
00228 }
00229 }
00230
00231 if(ener > 1000.) ener = 1000.;
00232 e.push_back(ener);
00233 a[length-1]->push_back(sum);
00234
00235 for (size_t j=0; j<length-1; j++) {
00236 G4double qRead;
00237 file_a >> qRead;
00238 a[j]->push_back(qRead);
00239 }
00240
00241 }
00242 } while (ener != (G4double)(-2));
00243
00244 file_a.close();
00245
00246 }
00247
00248
00249 G4double G4BremsstrahlungParameters::ParameterC(G4int id) const
00250 {
00251 G4int n = paramC.size();
00252 if (id < 0 || id >= n)
00253 {
00254 G4String stringConversion2(id);
00255 G4String ex = "Wrong id " + stringConversion2;
00256 G4Exception("G4BremsstrahlungParameters::ParameterC",
00257 "em1002",FatalException,ex);
00258
00259 }
00260
00261 return paramC[id];
00262 }
00263
00264
00265 void G4BremsstrahlungParameters::PrintData() const
00266 {
00267
00268 G4cout << G4endl;
00269 G4cout << "===== G4BremsstrahlungParameters =====" << G4endl;
00270 G4cout << G4endl;
00271 G4cout << "===== Parameters =====" << G4endl;
00272 G4cout << G4endl;
00273
00274 size_t nZ = activeZ.size();
00275 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
00276
00277 for (size_t j=0; j<nZ; j++) {
00278 G4int Z = (G4int)activeZ[j];
00279
00280 for (size_t i=0; i<length; i++) {
00281
00282 pos = param.find(Z*length + i);
00283 if (pos!= param.end()) {
00284
00285 G4cout << "===== Z= " << Z
00286 << " parameter[" << i << "] ====="
00287 << G4endl;
00288 G4VEMDataSet* dataSet = (*pos).second;
00289 dataSet->PrintData();
00290 }
00291 }
00292 }
00293
00294 G4cout << "==========================================" << G4endl;
00295 }
00296