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 #include "G4PenelopeCrossSection.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "G4PhysicsTable.hh"
00040 #include "G4PhysicsFreeVector.hh"
00041 #include "G4DataVector.hh"
00042
00043
00044 G4PenelopeCrossSection::G4PenelopeCrossSection(size_t nPointsE,size_t nShells) :
00045 numberOfEnergyPoints(nPointsE),numberOfShells(nShells),softCrossSections(0),
00046 hardCrossSections(0),shellCrossSections(0),shellNormalizedCrossSections(0)
00047 {
00048
00049 if (!numberOfEnergyPoints)
00050 {
00051 G4ExceptionDescription ed;
00052 ed << "G4PenelopeCrossSection: invalid number of energy points " << G4endl;
00053 G4Exception("G4PenelopeCrossSection::G4PenelopeCrossSection()",
00054 "em2017",FatalException,ed);
00055 }
00056
00057 isNormalized = false;
00058
00059
00060 softCrossSections = new G4PhysicsTable();
00061
00062
00063
00064
00065
00066
00067 for (size_t i=0;i<3;i++)
00068 softCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
00069
00070
00071 hardCrossSections = new G4PhysicsTable();
00072
00073
00074
00075
00076
00077
00078 for (size_t i=0;i<3;i++)
00079 hardCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
00080
00081
00082 if (numberOfShells)
00083 {
00084 shellCrossSections = new G4PhysicsTable();
00085 shellNormalizedCrossSections = new G4PhysicsTable();
00086
00087
00088 for (size_t i=0;i<numberOfShells;i++)
00089 {
00090 shellCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
00091 shellNormalizedCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
00092 }
00093 }
00094 }
00095
00096
00097 G4PenelopeCrossSection::~G4PenelopeCrossSection()
00098 {
00099
00100 if (shellCrossSections)
00101 {
00102 shellCrossSections->clearAndDestroy();
00103 delete shellCrossSections;
00104 }
00105 if (shellNormalizedCrossSections)
00106 {
00107 shellNormalizedCrossSections->clearAndDestroy();
00108 delete shellNormalizedCrossSections;
00109 }
00110 if (softCrossSections)
00111 {
00112 softCrossSections->clearAndDestroy();
00113 delete softCrossSections;
00114 }
00115 if (hardCrossSections)
00116 {
00117 hardCrossSections->clearAndDestroy();
00118 delete hardCrossSections;
00119 }
00120 }
00121
00122
00123 void G4PenelopeCrossSection::AddCrossSectionPoint(size_t binNumber,G4double energy,
00124 G4double XH0,
00125 G4double XH1, G4double XH2,
00126 G4double XS0, G4double XS1,
00127 G4double XS2)
00128 {
00129 if (!softCrossSections || !hardCrossSections)
00130 {
00131 G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
00132 G4endl;
00133 G4cout << "Trying to fill un-initialized tables" << G4endl;
00134 return;
00135 }
00136
00137
00138 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*softCrossSections)[0];
00139
00140 if (binNumber >= numberOfEnergyPoints)
00141 {
00142 G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
00143 G4endl;
00144 G4cout << "Trying to register more points than originally declared" << G4endl;
00145 return;
00146 }
00147 G4double logEne = std::log(energy);
00148
00149
00150 G4double val = std::log(std::max(XS0,1e-42*cm2));
00151 theVector->PutValue(binNumber,logEne,val);
00152
00153
00154 theVector = (G4PhysicsFreeVector*) (*softCrossSections)[1];
00155 val = std::log(std::max(XS1,1e-42*eV*cm2));
00156 theVector->PutValue(binNumber,logEne,val);
00157
00158
00159 theVector = (G4PhysicsFreeVector*) (*softCrossSections)[2];
00160 val = std::log(std::max(XS2,1e-42*eV*eV*cm2));
00161 theVector->PutValue(binNumber,logEne,val);
00162
00163
00164 theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
00165 val = std::log(std::max(XH0,1e-42*cm2));
00166 theVector->PutValue(binNumber,logEne,val);
00167
00168
00169 theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[1];
00170 val = std::log(std::max(XH1,1e-42*eV*cm2));
00171 theVector->PutValue(binNumber,logEne,val);
00172
00173
00174 theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[2];
00175 val = std::log(std::max(XH2,1e-42*eV*eV*cm2));
00176 theVector->PutValue(binNumber,logEne,val);
00177
00178 return;
00179 }
00180
00181
00182
00183 void G4PenelopeCrossSection::AddShellCrossSectionPoint(size_t binNumber,
00184 size_t shellID,
00185 G4double energy,
00186 G4double xs)
00187 {
00188 if (!shellCrossSections)
00189 {
00190 G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
00191 G4endl;
00192 G4cout << "Trying to fill un-initialized table" << G4endl;
00193 return;
00194 }
00195
00196 if (shellID >= numberOfShells)
00197 {
00198 G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
00199 G4endl;
00200 G4cout << "Trying to fill shell #" << shellID << " while the maximum is "
00201 << numberOfShells-1 << G4endl;
00202 return;
00203 }
00204
00205
00206 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
00207
00208 if (binNumber >= numberOfEnergyPoints)
00209 {
00210 G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
00211 G4endl;
00212 G4cout << "Trying to register more points than originally declared" << G4endl;
00213 return;
00214 }
00215 G4double logEne = std::log(energy);
00216 G4double val = std::log(std::max(xs,1e-42*cm2));
00217 theVector->PutValue(binNumber,logEne,val);
00218
00219 return;
00220 }
00221
00222
00223
00224 G4double G4PenelopeCrossSection::GetTotalCrossSection(G4double energy)
00225 {
00226 G4double result = 0;
00227
00228 if (!softCrossSections || !hardCrossSections)
00229 {
00230 G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
00231 G4endl;
00232 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
00233 return result;
00234 }
00235
00236
00237 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*softCrossSections)[0];
00238 if (theVector->GetVectorLength() < numberOfEnergyPoints)
00239 {
00240 G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
00241 G4endl;
00242 G4cout << "Soft cross section table looks not filled" << G4endl;
00243 return result;
00244 }
00245 G4double logene = std::log(energy);
00246 G4double logXS = theVector->Value(logene);
00247 G4double softXS = std::exp(logXS);
00248
00249
00250 theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
00251 if (theVector->GetVectorLength() < numberOfEnergyPoints)
00252 {
00253 G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
00254 G4endl;
00255 G4cout << "Hard cross section table looks not filled" << G4endl;
00256 return result;
00257 }
00258 logXS = theVector->Value(logene);
00259 G4double hardXS = std::exp(logXS);
00260
00261 result = hardXS + softXS;
00262 return result;
00263
00264 }
00265
00266
00267
00268 G4double G4PenelopeCrossSection::GetHardCrossSection(G4double energy)
00269 {
00270 G4double result = 0;
00271
00272 if (!hardCrossSections)
00273 {
00274 G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
00275 G4endl;
00276 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
00277 return result;
00278 }
00279
00280 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
00281 if (theVector->GetVectorLength() < numberOfEnergyPoints)
00282 {
00283 G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
00284 G4endl;
00285 G4cout << "Hard cross section table looks not filled" << G4endl;
00286 return result;
00287 }
00288 G4double logene = std::log(energy);
00289 G4double logXS = theVector->Value(logene);
00290 result = std::exp(logXS);
00291
00292 return result;
00293 }
00294
00295
00296
00297
00298 G4double G4PenelopeCrossSection::GetSoftStoppingPower(G4double energy)
00299 {
00300 G4double result = 0;
00301
00302 if (!softCrossSections)
00303 {
00304 G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
00305 G4endl;
00306 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
00307 return result;
00308 }
00309
00310 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*softCrossSections)[1];
00311 if (theVector->GetVectorLength() < numberOfEnergyPoints)
00312 {
00313 G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
00314 G4endl;
00315 G4cout << "Soft cross section table looks not filled" << G4endl;
00316 return result;
00317 }
00318 G4double logene = std::log(energy);
00319 G4double logXS = theVector->Value(logene);
00320 result = std::exp(logXS);
00321
00322 return result;
00323 }
00324
00325
00326
00327 G4double G4PenelopeCrossSection::GetShellCrossSection(size_t shellID,G4double energy)
00328 {
00329 G4double result = 0;
00330 if (!shellCrossSections)
00331 {
00332 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
00333 G4endl;
00334 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
00335 return result;
00336 }
00337 if (shellID >= numberOfShells)
00338 {
00339 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
00340 G4endl;
00341 G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
00342 << numberOfShells-1 << G4endl;
00343 return result;
00344 }
00345
00346 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
00347
00348 if (theVector->GetVectorLength() < numberOfEnergyPoints)
00349 {
00350 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
00351 G4endl;
00352 G4cout << "Shell cross section table looks not filled" << G4endl;
00353 return result;
00354 }
00355 G4double logene = std::log(energy);
00356 G4double logXS = theVector->Value(logene);
00357 result = std::exp(logXS);
00358
00359 return result;
00360 }
00361
00362
00363 G4double G4PenelopeCrossSection::GetNormalizedShellCrossSection(size_t shellID,G4double energy)
00364 {
00365 G4double result = 0;
00366 if (!shellNormalizedCrossSections)
00367 {
00368 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
00369 G4endl;
00370 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
00371 return result;
00372 }
00373
00374 if (!isNormalized)
00375 NormalizeShellCrossSections();
00376
00377 if (shellID >= numberOfShells)
00378 {
00379 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
00380 G4endl;
00381 G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
00382 << numberOfShells-1 << G4endl;
00383 return result;
00384 }
00385
00386 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellNormalizedCrossSections)[shellID];
00387
00388 if (theVector->GetVectorLength() < numberOfEnergyPoints)
00389 {
00390 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
00391 G4endl;
00392 G4cout << "Shell cross section table looks not filled" << G4endl;
00393 return result;
00394 }
00395 G4double logene = std::log(energy);
00396 G4double logXS = theVector->Value(logene);
00397 result = std::exp(logXS);
00398
00399 return result;
00400 }
00401
00402
00403
00404
00405
00406 void G4PenelopeCrossSection::NormalizeShellCrossSections()
00407 {
00408 if (isNormalized)
00409 {
00410 G4cout << "G4PenelopeCrossSection::NormalizeShellCrossSections()" << G4endl;
00411 G4cout << "already invoked. Ignore it" << G4endl;
00412 return;
00413 }
00414
00415 if (!shellNormalizedCrossSections)
00416 {
00417 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
00418 G4endl;
00419 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
00420 return;
00421 }
00422
00423 for (size_t i=0;i<numberOfEnergyPoints;i++)
00424 {
00425
00426
00427
00428
00429 G4double normFactor = 0.;
00430 for (size_t shellID=0;shellID<numberOfShells;shellID++)
00431 {
00432 G4PhysicsFreeVector* theVec =
00433 (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
00434
00435 normFactor += std::exp((*theVec)[i]);
00436 }
00437 G4double logNormFactor = std::log(normFactor);
00438
00439 for (size_t shellID=0;shellID<numberOfShells;shellID++)
00440 {
00441 G4PhysicsFreeVector* theVec =
00442 (G4PhysicsFreeVector*) (*shellNormalizedCrossSections)[shellID];
00443 G4PhysicsFreeVector* theFullVec =
00444 (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
00445 G4double previousValue = (*theFullVec)[i];
00446 G4double logEnergy = theFullVec->GetLowEdgeEnergy(i);
00447
00448 theVec->PutValue(i,logEnergy,previousValue-logNormFactor);
00449 }
00450 }
00451
00452 isNormalized = true;
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470 return;
00471 }