| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@
29 //
30 // History:
31 // -----------
32 // 16 Jun 2008 MGP Created; Cross section ma
33 // Documented in:
34 // M.G. Pia et al., PIXE Sim
35 // IEEE Trans. Nucl. Sci., v
36 //
37 // -------------------------------------------
38
39 #include "G4PixeCrossSectionHandler.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4IInterpolator.hh"
42 #include "G4LogLogInterpolator.hh"
43 #include "G4IDataSet.hh"
44 #include "G4DataSet.hh"
45 #include "G4CompositeDataSet.hh"
46 #include "G4PixeShellDataSet.hh"
47 #include "G4ProductionCutsTable.hh"
48 #include "G4Material.hh"
49 #include "G4Element.hh"
50 #include "Randomize.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4ParticleDefinition.hh"
53
54 #include <map>
55 #include <vector>
56 #include <fstream>
57 #include <sstream>
58 #include <memory>
59
60
61 G4PixeCrossSectionHandler::G4PixeCrossSectionH
62 {
63 crossSections = 0;
64 interpolation = 0;
65 // Initialise with default values
66 Initialise(0,"","","",1.*keV,0.1*GeV,200,MeV
67 ActiveElements();
68 }
69
70
71 G4PixeCrossSectionHandler::G4PixeCrossSectionH
72 const G4String& modelK,
73 const G4String& modelL,
74 const G4String& modelM,
75 G4double minE,
76 G4double maxE,
77 G4int bins,
78 G4double unitE,
79 G4double unitData,
80 G4int minZ,
81 G4int maxZ)
82 : interpolation(algorithm), eMin(minE), eMax
83 unit1(unitE), unit2(unitData), zMin(minZ),
84 {
85 crossSections = 0;
86
87 crossModel.push_back(modelK);
88 crossModel.push_back(modelL);
89 crossModel.push_back(modelM);
90
91 //std::cout << "PixeCrossSectionHandler cons
92 // << crossModel[0]
93 // << std::endl;
94
95 ActiveElements();
96 }
97
98 G4PixeCrossSectionHandler::~G4PixeCrossSection
99 {
100 delete interpolation;
101 interpolation = 0;
102 std::map<G4int,G4IDataSet*,std::less<G4int>
103
104 for (pos = dataMap.begin(); pos != dataMap.e
105 {
106 // The following is a workaround for STL
107 // which does not support the standard a
108 // the syntax pos->second
109 // G4IDataSet* dataSet = pos->second;
110 G4IDataSet* dataSet = (*pos).second;
111 delete dataSet;
112 }
113
114 if (crossSections != 0)
115 {
116 std::size_t n = crossSections->size();
117 for (std::size_t i=0; i<n; ++i)
118 {
119 delete (*crossSections)[i];
120 }
121 delete crossSections;
122 crossSections = 0;
123 }
124 }
125
126 void G4PixeCrossSectionHandler::Initialise(G4I
127 const G4String& modelK,
128 const G4String& modelL,
129 const G4String& modelM,
130 G4double minE, G4double maxE,
131 G4int numberOfBins,
132 G4double unitE, G4double unitData
133 G4int minZ, G4int maxZ)
134 {
135 if (algorithm != 0)
136 {
137 delete interpolation;
138 interpolation = algorithm;
139 }
140 else
141 {
142 interpolation = CreateInterpolation();
143 }
144
145 eMin = minE;
146 eMax = maxE;
147 nBins = numberOfBins;
148 unit1 = unitE;
149 unit2 = unitData;
150 zMin = minZ;
151 zMax = maxZ;
152
153 crossModel.push_back(modelK);
154 crossModel.push_back(modelL);
155 crossModel.push_back(modelM);
156
157 }
158
159 void G4PixeCrossSectionHandler::PrintData() co
160 {
161 std::map<G4int,G4IDataSet*,std::less<G4int>
162
163 for (pos = dataMap.begin(); pos != dataMap.e
164 {
165 // The following is a workaround for STL
166 // which does not support the standard a
167 // the syntax pos->first or pos->second
168 // G4int z = pos->first;
169 // G4IDataSet* dataSet = pos->second;
170 G4int z = (*pos).first;
171 G4IDataSet* dataSet = (*pos).second;
172 G4cout << "---- Data set for Z = "
173 << z
174 << G4endl;
175 dataSet->PrintData();
176 G4cout << "-----------------------------
177 }
178 }
179
180 void G4PixeCrossSectionHandler::LoadShellData(
181 {
182 std::size_t nZ = activeZ.size();
183 for (std::size_t i=0; i<nZ; ++i)
184 {
185 G4int Z = (G4int) activeZ[i];
186 G4IInterpolator* algo = interpolation->C
187 G4IDataSet* dataSet = new G4PixeShellDat
188
189 // Degug printing
190 //std::cout << "PixeCrossSectionHandler:
191 // << Z
192 // << ", modelK = "
193 // << crossModel[0]
194 // << " fileName = "
195 // << fileName
196 // << std::endl;
197
198 dataSet->LoadData(fileName);
199 dataMap[Z] = dataSet;
200 }
201
202 // Build cross sections for materials if not
203 if (! crossSections)
204 {
205 BuildForMaterials();
206 }
207
208 }
209
210 void G4PixeCrossSectionHandler::Clear()
211 {
212 // Reset the map of data sets: remove the da
213 std::map<G4int,G4IDataSet*,std::less<G4int>
214
215 if(! dataMap.empty())
216 {
217 for (pos = dataMap.begin(); pos != dataM
218 {
219 // The following is a workaround for STL O
220 // which does not support the standard and
221 // the syntax pos->first or pos->second
222 // G4IDataSet* dataSet = pos->second;
223 G4IDataSet* dataSet = (*pos).second;
224 delete dataSet;
225 dataSet = 0;
226 G4int i = (*pos).first;
227 dataMap[i] = 0;
228 }
229 dataMap.clear();
230 }
231
232 activeZ.clear();
233 ActiveElements();
234 }
235
236 G4double G4PixeCrossSectionHandler::FindValue(
237 {
238 G4double value = 0.;
239
240 std::map<G4int,G4IDataSet*,std::less<G4int>
241 pos = dataMap.find(Z);
242 if (pos!= dataMap.end())
243 {
244 // The following is a workaround for STL
245 // which does not support the standard a
246 // the syntax pos->first or pos->second
247 // G4IDataSet* dataSet = pos->second;
248 G4IDataSet* dataSet = (*pos).second;
249 value = dataSet->FindValue(energy);
250 }
251 else
252 {
253 G4cout << "WARNING: G4PixeCrossSectionHa
254 << Z << G4endl;
255 }
256 return value;
257 }
258
259 G4double G4PixeCrossSectionHandler::FindValue(
260 G4int shellIndex) const
261 {
262 G4double value = 0.;
263
264 std::map<G4int,G4IDataSet*,std::less<G4int>
265 pos = dataMap.find(Z);
266 if (pos!= dataMap.end())
267 {
268 // The following is a workaround for STL
269 // which does not support the standard a
270 // the syntax pos->first or pos->second
271 // G4IDataSet* dataSet = pos->second;
272 G4IDataSet* dataSet = (*pos).second;
273 if (shellIndex >= 0)
274 {
275 G4int nComponents = (G4int)dataSet->Number
276 if(shellIndex < nComponents)
277 // The value is the cross section for sh
278 value = dataSet->GetComponent(shellIndex
279 else
280 {
281 G4cout << "WARNING: G4PixeCrossSection
282 << " shellIndex= " << shellIndex
283 << " for Z= "
284 << Z << G4endl;
285 }
286 } else {
287 value = dataSet->FindValue(energy);
288 }
289 }
290 else
291 {
292 G4cout << "WARNING: G4PixeCrossSectionHa
293 << Z << G4endl;
294 }
295 return value;
296 }
297
298
299 G4double G4PixeCrossSectionHandler::ValueForMa
300 G4double energy) const
301 {
302 G4double value = 0.;
303
304 const G4ElementVector* elementVector = mater
305 const G4double* nAtomsPerVolume = material->
306 std::size_t nElements = material->GetNumberO
307
308 for (std::size_t i=0 ; i<nElements ; ++i)
309 {
310 G4int Z = (G4int) (*elementVector)[i]->G
311 G4double elementValue = FindValue(Z,ener
312 G4double nAtomsVol = nAtomsPerVolume[i];
313 value += nAtomsVol * elementValue;
314 }
315
316 return value;
317 }
318
319 /*
320 G4IDataSet* G4PixeCrossSectionHandler::Build
321 {
322 // Builds a CompositeDataSet containing the
323 // in the material table
324
325 G4DataVector energyVector;
326 G4double dBin = std::log10(eMax/eMin) / nBin
327
328 for (G4int i=0; i<nBins+1; i++)
329 {
330 energyVector.push_back(std::pow(10., std::lo
331 }
332
333 // Factory method to build cross sections in
334 // related to the type of physics process
335
336 if (crossSections != 0)
337 { // Reset the list of cross sections
338 std::vector<G4IDataSet*>::iterator mat;
339 if (! crossSections->empty())
340 {
341 for (mat = crossSections->begin(); mat!= cro
342 {
343 G4IDataSet* set = *mat;
344 delete set;
345 set = 0;
346 }
347 crossSections->clear();
348 delete crossSections;
349 crossSections = 0;
350 }
351 }
352
353 crossSections = BuildCrossSectionsForMateria
354
355 if (crossSections == 0)
356 G4Exception("G4PixeCrossSectionHandler::Buil
357 "pii00000201",
358 FatalException,
359 "crossSections = 0");
360
361 G4IInterpolator* algo = CreateInterpolation(
362 G4IDataSet* materialSet = new G4CompositeDat
363
364 G4DataVector* energies;
365 G4DataVector* data;
366
367 const G4ProductionCutsTable* theCoupleTable=
368 G4ProductionCutsTable::GetProductionCutsTabl
369 std::size_t numOfCouples = theCoupleTable->G
370
371
372 for (std::size_t m=0; m<numOfCouples; m++)
373 {
374 energies = new G4DataVector;
375 data = new G4DataVector;
376 for (G4int bin=0; bin<nBins; bin++)
377 {
378 G4double energy = energyVector[bin];
379 energies->push_back(energy);
380 G4IDataSet* matCrossSet = (*crossSections)[m
381 G4double materialCrossSection = 0.0;
382 G4int nElm = matCrossSet->NumberOfComponents
383 for(G4int j=0; j<nElm; j++) {
384 materialCrossSection += matCrossSet->GetComp
385 }
386
387 if (materialCrossSection > 0.)
388 {
389 data->push_back(1./materialCrossSection);
390 }
391 else
392 {
393 data->push_back(DBL_MAX);
394 }
395 }
396 G4IInterpolator* algo = CreateInterpolation(
397 G4IDataSet* dataSet = new G4DataSet(m,energi
398 materialSet->AddComponent(dataSet);
399 }
400
401 return materialSet;
402 }
403
404 */
405
406 void G4PixeCrossSectionHandler::BuildForMateri
407 {
408 // Builds a CompositeDataSet containing the
409 // in the material table
410
411 G4DataVector energyVector;
412 G4double dBin = std::log10(eMax/eMin) / nBin
413
414 for (G4int i=0; i<nBins+1; i++)
415 {
416 energyVector.push_back(std::pow(10., std
417 }
418
419 if (crossSections != 0)
420 { // Reset the list of cross sections
421 std::vector<G4IDataSet*>::iterator mat;
422 if (! crossSections->empty())
423 {
424 for (mat = crossSections->begin(); mat!= c
425 {
426 G4IDataSet* set = *mat;
427 delete set;
428 set = 0;
429 }
430 crossSections->clear();
431 delete crossSections;
432 crossSections = 0;
433 }
434 }
435
436 crossSections = BuildCrossSectionsForMateria
437
438 if (crossSections == 0)
439 G4Exception("G4PixeCrossSectionHandler::Bu
440 "pii00000210",
441 FatalException,
442 ", crossSections = 0");
443
444 return;
445 }
446
447
448 G4int G4PixeCrossSectionHandler::SelectRandomA
449 G4double e) const
450 {
451 // Select randomly an element within the mat
452 // determined by the cross sections in the d
453
454 G4int nElements = (G4int)material->GetNumber
455
456 // Special case: the material consists of on
457 if (nElements == 1)
458 {
459 G4int Z = (G4int) material->GetZ();
460 return Z;
461 }
462
463 // Composite material
464
465 const G4ElementVector* elementVector = mater
466 std::size_t materialIndex = material->GetInd
467
468 G4IDataSet* materialSet = (*crossSections)[m
469 G4double materialCrossSection0 = 0.0;
470 G4DataVector cross;
471 cross.clear();
472 for ( G4int i=0; i < nElements; ++i )
473 {
474 G4double cr = materialSet->GetComponent(
475 materialCrossSection0 += cr;
476 cross.push_back(materialCrossSection0);
477 }
478
479 G4double random = G4UniformRand() * material
480
481 for (G4int k=0 ; k < nElements ; ++k )
482 {
483 if (random <= cross[k]) return (G4int) (
484 }
485 // It should never get here
486 return 0;
487 }
488
489 /*
490 const G4Element* G4PixeCrossSectionHandler::
491 G4double e) const
492 {
493 // Select randomly an element within the mat
494 // by the cross sections in the data set
495
496 const G4Material* material = couple->GetMate
497 G4Element* nullElement = 0;
498 G4int nElements = material->GetNumberOfEleme
499 const G4ElementVector* elementVector = mater
500
501 // Special case: the material consists of on
502 if (nElements == 1)
503 {
504 G4Element* element = (*elementVector)[0];
505 return element;
506 }
507 else
508 {
509 // Composite material
510
511 std::size_t materialIndex = couple->GetIndex
512
513 G4IDataSet* materialSet = (*crossSections)[m
514 G4double materialCrossSection0 = 0.0;
515 G4DataVector cross;
516 cross.clear();
517 for (G4int i=0; i<nElements; i++)
518 {
519 G4double cr = materialSet->GetComponent(i)->
520 materialCrossSection0 += cr;
521 cross.push_back(materialCrossSection0);
522 }
523
524 G4double random = G4UniformRand() * material
525
526 for (G4int k=0 ; k < nElements ; k++ )
527 {
528 if (random <= cross[k]) return (*elementVect
529 }
530 // It should never end up here
531 G4cout << "G4PixeCrossSectionHandler::Select
532 return nullElement;
533 }
534 }
535 */
536
537
538 G4int G4PixeCrossSectionHandler::SelectRandomS
539 {
540 // Select randomly a shell, according to the
541 // in the data set
542
543 // Note for later improvement: it would be u
544 // used shells to improve performance
545
546 G4int shell = 0;
547
548 G4double totCrossSection = FindValue(Z,e);
549 G4double random = G4UniformRand() * totCross
550 G4double partialSum = 0.;
551
552 G4IDataSet* dataSet = nullptr;
553 auto pos = dataMap.find(Z);
554 // The following is a workaround for STL Obj
555 // which does not support the standard and d
556 // the syntax pos->first or pos->second
557 // if (pos != dataMap.end()) dataSet = pos->
558 if (pos != dataMap.end()) dataSet = (*pos).s
559
560 if (dataSet != nullptr) {
561 G4int nShells = (G4int)dataSet->NumberOfCo
562 for (G4int i=0; i<nShells; ++i)
563 {
564 const G4IDataSet* shellDataSet = dataSet
565 if (shellDataSet != 0)
566 {
567 G4double value = shellDataSet->FindVal
568 partialSum += value;
569 if (random <= partialSum) return i;
570 }
571 }
572 }
573 // It should never get here
574 return shell;
575 }
576
577 void G4PixeCrossSectionHandler::ActiveElements
578 {
579 const G4MaterialTable* materialTable = G4Mat
580 if (materialTable == 0)
581 G4Exception("G4PixeCrossSectionHandler::Ac
582 "pii00000220",
583 FatalException,
584 "no MaterialTable found");
585
586 std::size_t nMaterials = G4Material::GetNumb
587
588 for (std::size_t mat=0; mat<nMaterials; ++ma
589 {
590 const G4Material* material= (*materialTa
591 const G4ElementVector* elementVector = m
592 const std::size_t nElements = material->
593
594 for (std::size_t iEl=0; iEl<nElements; +
595 {
596 G4double Z = (*elementVector)[iEl]->GetZ()
597 if (!(activeZ.contains(Z)) && Z >= zMin &&
598 {
599 activeZ.push_back(Z);
600 }
601 }
602 }
603 }
604
605 G4IInterpolator* G4PixeCrossSectionHandler::Cr
606 {
607 G4IInterpolator* algorithm = new G4LogLogInt
608 return algorithm;
609 }
610
611 G4int G4PixeCrossSectionHandler::NumberOfCompo
612 {
613 G4int n = 0;
614
615 std::map<G4int,G4IDataSet*,std::less<G4int>
616 pos = dataMap.find(Z);
617 if (pos!= dataMap.end())
618 {
619 G4IDataSet* dataSet = (*pos).second;
620 n = (G4int)dataSet->NumberOfComponents()
621 }
622 else
623 {
624 G4cout << "WARNING: G4PixeCrossSectionHa
625 << "find Z = "
626 << Z << G4endl;
627 }
628 return n;
629 }
630
631
632 std::vector<G4IDataSet*>*
633 G4PixeCrossSectionHandler::BuildCrossSectionsF
634 {
635 G4DataVector* energies;
636 G4DataVector* data;
637
638 std::vector<G4IDataSet*>* matCrossSections =
639
640 //const G4ProductionCutsTable* theCoupleTabl
641 //std::size_t numOfCouples = theCoupleTable-
642
643 std::size_t nOfBins = energyVector.size();
644 const auto interpolationAlgo = std::unique_p
645
646 const G4MaterialTable* materialTable = G4Mat
647 if (materialTable == 0)
648 G4Exception("G4PixeCrossSectionHandler::Bu
649 "pii00000230",
650 FatalException,
651 "no MaterialTable found");
652
653 std::size_t nMaterials = G4Material::GetNumb
654
655 for (std::size_t mat=0; mat<nMaterials; ++ma
656 {
657 const G4Material* material = (*materialT
658 G4int nElements = (G4int)material->GetNu
659 const G4ElementVector* elementVector = m
660 const G4double* nAtomsPerVolume = materi
661
662 G4IInterpolator* algo = interpolationAlg
663
664 G4IDataSet* setForMat = new G4CompositeD
665
666 for (G4int i=0; i<nElements; ++i) {
667
668 G4int Z = (G4int) (*elementVector)[i]-
669 G4double density = nAtomsPerVolume[i];
670
671 energies = new G4DataVector;
672 data = new G4DataVector;
673
674
675 for (std::size_t bin=0; bin<nOfBins; +
676 {
677 G4double e = energyVector[bin];
678 energies->push_back(e);
679 G4double cross = 0.;
680 if (Z >= zMin && Z <= zMax) cross = dens
681 data->push_back(cross);
682 }
683
684 G4IInterpolator* algo1 = interpolation
685 G4IDataSet* elSet = new G4DataSet(i,en
686 setForMat->AddComponent(elSet);
687 }
688
689 matCrossSections->push_back(setForMat);
690 }
691 return matCrossSections;
692 }
693
694
695 G4double G4PixeCrossSectionHandler::Microscopi
696 G4double kineticEnergy,
697 G4double Z,
698 G4double deltaCut) const
699 {
700 // Cross section formula is OK for spin=0, 1
701 // Calculates the microscopic cross section
702 // Formula documented in Geant4 Phys. Ref. M
703 // ( it is called for elements, AtomicNumber
704
705 G4double cross = 0.;
706
707 // Particle mass and energy
708 G4double particleMass = particleDef->GetPD
709 G4double energy = kineticEnergy + particle
710
711 // Some kinematics
712 G4double gamma = energy / particleMass;
713 G4double beta2 = 1. - 1. / (gamma * gamma);
714 G4double var = electron_mass_c2 / particleMa
715 G4double tMax = 2. * electron_mass_c2 * (gam
716
717 // Calculate the total cross section
718
719 if ( tMax > deltaCut )
720 {
721 var = deltaCut / tMax;
722 cross = (1. - var * (1. - beta2 * std::l
723
724 G4double spin = particleDef->GetPDGSpin(
725
726 // +term for spin=1/2 particle
727 if (spin == 0.5)
728 {
729 cross += 0.5 * (tMax - deltaCut) / (energ
730 }
731 // +term for spin=1 particle
732 else if (spin > 0.9 )
733 {
734 cross += -std::log(var) / (3.*deltaCut) +
735 ((5.+1./var)*0.25 /(energy*energy) - bet
736 }
737 cross *= twopi_mc2_rcl2 * Z / beta2 ;
738 }
739
740 //std::cout << "Microscopic = " << cross/bar
741 // << ", e = " << kineticEnergy/MeV <<std
742
743 return cross;
744 }
745
746