| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 //
30 // History:
31 // -----------
32 // 1 Aug 2001 MGP Created
33 // 09 Oct 2001 VI Add FindValue with 3 parameters
34 // + NumberOfComponents
35 // 19 Jul 2002 VI Create composite data set for material
36 // 21 Jan 2003 VI Cut per region
37 //
38 // -------------------------------------------------------------------
39
40 #include "G4RDVCrossSectionHandler.hh"
41 #include "G4RDVDataSetAlgorithm.hh"
42 #include "G4RDLogLogInterpolation.hh"
43 #include "G4RDVEMDataSet.hh"
44 #include "G4RDEMDataSet.hh"
45 #include "G4RDCompositeEMDataSet.hh"
46 #include "G4RDShellEMDataSet.hh"
47 #include "G4ProductionCutsTable.hh"
48 #include "G4Material.hh"
49 #include "G4Element.hh"
50 #include "Randomize.hh"
51 #include <map>
52 #include <vector>
53 #include <fstream>
54 #include <sstream>
55
56
57 G4RDVCrossSectionHandler::G4RDVCrossSectionHandler()
58 {
59 crossSections = 0;
60 interpolation = 0;
61 Initialise();
62 ActiveElements();
63 }
64
65
66 G4RDVCrossSectionHandler::G4RDVCrossSectionHandler(G4RDVDataSetAlgorithm* algorithm,
67 G4double minE,
68 G4double maxE,
69 G4int bins,
70 G4double unitE,
71 G4double unitData,
72 G4int minZ,
73 G4int maxZ)
74 : interpolation(algorithm), eMin(minE), eMax(maxE), nBins(bins),
75 unit1(unitE), unit2(unitData), zMin(minZ), zMax(maxZ)
76 {
77 crossSections = 0;
78 ActiveElements();
79 }
80
81 G4RDVCrossSectionHandler::~G4RDVCrossSectionHandler()
82 {
83 delete interpolation;
84 interpolation = 0;
85 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
86
87 for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
88 {
89 // The following is a workaround for STL ObjectSpace implementation,
90 // which does not support the standard and does not accept
91 // the syntax pos->second
92 // G4RDVEMDataSet* dataSet = pos->second;
93 G4RDVEMDataSet* dataSet = (*pos).second;
94 delete dataSet;
95 }
96
97 if (crossSections != 0)
98 {
99 size_t n = crossSections->size();
100 for (size_t i=0; i<n; i++)
101 {
102 delete (*crossSections)[i];
103 }
104 delete crossSections;
105 crossSections = 0;
106 }
107 }
108
109 void G4RDVCrossSectionHandler::Initialise(G4RDVDataSetAlgorithm* algorithm,
110 G4double minE, G4double maxE,
111 G4int numberOfBins,
112 G4double unitE, G4double unitData,
113 G4int minZ, G4int maxZ)
114 {
115 if (algorithm != 0)
116 {
117 delete interpolation;
118 interpolation = algorithm;
119 }
120 else
121 {
122 interpolation = CreateInterpolation();
123 }
124
125 eMin = minE;
126 eMax = maxE;
127 nBins = numberOfBins;
128 unit1 = unitE;
129 unit2 = unitData;
130 zMin = minZ;
131 zMax = maxZ;
132 }
133
134 void G4RDVCrossSectionHandler::PrintData() const
135 {
136 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
137
138 for (pos = dataMap.begin(); pos != dataMap.end(); pos++)
139 {
140 // The following is a workaround for STL ObjectSpace implementation,
141 // which does not support the standard and does not accept
142 // the syntax pos->first or pos->second
143 // G4int z = pos->first;
144 // G4RDVEMDataSet* dataSet = pos->second;
145 G4int z = (*pos).first;
146 G4RDVEMDataSet* dataSet = (*pos).second;
147 G4cout << "---- Data set for Z = "
148 << z
149 << G4endl;
150 dataSet->PrintData();
151 G4cout << "--------------------------------------------------" << G4endl;
152 }
153 }
154
155 void G4RDVCrossSectionHandler::LoadData(const G4String& fileName)
156 {
157 size_t nZ = activeZ.size();
158 for (size_t i=0; i<nZ; i++)
159 {
160 G4int Z = (G4int) activeZ[i];
161
162 // Build the complete string identifying the file with the data set
163
164 const char* path = G4FindDataDir("G4LEDATA");
165 if (!path)
166 {
167 G4String excep = "G4LEDATA environment variable not set!";
168 G4Exception("G4RDVCrossSectionHandler::LoadData()",
169 "InvalidSetup", FatalException, excep);
170 }
171
172 std::ostringstream ost;
173 ost << path << '/' << fileName << Z << ".dat";
174 std::ifstream file(ost.str().c_str());
175 std::filebuf* lsdp = file.rdbuf();
176
177 if (! (lsdp->is_open()) )
178 {
179 G4String excep = "Data file: " + ost.str() + " not found!";
180 G4Exception("G4RDVCrossSectionHandler::LoadData()",
181 "DataNotFound", FatalException, excep);
182 }
183 G4double a = 0;
184 G4int k = 1;
185 G4DataVector* energies = new G4DataVector;
186 G4DataVector* data = new G4DataVector;
187 do
188 {
189 file >> a;
190 G4int nColumns = 2;
191 // The file is organized into two columns:
192 // 1st column is the energy
193 // 2nd column is the corresponding value
194 // The file terminates with the pattern: -1 -1
195 // -2 -2
196 if (a == -1 || a == -2)
197 {
198 }
199 else
200 {
201 if (k%nColumns != 0)
202 {
203 G4double e = a * unit1;
204 energies->push_back(e);
205 k++;
206 }
207 else if (k%nColumns == 0)
208 {
209 G4double value = a * unit2;
210 data->push_back(value);
211 k = 1;
212 }
213 }
214 } while (a != -2); // end of file
215
216 file.close();
217 G4RDVDataSetAlgorithm* algo = interpolation->Clone();
218 G4RDVEMDataSet* dataSet = new G4RDEMDataSet(Z,energies,data,algo);
219 dataMap[Z] = dataSet;
220 }
221 }
222
223 void G4RDVCrossSectionHandler::LoadShellData(const G4String& fileName)
224 {
225 size_t nZ = activeZ.size();
226 for (size_t i=0; i<nZ; i++)
227 {
228 G4int Z = (G4int) activeZ[i];
229
230 // Riccardo Capra <capra@ge.infn.it>: PLEASE CHECK THE FOLLOWING PIECE OF CODE
231 // "energies" AND "data" G4DataVector ARE ALLOCATED, FILLED IN AND NEVER USED OR
232 // DELETED. WHATSMORE LOADING FILE OPERATIONS WERE DONE BY G4RDShellEMDataSet
233 // EVEN BEFORE THE CHANGES I DID ON THIS FILE. SO THE FOLLOWING CODE IN MY
234 // OPINION SHOULD BE USELESS AND SHOULD PRODUCE A MEMORY LEAK.
235
236 // Build the complete string identifying the file with the data set
237
238 const char* path = G4FindDataDir("G4LEDATA");
239 if (!path)
240 {
241 G4String excep = "G4LEDATA environment variable not set!";
242 G4Exception("G4RDVCrossSectionHandler::LoadShellData()",
243 "InvalidSetup", FatalException, excep);
244 }
245
246 std::ostringstream ost;
247
248 ost << path << '/' << fileName << Z << ".dat";
249
250 std::ifstream file(ost.str().c_str());
251 std::filebuf* lsdp = file.rdbuf();
252
253 if (! (lsdp->is_open()) )
254 {
255 G4String excep = "Data file: " + ost.str() + " not found!";
256 G4Exception("G4RDVCrossSectionHandler::LoadShellData()",
257 "DataNotFound", FatalException, excep);
258 }
259 G4double a = 0;
260 G4int k = 1;
261 G4DataVector* energies = new G4DataVector;
262 G4DataVector* data = new G4DataVector;
263 do
264 {
265 file >> a;
266 G4int nColumns = 2;
267 // The file is organized into two columns:
268 // 1st column is the energy
269 // 2nd column is the corresponding value
270 // The file terminates with the pattern: -1 -1
271 // -2 -2
272 if (a == -1 || a == -2)
273 {
274 }
275 else
276 {
277 if (k%nColumns != 0)
278 {
279 G4double e = a * unit1;
280 energies->push_back(e);
281 k++;
282 }
283 else if (k%nColumns == 0)
284 {
285 G4double value = a * unit2;
286 data->push_back(value);
287 k = 1;
288 }
289 }
290 } while (a != -2); // end of file
291
292 file.close();
293
294 // Riccardo Capra <capra@ge.infn.it>: END OF CODE THAT IN MY OPINION SHOULD BE
295 // REMOVED.
296
297 G4RDVDataSetAlgorithm* algo = interpolation->Clone();
298 G4RDVEMDataSet* dataSet = new G4RDShellEMDataSet(Z, algo);
299 dataSet->LoadData(fileName);
300 dataMap[Z] = dataSet;
301 }
302 }
303
304 void G4RDVCrossSectionHandler::Clear()
305 {
306 // Reset the map of data sets: remove the data sets from the map
307 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
308
309 if(! dataMap.empty())
310 {
311 for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
312 {
313 // The following is a workaround for STL ObjectSpace implementation,
314 // which does not support the standard and does not accept
315 // the syntax pos->first or pos->second
316 // G4RDVEMDataSet* dataSet = pos->second;
317 G4RDVEMDataSet* dataSet = (*pos).second;
318 delete dataSet;
319 dataSet = 0;
320 G4int i = (*pos).first;
321 dataMap[i] = 0;
322 }
323 dataMap.clear();
324 }
325
326 activeZ.clear();
327 ActiveElements();
328 }
329
330 G4double G4RDVCrossSectionHandler::FindValue(G4int Z, G4double energy) const
331 {
332 G4double value = 0.;
333
334 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
335 pos = dataMap.find(Z);
336 if (pos!= dataMap.end())
337 {
338 // The following is a workaround for STL ObjectSpace implementation,
339 // which does not support the standard and does not accept
340 // the syntax pos->first or pos->second
341 // G4RDVEMDataSet* dataSet = pos->second;
342 G4RDVEMDataSet* dataSet = (*pos).second;
343 value = dataSet->FindValue(energy);
344 }
345 else
346 {
347 G4cout << "WARNING: G4RDVCrossSectionHandler::FindValue did not find Z = "
348 << Z << G4endl;
349 }
350 return value;
351 }
352
353 G4double G4RDVCrossSectionHandler::FindValue(G4int Z, G4double energy,
354 G4int shellIndex) const
355 {
356 G4double value = 0.;
357
358 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
359 pos = dataMap.find(Z);
360 if (pos!= dataMap.end())
361 {
362 // The following is a workaround for STL ObjectSpace implementation,
363 // which does not support the standard and does not accept
364 // the syntax pos->first or pos->second
365 // G4RDVEMDataSet* dataSet = pos->second;
366 G4RDVEMDataSet* dataSet = (*pos).second;
367 if (shellIndex >= 0)
368 {
369 G4int nComponents = dataSet->NumberOfComponents();
370 if(shellIndex < nComponents)
371 // - MGP - Why doesn't it use G4RDVEMDataSet::FindValue directly?
372 value = dataSet->GetComponent(shellIndex)->FindValue(energy);
373 else
374 {
375 G4cout << "WARNING: G4RDVCrossSectionHandler::FindValue did not find"
376 << " shellIndex= " << shellIndex
377 << " for Z= "
378 << Z << G4endl;
379 }
380 } else {
381 value = dataSet->FindValue(energy);
382 }
383 }
384 else
385 {
386 G4cout << "WARNING: G4RDVCrossSectionHandler::FindValue did not find Z = "
387 << Z << G4endl;
388 }
389 return value;
390 }
391
392
393 G4double G4RDVCrossSectionHandler::ValueForMaterial(const G4Material* material,
394 G4double energy) const
395 {
396 G4double value = 0.;
397
398 const G4ElementVector* elementVector = material->GetElementVector();
399 const G4double* nAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
400 G4int nElements = material->GetNumberOfElements();
401
402 for (G4int i=0 ; i<nElements ; i++)
403 {
404 G4int Z = (G4int) (*elementVector)[i]->GetZ();
405 G4double elementValue = FindValue(Z,energy);
406 G4double nAtomsVol = nAtomsPerVolume[i];
407 value += nAtomsVol * elementValue;
408 }
409
410 return value;
411 }
412
413
414 G4RDVEMDataSet* G4RDVCrossSectionHandler::BuildMeanFreePathForMaterials(const G4DataVector* energyCuts)
415 {
416 // Builds a CompositeDataSet containing the mean free path for each material
417 // in the material table
418
419 G4DataVector energyVector;
420 G4double dBin = std::log10(eMax/eMin) / nBins;
421
422 for (G4int i=0; i<nBins+1; i++)
423 {
424 energyVector.push_back(std::pow(10., std::log10(eMin)+i*dBin));
425 }
426
427 // Factory method to build cross sections in derived classes,
428 // related to the type of physics process
429
430 if (crossSections != 0)
431 { // Reset the list of cross sections
432 std::vector<G4RDVEMDataSet*>::iterator mat;
433 if (! crossSections->empty())
434 {
435 for (mat = crossSections->begin(); mat!= crossSections->end(); ++mat)
436 {
437 G4RDVEMDataSet* set = *mat;
438 delete set;
439 set = 0;
440 }
441 crossSections->clear();
442 delete crossSections;
443 crossSections = 0;
444 }
445 }
446
447 crossSections = BuildCrossSectionsForMaterials(energyVector,energyCuts);
448
449 if (crossSections == 0)
450 G4Exception("G4RDVCrossSectionHandler::BuildMeanFreePathForMaterials()",
451 "InvalidCondition", FatalException, "CrossSections = 0!");
452
453 G4RDVDataSetAlgorithm* algo = CreateInterpolation();
454 G4RDVEMDataSet* materialSet = new G4RDCompositeEMDataSet(algo);
455
456 G4DataVector* energies;
457 G4DataVector* data;
458
459 const G4ProductionCutsTable* theCoupleTable=
460 G4ProductionCutsTable::GetProductionCutsTable();
461 size_t numOfCouples = theCoupleTable->GetTableSize();
462
463
464 for (size_t m=0; m<numOfCouples; m++)
465 {
466 energies = new G4DataVector;
467 data = new G4DataVector;
468 for (G4int bin=0; bin<nBins; bin++)
469 {
470 G4double energy = energyVector[bin];
471 energies->push_back(energy);
472 G4RDVEMDataSet* matCrossSet = (*crossSections)[m];
473 G4double materialCrossSection = 0.0;
474 G4int nElm = matCrossSet->NumberOfComponents();
475 for(G4int j=0; j<nElm; j++) {
476 materialCrossSection += matCrossSet->GetComponent(j)->FindValue(energy);
477 }
478
479 if (materialCrossSection > 0.)
480 {
481 data->push_back(1./materialCrossSection);
482 }
483 else
484 {
485 data->push_back(DBL_MAX);
486 }
487 }
488 G4RDVDataSetAlgorithm* algo = CreateInterpolation();
489 G4RDVEMDataSet* dataSet = new G4RDEMDataSet(m,energies,data,algo,1.,1.);
490 materialSet->AddComponent(dataSet);
491 }
492
493 return materialSet;
494 }
495
496 G4int G4RDVCrossSectionHandler::SelectRandomAtom(const G4MaterialCutsCouple* couple,
497 G4double e) const
498 {
499 // Select randomly an element within the material, according to the weight
500 // determined by the cross sections in the data set
501
502 const G4Material* material = couple->GetMaterial();
503 G4int nElements = material->GetNumberOfElements();
504
505 // Special case: the material consists of one element
506 if (nElements == 1)
507 {
508 G4int Z = (G4int) material->GetZ();
509 return Z;
510 }
511
512 // Composite material
513
514 const G4ElementVector* elementVector = material->GetElementVector();
515 size_t materialIndex = couple->GetIndex();
516
517 G4RDVEMDataSet* materialSet = (*crossSections)[materialIndex];
518 G4double materialCrossSection0 = 0.0;
519 G4DataVector cross;
520 cross.clear();
521 for ( G4int i=0; i < nElements; i++ )
522 {
523 G4double cr = materialSet->GetComponent(i)->FindValue(e);
524 materialCrossSection0 += cr;
525 cross.push_back(materialCrossSection0);
526 }
527
528 G4double random = G4UniformRand() * materialCrossSection0;
529
530 for (G4int k=0 ; k < nElements ; k++ )
531 {
532 if (random <= cross[k]) return (G4int) (*elementVector)[k]->GetZ();
533 }
534 // It should never get here
535 return 0;
536 }
537
538 const G4Element* G4RDVCrossSectionHandler::SelectRandomElement(const G4MaterialCutsCouple* couple,
539 G4double e) const
540 {
541 // Select randomly an element within the material, according to the weight determined
542 // by the cross sections in the data set
543
544 const G4Material* material = couple->GetMaterial();
545 G4Element* nullElement = 0;
546 G4int nElements = material->GetNumberOfElements();
547 const G4ElementVector* elementVector = material->GetElementVector();
548
549 // Special case: the material consists of one element
550 if (nElements == 1)
551 {
552 G4Element* element = (*elementVector)[0];
553 return element;
554 }
555 else
556 {
557 // Composite material
558
559 size_t materialIndex = couple->GetIndex();
560
561 G4RDVEMDataSet* materialSet = (*crossSections)[materialIndex];
562 G4double materialCrossSection0 = 0.0;
563 G4DataVector cross;
564 cross.clear();
565 for (G4int i=0; i<nElements; i++)
566 {
567 G4double cr = materialSet->GetComponent(i)->FindValue(e);
568 materialCrossSection0 += cr;
569 cross.push_back(materialCrossSection0);
570 }
571
572 G4double random = G4UniformRand() * materialCrossSection0;
573
574 for (G4int k=0 ; k < nElements ; k++ )
575 {
576 if (random <= cross[k]) return (*elementVector)[k];
577 }
578 // It should never end up here
579 G4cout << "G4RDVCrossSectionHandler::SelectRandomElement - no element found" << G4endl;
580 return nullElement;
581 }
582 }
583
584 G4int G4RDVCrossSectionHandler::SelectRandomShell(G4int Z, G4double e) const
585 {
586 // Select randomly a shell, according to the weight determined by the cross sections
587 // in the data set
588
589 // Note for later improvement: it would be useful to add a cache mechanism for already
590 // used shells to improve performance
591
592 G4int shell = 0;
593
594 G4double totCrossSection = FindValue(Z,e);
595 G4double random = G4UniformRand() * totCrossSection;
596 G4double partialSum = 0.;
597
598 G4RDVEMDataSet* dataSet = 0;
599 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
600 pos = dataMap.find(Z);
601 // The following is a workaround for STL ObjectSpace implementation,
602 // which does not support the standard and does not accept
603 // the syntax pos->first or pos->second
604 // if (pos != dataMap.end()) dataSet = pos->second;
605 if (pos != dataMap.end()) dataSet = (*pos).second;
606
607 size_t nShells = dataSet->NumberOfComponents();
608 for (size_t i=0; i<nShells; i++)
609 {
610 const G4RDVEMDataSet* shellDataSet = dataSet->GetComponent(i);
611 if (shellDataSet != 0)
612 {
613 G4double value = shellDataSet->FindValue(e);
614 partialSum += value;
615 if (random <= partialSum) return i;
616 }
617 }
618 // It should never get here
619 return shell;
620 }
621
622 void G4RDVCrossSectionHandler::ActiveElements()
623 {
624 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
625 if (materialTable == 0)
626 G4Exception("G4RDVCrossSectionHandler::ActiveElements",
627 "InvalidSetup", FatalException, "No MaterialTable found!");
628
629 G4int nMaterials = G4Material::GetNumberOfMaterials();
630
631 for (G4int m=0; m<nMaterials; m++)
632 {
633 const G4Material* material= (*materialTable)[m];
634 const G4ElementVector* elementVector = material->GetElementVector();
635 const G4int nElements = material->GetNumberOfElements();
636
637 for (G4int iEl=0; iEl<nElements; iEl++)
638 {
639 G4Element* element = (*elementVector)[iEl];
640 G4double Z = element->GetZ();
641 if (!(activeZ.contains(Z)) && Z >= zMin && Z <= zMax)
642 {
643 activeZ.push_back(Z);
644 }
645 }
646 }
647 }
648
649 G4RDVDataSetAlgorithm* G4RDVCrossSectionHandler::CreateInterpolation()
650 {
651 G4RDVDataSetAlgorithm* algorithm = new G4RDLogLogInterpolation;
652 return algorithm;
653 }
654
655 G4int G4RDVCrossSectionHandler::NumberOfComponents(G4int Z) const
656 {
657 G4int n = 0;
658
659 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
660 pos = dataMap.find(Z);
661 if (pos!= dataMap.end())
662 {
663 G4RDVEMDataSet* dataSet = (*pos).second;
664 n = dataSet->NumberOfComponents();
665 }
666 else
667 {
668 G4cout << "WARNING: G4RDVCrossSectionHandler::NumberOfComponents did not "
669 << "find Z = "
670 << Z << G4endl;
671 }
672 return n;
673 }
674
675
676