| 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 // 1 Aug 2001 MGP Created
33 // 09 Oct 2001 VI Add FindValue with
34 // + NumberOfComponen
35 // 19 Jul 2002 VI Create composite d
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::G4RDVCrossSectionHan
58 {
59 crossSections = 0;
60 interpolation = 0;
61 Initialise();
62 ActiveElements();
63 }
64
65
66 G4RDVCrossSectionHandler::G4RDVCrossSectionHan
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
75 unit1(unitE), unit2(unitData), zMin(minZ),
76 {
77 crossSections = 0;
78 ActiveElements();
79 }
80
81 G4RDVCrossSectionHandler::~G4RDVCrossSectionHa
82 {
83 delete interpolation;
84 interpolation = 0;
85 std::map<G4int,G4RDVEMDataSet*,std::less<G4i
86
87 for (pos = dataMap.begin(); pos != dataMap.e
88 {
89 // The following is a workaround for STL
90 // which does not support the standard a
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(G4RD
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() con
135 {
136 std::map<G4int,G4RDVEMDataSet*,std::less<G4i
137
138 for (pos = dataMap.begin(); pos != dataMap.e
139 {
140 // The following is a workaround for STL
141 // which does not support the standard a
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 << "-----------------------------
152 }
153 }
154
155 void G4RDVCrossSectionHandler::LoadData(const
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
163
164 const char* path = G4FindDataDir("G4LEDA
165 if (!path)
166 {
167 G4String excep = "G4LEDATA environment var
168 G4Exception("G4RDVCrossSectionHandler::Loa
169 "InvalidSetup", FatalExc
170 }
171
172 std::ostringstream ost;
173 ost << path << '/' << fileName << Z << "
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()
180 G4Exception("G4RDVCrossSectionHandler::Loa
181 "DataNotFound", FatalExc
182 }
183 G4double a = 0;
184 G4int k = 1;
185 G4DataVector* energies = new G4DataVecto
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: -
195 // -
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 = interpolat
218 G4RDVEMDataSet* dataSet = new G4RDEMData
219 dataMap[Z] = dataSet;
220 }
221 }
222
223 void G4RDVCrossSectionHandler::LoadShellData(c
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>: PL
231 // "energies" AND "data" G4DataVector AR
232 // DELETED. WHATSMORE LOADING FILE OPERA
233 // EVEN BEFORE THE CHANGES I DID ON THIS
234 // OPINION SHOULD BE USELESS AND SHOULD
235
236 // Build the complete string identifying
237
238 const char* path = G4FindDataDir("G4LEDA
239 if (!path)
240 {
241 G4String excep = "G4LEDATA environment var
242 G4Exception("G4RDVCrossSectionHandler::Loa
243 "InvalidSetup", FatalExc
244 }
245
246 std::ostringstream ost;
247
248 ost << path << '/' << fileName << Z << "
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()
256 G4Exception("G4RDVCrossSectionHandler::Loa
257 "DataNotFound", FatalExc
258 }
259 G4double a = 0;
260 G4int k = 1;
261 G4DataVector* energies = new G4DataVecto
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: -
271 // -
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>: EN
295 // REMOVED.
296
297 G4RDVDataSetAlgorithm* algo = interpolat
298 G4RDVEMDataSet* dataSet = new G4RDShellE
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 da
307 std::map<G4int,G4RDVEMDataSet*,std::less<G4i
308
309 if(! dataMap.empty())
310 {
311 for (pos = dataMap.begin(); pos != dat
312 {
313 // The following is a workaround for STL O
314 // which does not support the standard and
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(G
331 {
332 G4double value = 0.;
333
334 std::map<G4int,G4RDVEMDataSet*,std::less<G4i
335 pos = dataMap.find(Z);
336 if (pos!= dataMap.end())
337 {
338 // The following is a workaround for STL
339 // which does not support the standard a
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: G4RDVCrossSectionHan
348 << Z << G4endl;
349 }
350 return value;
351 }
352
353 G4double G4RDVCrossSectionHandler::FindValue(G
354 G4i
355 {
356 G4double value = 0.;
357
358 std::map<G4int,G4RDVEMDataSet*,std::less<G4i
359 pos = dataMap.find(Z);
360 if (pos!= dataMap.end())
361 {
362 // The following is a workaround for STL
363 // which does not support the standard a
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->NumberOfCompo
370 if(shellIndex < nComponents)
371 // - MGP - Why doesn't it use G4RDVEMDat
372 value = dataSet->GetComponent(shellIndex
373 else
374 {
375 G4cout << "WARNING: G4RDVCrossSectionH
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: G4RDVCrossSectionHan
387 << Z << G4endl;
388 }
389 return value;
390 }
391
392
393 G4double G4RDVCrossSectionHandler::ValueForMat
394 G4double energy) const
395 {
396 G4double value = 0.;
397
398 const G4ElementVector* elementVector = mater
399 const G4double* nAtomsPerVolume = material->
400 G4int nElements = material->GetNumberOfEleme
401
402 for (G4int i=0 ; i<nElements ; i++)
403 {
404 G4int Z = (G4int) (*elementVector)[i]->G
405 G4double elementValue = FindValue(Z,ener
406 G4double nAtomsVol = nAtomsPerVolume[i];
407 value += nAtomsVol * elementValue;
408 }
409
410 return value;
411 }
412
413
414 G4RDVEMDataSet* G4RDVCrossSectionHandler::Buil
415 {
416 // Builds a CompositeDataSet containing the
417 // in the material table
418
419 G4DataVector energyVector;
420 G4double dBin = std::log10(eMax/eMin) / nBin
421
422 for (G4int i=0; i<nBins+1; i++)
423 {
424 energyVector.push_back(std::pow(10., std
425 }
426
427 // Factory method to build cross sections in
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 m
433 if (! crossSections->empty())
434 {
435 for (mat = crossSections->begin(); mat!= c
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 = BuildCrossSectionsForMateria
448
449 if (crossSections == 0)
450 G4Exception("G4RDVCrossSectionHandler::Bui
451 "InvalidCondition", FatalExcep
452
453 G4RDVDataSetAlgorithm* algo = CreateInterpol
454 G4RDVEMDataSet* materialSet = new G4RDCompos
455
456 G4DataVector* energies;
457 G4DataVector* data;
458
459 const G4ProductionCutsTable* theCoupleTable=
460 G4ProductionCutsTable::GetProductionCu
461 size_t numOfCouples = theCoupleTable->GetTab
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 = (*crossSecti
473 G4double materialCrossSection = 0.0;
474 G4int nElm = matCrossSet->NumberOfCo
475 for(G4int j=0; j<nElm; j++) {
476 materialCrossSection += matCrossSe
477 }
478
479 if (materialCrossSection > 0.)
480 {
481 data->push_back(1./materialCrossSectio
482 }
483 else
484 {
485 data->push_back(DBL_MAX);
486 }
487 }
488 G4RDVDataSetAlgorithm* algo = CreateInte
489 G4RDVEMDataSet* dataSet = new G4RDEMData
490 materialSet->AddComponent(dataSet);
491 }
492
493 return materialSet;
494 }
495
496 G4int G4RDVCrossSectionHandler::SelectRandomAt
497
498 {
499 // Select randomly an element within the mat
500 // determined by the cross sections in the d
501
502 const G4Material* material = couple->GetMate
503 G4int nElements = material->GetNumberOfEleme
504
505 // Special case: the material consists of on
506 if (nElements == 1)
507 {
508 G4int Z = (G4int) material->GetZ();
509 return Z;
510 }
511
512 // Composite material
513
514 const G4ElementVector* elementVector = mater
515 size_t materialIndex = couple->GetIndex();
516
517 G4RDVEMDataSet* materialSet = (*crossSection
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(
524 materialCrossSection0 += cr;
525 cross.push_back(materialCrossSection0);
526 }
527
528 G4double random = G4UniformRand() * material
529
530 for (G4int k=0 ; k < nElements ; k++ )
531 {
532 if (random <= cross[k]) return (G4int) (
533 }
534 // It should never get here
535 return 0;
536 }
537
538 const G4Element* G4RDVCrossSectionHandler::Sel
539 G4double e) const
540 {
541 // Select randomly an element within the mat
542 // by the cross sections in the data set
543
544 const G4Material* material = couple->GetMate
545 G4Element* nullElement = 0;
546 G4int nElements = material->GetNumberOfEleme
547 const G4ElementVector* elementVector = mater
548
549 // Special case: the material consists of on
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 = (*crossSec
562 G4double materialCrossSection0 = 0.0;
563 G4DataVector cross;
564 cross.clear();
565 for (G4int i=0; i<nElements; i++)
566 {
567 G4double cr = materialSet->GetCompon
568 materialCrossSection0 += cr;
569 cross.push_back(materialCrossSection
570 }
571
572 G4double random = G4UniformRand() * mate
573
574 for (G4int k=0 ; k < nElements ; k++ )
575 {
576 if (random <= cross[k]) return (*ele
577 }
578 // It should never end up here
579 G4cout << "G4RDVCrossSectionHandler::Sel
580 return nullElement;
581 }
582 }
583
584 G4int G4RDVCrossSectionHandler::SelectRandomSh
585 {
586 // Select randomly a shell, according to the
587 // in the data set
588
589 // Note for later improvement: it would be u
590 // used shells to improve performance
591
592 G4int shell = 0;
593
594 G4double totCrossSection = FindValue(Z,e);
595 G4double random = G4UniformRand() * totCross
596 G4double partialSum = 0.;
597
598 G4RDVEMDataSet* dataSet = 0;
599 std::map<G4int,G4RDVEMDataSet*,std::less<G4i
600 pos = dataMap.find(Z);
601 // The following is a workaround for STL Obj
602 // which does not support the standard and d
603 // the syntax pos->first or pos->second
604 // if (pos != dataMap.end()) dataSet = pos->
605 if (pos != dataMap.end()) dataSet = (*pos).s
606
607 size_t nShells = dataSet->NumberOfComponents
608 for (size_t i=0; i<nShells; i++)
609 {
610 const G4RDVEMDataSet* shellDataSet = dat
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 = G4Mat
625 if (materialTable == 0)
626 G4Exception("G4RDVCrossSectionHandler::Act
627 "InvalidSetup", FatalException
628
629 G4int nMaterials = G4Material::GetNumberOfMa
630
631 for (G4int m=0; m<nMaterials; m++)
632 {
633 const G4Material* material= (*materialTa
634 const G4ElementVector* elementVector = m
635 const G4int nElements = material->GetNum
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 &&
642 {
643 activeZ.push_back(Z);
644 }
645 }
646 }
647 }
648
649 G4RDVDataSetAlgorithm* G4RDVCrossSectionHandle
650 {
651 G4RDVDataSetAlgorithm* algorithm = new G4RDL
652 return algorithm;
653 }
654
655 G4int G4RDVCrossSectionHandler::NumberOfCompon
656 {
657 G4int n = 0;
658
659 std::map<G4int,G4RDVEMDataSet*,std::less<G4i
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: G4RDVCrossSectionHan
669 << "find Z = "
670 << Z << G4endl;
671 }
672 return n;
673 }
674
675
676