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 #include "G4PixeCrossSectionHandler.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4IInterpolator.hh"
00043 #include "G4LogLogInterpolator.hh"
00044 #include "G4IDataSet.hh"
00045 #include "G4DataSet.hh"
00046 #include "G4CompositeDataSet.hh"
00047 #include "G4PixeShellDataSet.hh"
00048 #include "G4ProductionCutsTable.hh"
00049 #include "G4Material.hh"
00050 #include "G4Element.hh"
00051 #include "Randomize.hh"
00052 #include "G4SystemOfUnits.hh"
00053 #include "G4ParticleDefinition.hh"
00054
00055 #include <map>
00056 #include <vector>
00057 #include <fstream>
00058 #include <sstream>
00059
00060
00061 G4PixeCrossSectionHandler::G4PixeCrossSectionHandler()
00062 {
00063 crossSections = 0;
00064 interpolation = 0;
00065
00066 Initialise(0,"","","",1.*keV,0.1*GeV,200,MeV,barn,6,92);
00067 ActiveElements();
00068 }
00069
00070
00071 G4PixeCrossSectionHandler::G4PixeCrossSectionHandler(G4IInterpolator* algorithm,
00072 const G4String& modelK,
00073 const G4String& modelL,
00074 const G4String& modelM,
00075 G4double minE,
00076 G4double maxE,
00077 G4int bins,
00078 G4double unitE,
00079 G4double unitData,
00080 G4int minZ,
00081 G4int maxZ)
00082 : interpolation(algorithm), eMin(minE), eMax(maxE), nBins(bins),
00083 unit1(unitE), unit2(unitData), zMin(minZ), zMax(maxZ)
00084 {
00085 crossSections = 0;
00086
00087 crossModel.push_back(modelK);
00088 crossModel.push_back(modelL);
00089 crossModel.push_back(modelM);
00090
00091
00092
00093
00094
00095 ActiveElements();
00096 }
00097
00098 G4PixeCrossSectionHandler::~G4PixeCrossSectionHandler()
00099 {
00100 delete interpolation;
00101 interpolation = 0;
00102 std::map<G4int,G4IDataSet*,std::less<G4int> >::iterator pos;
00103
00104 for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
00105 {
00106
00107
00108
00109
00110 G4IDataSet* dataSet = (*pos).second;
00111 delete dataSet;
00112 }
00113
00114 if (crossSections != 0)
00115 {
00116 size_t n = crossSections->size();
00117 for (size_t i=0; i<n; i++)
00118 {
00119 delete (*crossSections)[i];
00120 }
00121 delete crossSections;
00122 crossSections = 0;
00123 }
00124 }
00125
00126 void G4PixeCrossSectionHandler::Initialise(G4IInterpolator* algorithm,
00127 const G4String& modelK,
00128 const G4String& modelL,
00129 const G4String& modelM,
00130 G4double minE, G4double maxE,
00131 G4int numberOfBins,
00132 G4double unitE, G4double unitData,
00133 G4int minZ, G4int maxZ)
00134 {
00135 if (algorithm != 0)
00136 {
00137 delete interpolation;
00138 interpolation = algorithm;
00139 }
00140 else
00141 {
00142 interpolation = CreateInterpolation();
00143 }
00144
00145 eMin = minE;
00146 eMax = maxE;
00147 nBins = numberOfBins;
00148 unit1 = unitE;
00149 unit2 = unitData;
00150 zMin = minZ;
00151 zMax = maxZ;
00152
00153 crossModel.push_back(modelK);
00154 crossModel.push_back(modelL);
00155 crossModel.push_back(modelM);
00156
00157 }
00158
00159 void G4PixeCrossSectionHandler::PrintData() const
00160 {
00161 std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
00162
00163 for (pos = dataMap.begin(); pos != dataMap.end(); pos++)
00164 {
00165
00166
00167
00168
00169
00170 G4int z = (*pos).first;
00171 G4IDataSet* dataSet = (*pos).second;
00172 G4cout << "---- Data set for Z = "
00173 << z
00174 << G4endl;
00175 dataSet->PrintData();
00176 G4cout << "--------------------------------------------------" << G4endl;
00177 }
00178 }
00179
00180 void G4PixeCrossSectionHandler::LoadShellData(const G4String& fileName)
00181 {
00182 size_t nZ = activeZ.size();
00183 for (size_t i=0; i<nZ; i++)
00184 {
00185 G4int Z = (G4int) activeZ[i];
00186 G4IInterpolator* algo = interpolation->Clone();
00187 G4IDataSet* dataSet = new G4PixeShellDataSet(Z, algo,crossModel[0],crossModel[1],crossModel[2]);
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198 dataSet->LoadData(fileName);
00199 dataMap[Z] = dataSet;
00200 }
00201
00202
00203 if (! crossSections)
00204 {
00205 BuildForMaterials();
00206 }
00207
00208 }
00209
00210 void G4PixeCrossSectionHandler::Clear()
00211 {
00212
00213 std::map<G4int,G4IDataSet*,std::less<G4int> >::iterator pos;
00214
00215 if(! dataMap.empty())
00216 {
00217 for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
00218 {
00219
00220
00221
00222
00223 G4IDataSet* dataSet = (*pos).second;
00224 delete dataSet;
00225 dataSet = 0;
00226 G4int i = (*pos).first;
00227 dataMap[i] = 0;
00228 }
00229 dataMap.clear();
00230 }
00231
00232 activeZ.clear();
00233 ActiveElements();
00234 }
00235
00236 G4double G4PixeCrossSectionHandler::FindValue(G4int Z, G4double energy) const
00237 {
00238 G4double value = 0.;
00239
00240 std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
00241 pos = dataMap.find(Z);
00242 if (pos!= dataMap.end())
00243 {
00244
00245
00246
00247
00248 G4IDataSet* dataSet = (*pos).second;
00249 value = dataSet->FindValue(energy);
00250 }
00251 else
00252 {
00253 G4cout << "WARNING: G4PixeCrossSectionHandler::FindValue(Z,e) did not find Z = "
00254 << Z << G4endl;
00255 }
00256 return value;
00257 }
00258
00259 G4double G4PixeCrossSectionHandler::FindValue(G4int Z, G4double energy,
00260 G4int shellIndex) const
00261 {
00262 G4double value = 0.;
00263
00264 std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
00265 pos = dataMap.find(Z);
00266 if (pos!= dataMap.end())
00267 {
00268
00269
00270
00271
00272 G4IDataSet* dataSet = (*pos).second;
00273 if (shellIndex >= 0)
00274 {
00275 G4int nComponents = dataSet->NumberOfComponents();
00276 if(shellIndex < nComponents)
00277
00278 value = dataSet->GetComponent(shellIndex)->FindValue(energy);
00279 else
00280 {
00281 G4cout << "WARNING: G4PixeCrossSectionHandler::FindValue(Z,e,shell) did not find"
00282 << " shellIndex= " << shellIndex
00283 << " for Z= "
00284 << Z << G4endl;
00285 }
00286 } else {
00287 value = dataSet->FindValue(energy);
00288 }
00289 }
00290 else
00291 {
00292 G4cout << "WARNING: G4PixeCrossSectionHandler::FindValue did not find Z = "
00293 << Z << G4endl;
00294 }
00295 return value;
00296 }
00297
00298
00299 G4double G4PixeCrossSectionHandler::ValueForMaterial(const G4Material* material,
00300 G4double energy) const
00301 {
00302 G4double value = 0.;
00303
00304 const G4ElementVector* elementVector = material->GetElementVector();
00305 const G4double* nAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
00306 G4int nElements = material->GetNumberOfElements();
00307
00308 for (G4int i=0 ; i<nElements ; i++)
00309 {
00310 G4int Z = (G4int) (*elementVector)[i]->GetZ();
00311 G4double elementValue = FindValue(Z,energy);
00312 G4double nAtomsVol = nAtomsPerVolume[i];
00313 value += nAtomsVol * elementValue;
00314 }
00315
00316 return value;
00317 }
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406 void G4PixeCrossSectionHandler::BuildForMaterials()
00407 {
00408
00409
00410
00411 G4DataVector energyVector;
00412 G4double dBin = std::log10(eMax/eMin) / nBins;
00413
00414 for (G4int i=0; i<nBins+1; i++)
00415 {
00416 energyVector.push_back(std::pow(10., std::log10(eMin)+i*dBin));
00417 }
00418
00419 if (crossSections != 0)
00420 {
00421 std::vector<G4IDataSet*>::iterator mat;
00422 if (! crossSections->empty())
00423 {
00424 for (mat = crossSections->begin(); mat!= crossSections->end(); ++mat)
00425 {
00426 G4IDataSet* set = *mat;
00427 delete set;
00428 set = 0;
00429 }
00430 crossSections->clear();
00431 delete crossSections;
00432 crossSections = 0;
00433 }
00434 }
00435
00436 crossSections = BuildCrossSectionsForMaterials(energyVector);
00437
00438 if (crossSections == 0)
00439 G4Exception("G4PixeCrossSectionHandler::BuildForMaterials",
00440 "pii00000210",
00441 FatalException,
00442 ", crossSections = 0");
00443
00444 return;
00445 }
00446
00447
00448 G4int G4PixeCrossSectionHandler::SelectRandomAtom(const G4Material* material,
00449 G4double e) const
00450 {
00451
00452
00453
00454 G4int nElements = material->GetNumberOfElements();
00455
00456
00457 if (nElements == 1)
00458 {
00459 G4int Z = (G4int) material->GetZ();
00460 return Z;
00461 }
00462
00463
00464
00465 const G4ElementVector* elementVector = material->GetElementVector();
00466 size_t materialIndex = material->GetIndex();
00467
00468 G4IDataSet* materialSet = (*crossSections)[materialIndex];
00469 G4double materialCrossSection0 = 0.0;
00470 G4DataVector cross;
00471 cross.clear();
00472 for ( G4int i=0; i < nElements; i++ )
00473 {
00474 G4double cr = materialSet->GetComponent(i)->FindValue(e);
00475 materialCrossSection0 += cr;
00476 cross.push_back(materialCrossSection0);
00477 }
00478
00479 G4double random = G4UniformRand() * materialCrossSection0;
00480
00481 for (G4int k=0 ; k < nElements ; k++ )
00482 {
00483 if (random <= cross[k]) return (G4int) (*elementVector)[k]->GetZ();
00484 }
00485
00486 return 0;
00487 }
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538 G4int G4PixeCrossSectionHandler::SelectRandomShell(G4int Z, G4double e) const
00539 {
00540
00541
00542
00543
00544
00545
00546 G4int shell = 0;
00547
00548 G4double totCrossSection = FindValue(Z,e);
00549 G4double random = G4UniformRand() * totCrossSection;
00550 G4double partialSum = 0.;
00551
00552 G4IDataSet* dataSet = 0;
00553 std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
00554 pos = dataMap.find(Z);
00555
00556
00557
00558
00559 if (pos != dataMap.end()) dataSet = (*pos).second;
00560
00561 size_t nShells = dataSet->NumberOfComponents();
00562 for (size_t i=0; i<nShells; i++)
00563 {
00564 const G4IDataSet* shellDataSet = dataSet->GetComponent(i);
00565 if (shellDataSet != 0)
00566 {
00567 G4double value = shellDataSet->FindValue(e);
00568 partialSum += value;
00569 if (random <= partialSum) return i;
00570 }
00571 }
00572
00573 return shell;
00574 }
00575
00576 void G4PixeCrossSectionHandler::ActiveElements()
00577 {
00578 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
00579 if (materialTable == 0)
00580 G4Exception("G4PixeCrossSectionHandler::ActiveElements",
00581 "pii00000220",
00582 FatalException,
00583 "no MaterialTable found");
00584
00585 G4int nMaterials = G4Material::GetNumberOfMaterials();
00586
00587 for (G4int mat=0; mat<nMaterials; mat++)
00588 {
00589 const G4Material* material= (*materialTable)[mat];
00590 const G4ElementVector* elementVector = material->GetElementVector();
00591 const G4int nElements = material->GetNumberOfElements();
00592
00593 for (G4int iEl=0; iEl<nElements; iEl++)
00594 {
00595 G4Element* element = (*elementVector)[iEl];
00596 G4double Z = element->GetZ();
00597 if (!(activeZ.contains(Z)) && Z >= zMin && Z <= zMax)
00598 {
00599 activeZ.push_back(Z);
00600 }
00601 }
00602 }
00603 }
00604
00605 G4IInterpolator* G4PixeCrossSectionHandler::CreateInterpolation()
00606 {
00607 G4IInterpolator* algorithm = new G4LogLogInterpolator;
00608 return algorithm;
00609 }
00610
00611 G4int G4PixeCrossSectionHandler::NumberOfComponents(G4int Z) const
00612 {
00613 G4int n = 0;
00614
00615 std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
00616 pos = dataMap.find(Z);
00617 if (pos!= dataMap.end())
00618 {
00619 G4IDataSet* dataSet = (*pos).second;
00620 n = dataSet->NumberOfComponents();
00621 }
00622 else
00623 {
00624 G4cout << "WARNING: G4PixeCrossSectionHandler::NumberOfComponents did not "
00625 << "find Z = "
00626 << Z << G4endl;
00627 }
00628 return n;
00629 }
00630
00631
00632 std::vector<G4IDataSet*>*
00633 G4PixeCrossSectionHandler::BuildCrossSectionsForMaterials(const G4DataVector& energyVector)
00634 {
00635 G4DataVector* energies;
00636 G4DataVector* data;
00637
00638 std::vector<G4IDataSet*>* matCrossSections = new std::vector<G4IDataSet*>;
00639
00640
00641
00642
00643 size_t nOfBins = energyVector.size();
00644 const G4IInterpolator* interpolationAlgo = CreateInterpolation();
00645
00646 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
00647 if (materialTable == 0)
00648 G4Exception("G4PixeCrossSectionHandler::BuildCrossSectionsForMaterials",
00649 "pii00000230",
00650 FatalException,
00651 "no MaterialTable found");
00652
00653 G4int nMaterials = G4Material::GetNumberOfMaterials();
00654
00655 for (G4int mat=0; mat<nMaterials; mat++)
00656 {
00657 const G4Material* material = (*materialTable)[mat];
00658 G4int nElements = material->GetNumberOfElements();
00659 const G4ElementVector* elementVector = material->GetElementVector();
00660 const G4double* nAtomsPerVolume = material->GetAtomicNumDensityVector();
00661
00662 G4IInterpolator* algo = interpolationAlgo->Clone();
00663
00664 G4IDataSet* setForMat = new G4CompositeDataSet(algo,1.,1.);
00665
00666 for (G4int i=0; i<nElements; i++) {
00667
00668 G4int Z = (G4int) (*elementVector)[i]->GetZ();
00669 G4double density = nAtomsPerVolume[i];
00670
00671 energies = new G4DataVector;
00672 data = new G4DataVector;
00673
00674
00675 for (size_t bin=0; bin<nOfBins; bin++)
00676 {
00677 G4double e = energyVector[bin];
00678 energies->push_back(e);
00679 G4double cross = 0.;
00680 if (Z >= zMin && Z <= zMax) cross = density*FindValue(Z,e);
00681 data->push_back(cross);
00682 }
00683
00684 G4IInterpolator* algo1 = interpolationAlgo->Clone();
00685 G4IDataSet* elSet = new G4DataSet(i,energies,data,algo1,1.,1.);
00686 setForMat->AddComponent(elSet);
00687 }
00688
00689 matCrossSections->push_back(setForMat);
00690 }
00691 return matCrossSections;
00692 }
00693
00694
00695 G4double G4PixeCrossSectionHandler::MicroscopicCrossSection(const G4ParticleDefinition* particleDef,
00696 G4double kineticEnergy,
00697 G4double Z,
00698 G4double deltaCut) const
00699 {
00700
00701
00702
00703
00704
00705 G4double cross = 0.;
00706
00707
00708 G4double particleMass = particleDef->GetPDGMass();
00709 G4double energy = kineticEnergy + particleMass;
00710
00711
00712 G4double gamma = energy / particleMass;
00713 G4double beta2 = 1. - 1. / (gamma * gamma);
00714 G4double var = electron_mass_c2 / particleMass;
00715 G4double tMax = 2. * electron_mass_c2 * (gamma*gamma - 1.) / (1. + 2.*gamma*var + var*var);
00716
00717
00718
00719 if ( tMax > deltaCut )
00720 {
00721 var = deltaCut / tMax;
00722 cross = (1. - var * (1. - beta2 * std::log(var))) / deltaCut;
00723
00724 G4double spin = particleDef->GetPDGSpin() ;
00725
00726
00727 if (spin == 0.5)
00728 {
00729 cross += 0.5 * (tMax - deltaCut) / (energy*energy);
00730 }
00731
00732 else if (spin > 0.9 )
00733 {
00734 cross += -std::log(var) / (3.*deltaCut) + (tMax-deltaCut) *
00735 ((5.+1./var)*0.25 /(energy*energy) - beta2 / (tMax*deltaCut))/3.;
00736 }
00737 cross *= twopi_mc2_rcl2 * Z / beta2 ;
00738 }
00739
00740
00741
00742
00743 return cross;
00744 }
00745