source: trunk/source/processes/electromagnetic/lowenergy/src/G4VCrossSectionHandler.cc@ 1058

Last change on this file since 1058 was 1007, checked in by garnier, 17 years ago

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