source: trunk/source/materials/test/G4SandiaTableTest.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
File size: 17.1 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: G4SandiaTableTest.cc,v 1.9 2006/06/29 19:13:17 gunter Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30//
31//
32////////////////////////////////////////////////////////////////////////
33//
34//  This program illustrates the different ways to define photoabsorption
35//  cross section according G4Sandiatable
36//
37// History:
38//
39// 15.09.99 V.Grichine, start from G4MaterialTest.cc
40//
41
42#include "G4ios.hh"
43#include <iomanip>
44#include "globals.hh"
45#include "G4UnitsTable.hh"
46
47#include "G4Material.hh"
48#include "G4SandiaTable.hh"
49
50int main() 
51{
52  // set output format
53
54  G4cout.setf( std::ios::scientific, std::ios::floatfield );
55
56  G4String name, symbol;             // a=mass of a mole;
57  G4double a, z, density;            // z=mean number of protons; 
58  G4int iz, n;                       // iz=nb of protons  in an isotope;
59                                     // n=nb of nucleons in an isotope;
60
61  G4int ncomponents, natoms, nel ;
62  G4double abundance, fractionmass ;
63  G4double temperature, pressure ;
64
65  G4UnitDefinition::BuildUnitsTable();
66
67//
68// define Elements
69//
70
71  a = 1.01*g/mole;
72  G4Element* elHold  = new G4Element(name="Hydrogen",symbol="H" , z= 1., a);
73
74  a = 1.01*g/mole;
75  G4Isotope* ih1 = new G4Isotope("Hydrogen",iz=1,n=1,a);
76
77  a = 2.01*g/mole;
78  G4Isotope* ih2 = new G4Isotope("Deuterium",iz=1,n=2,a);
79
80  G4Element* elH = new G4Element(name="Hydrogen",symbol="H",2);
81  elH->AddIsotope(ih1,.999);
82  elH->AddIsotope(ih2,.001);
83
84
85  a = 12.01*g/mole;
86  G4Element* elC  = new G4Element(name="Carbon"  ,symbol="C" , z= 6., a);
87
88  a = 14.01*g/mole;
89  G4Element* elN  = new G4Element(name="Nitrogen",symbol="N" , z= 7., a);
90
91  a = 16.00*g/mole;
92  G4Element* elO  = new G4Element(name="Oxygen"  ,symbol="O" , z= 8., a);
93
94  a = 28.09*g/mole;
95  G4Element* elSi = new G4Element(name="Silicon",symbol="Si" , z= 14., a);
96
97  a = 55.85*g/mole;
98  G4Element* elFe = new G4Element(name="Iron"    ,symbol="Fe", z=26., a);
99
100  a = 131.29*g/mole;
101  G4Element* elXe = new G4Element(name="Xenon", symbol="Xe", z=54., a);
102 
103  a = 39.948*g/mole;
104  G4Element* elAr = new G4Element(name="Argon", symbol="Ar", z=18., a);
105
106
107  a = 19.00*g/mole;
108  G4Element* elF  = new G4Element(name="Fluorine", symbol="F", z=9., a);
109
110  a = 69.723*g/mole;
111  G4Element* elGa  = new G4Element(name="Ga", symbol="Ga", z=31., a);
112
113  a = 74.9216*g/mole;
114  G4Element* elAs  = new G4Element(name="As", symbol="As", z=33., a);
115
116
117//
118// define an Element from isotopes, by relative abundance
119//
120
121  G4Isotope* U5 = new G4Isotope(name="U235", iz=92, n=235, a=235.01*g/mole);
122  G4Isotope* U8 = new G4Isotope(name="U238", iz=92, n=238, a=238.03*g/mole);
123
124  G4Element* elU  = new G4Element(name="enriched Uranium", 
125                                symbol="U", ncomponents=2);
126  elU->AddIsotope(U5, abundance= 90.*perCent);
127  elU->AddIsotope(U8, abundance= 10.*perCent);
128
129
130// G4cout << *(G4Isotope::GetIsotopeTable()) << G4endl;
131
132// G4cout << *(G4Element::GetElementTable()) << G4endl;
133
134//
135// define simple materials
136//
137
138  density = 2.700*g/cm3;
139  a = 26.98*g/mole;
140  G4Material* Al = new G4Material(name="Aluminium", z=13., a, density);
141
142  density = 1.390*g/cm3;
143  a = 39.95*g/mole;
144  G4Material* lAr = new G4Material(name="liquidArgon", z=18., a, density);
145
146  density = 8.960*g/cm3;
147  a = 63.55*g/mole;
148  G4Material* Cu = new G4Material(name="Copper"   , z=29., a, density);
149
150  density = 11.35*g/cm3;
151  a = 207.19*g/mole;
152  G4Material* Pb = new G4Material(name="Lead "     , z=82., a, density);
153
154//
155// define a material from elements.   case 1: chemical molecule
156//
157 
158  density = 1.000*g/cm3;
159  G4Material* H2O = new G4Material(name="Water", density, ncomponents=2);
160  H2O->AddElement(elH, natoms=2);
161  H2O->AddElement(elO, natoms=1);
162
163  density = 1.032*g/cm3;
164  G4Material* Sci = new G4Material(name="Scintillator", density, ncomponents=2);
165  Sci->AddElement(elC, natoms=9);
166  Sci->AddElement(elH, natoms=10);
167
168  density = 2.200*g/cm3;
169  G4Material* SiO2 = new G4Material(name="quartz", density, ncomponents=2);
170  SiO2->AddElement(elSi, natoms=1);
171  SiO2->AddElement(elO , natoms=2);
172
173//
174// define a material from elements.   case 2: mixture by fractional mass
175//
176
177  density = 1.290*mg/cm3;
178  G4Material* oldAir = new G4Material(name="Air  "  , density, ncomponents=2);
179  oldAir->AddElement(elN, fractionmass=0.7);
180  oldAir->AddElement(elO, fractionmass=0.3);
181
182//
183// define a material from elements and/or others materials (mixture of mixtures)
184//
185
186  density = 0.200*g/cm3;
187  G4Material* Aerog = new G4Material(name="Aerogel", density, ncomponents=3);
188  Aerog->AddMaterial(SiO2, fractionmass=0.625);
189  Aerog->AddMaterial(H2O , fractionmass=0.374);
190  Aerog->AddElement (elC , fractionmass=0.1*perCent);
191
192//
193// examples of gas in non STP conditions
194//
195
196  density     = 27.*mg/cm3;
197  pressure    = 50.*atmosphere;
198  temperature = 325.*kelvin;
199  G4Material* CO2 = new G4Material(name="Carbonic gas", density, ncomponents=2,
200                                     kStateGas,temperature,pressure);
201  CO2->AddElement(elC, natoms=1);
202  CO2->AddElement(elO, natoms=2);
203 
204  density     = 0.3*mg/cm3;
205  pressure    = 2.*atmosphere;
206  temperature = 500.*kelvin;
207  G4Material* steam = new G4Material(name="Water steam ", density, ncomponents=1,
208                                      kStateGas,temperature,pressure);
209  steam->AddMaterial(H2O, fractionmass=1.);
210
211//
212// examples of vacuum
213//
214
215  density     = universe_mean_density;            //from PhysicalConstants.h
216  pressure    = 3.e-18*pascal;
217  temperature = 2.73*kelvin;
218  new G4Material(name="Galactic", z=1., a=1.01*g/mole, density,
219                   kStateGas,temperature,pressure);
220
221  density     = 1.e-5*g/cm3;
222  pressure    = 2.e-2*bar;
223  temperature = STP_Temperature;                      //from PhysicalConstants.h
224
225  G4Material* beam = new G4Material(name="Beam ", density, ncomponents=1,
226                                      kStateGas,temperature,pressure);
227  beam->AddMaterial(oldAir, fractionmass=1.);
228
229  // maylar
230
231  density = 1.39*g/cm3;
232  G4Material* Maylar = new G4Material(name="Maylar", density, nel=3);
233  Maylar->AddElement(elO,2);
234  Maylar->AddElement(elC,5);
235  Maylar->AddElement(elH,4);
236
237  // Kapton Dupont de Nemur (density: 1.396-1.430, get middle )
238
239  density = 1.413*g/cm3;
240  G4Material* Kapton = new G4Material(name="Kapton", density, nel=4);
241  Kapton->AddElement(elO,5);
242  Kapton->AddElement(elC,22);
243  Kapton->AddElement(elN,2);
244  Kapton->AddElement(elH,10);
245
246  // Germanium as detector material
247
248  density = 5.323*g/cm3;
249  a = 72.59*g/mole;
250  G4Material* Ge = new G4Material(name="Ge", z=32., a, density);
251
252  // GaAs detectors
253
254  density = 5.32*g/cm3;
255  G4Material* GaAs = new G4Material(name="GaAs",density, nel=2);
256  GaAs->AddElement(elGa,1);
257  GaAs->AddElement(elAs,1);
258
259  // Diamond detectors
260
261  density = 3.5*g/cm3;
262  G4Material* Diamond = new G4Material(name="Diamond",density, nel=1);
263  Diamond->AddElement(elC,1);
264
265
266  G4double TRT_Xe_density = 5.485*mg/cm3;
267  G4Material* TRT_Xe = new G4Material(name="TRT_Xe", TRT_Xe_density, nel=1,
268                                      kStateGas,293.15*kelvin,1.*atmosphere);
269  TRT_Xe->AddElement(elXe,1);
270
271  G4double TRT_CO2_density = 1.842*mg/cm3;
272  G4Material* TRT_CO2 = new G4Material(name="TRT_CO2", TRT_CO2_density, nel=2,
273                                       kStateGas,293.15*kelvin,1.*atmosphere);
274  TRT_CO2->AddElement(elC,1);
275  TRT_CO2->AddElement(elO,2);
276
277  G4double TRT_CF4_density = 3.9*mg/cm3;
278  G4Material* TRT_CF4 = new G4Material(name="TRT_CF4", TRT_CF4_density, nel=2,
279                                       kStateGas,293.15*kelvin,1.*atmosphere);
280  TRT_CF4->AddElement(elC,1);
281  TRT_CF4->AddElement(elF,4);
282
283  G4double XeCO2CF4_density = 4.76*mg/cm3;
284  G4Material* XeCO2CF4 = new G4Material(name="XeCO2CF4", XeCO2CF4_density,
285                                        ncomponents=3,
286                                        kStateGas,293.15*kelvin,1.*atmosphere);
287  XeCO2CF4->AddMaterial(TRT_Xe,0.807);
288  XeCO2CF4->AddMaterial(TRT_CO2,0.039);
289  XeCO2CF4->AddMaterial(TRT_CF4,0.154);
290     
291  density = 0.935*g/cm3;
292  G4Material* TRT_CH2 = new G4Material(name="TRT_CH2",density, nel=2);
293  TRT_CH2->AddElement(elC,1);
294  TRT_CH2->AddElement(elH,2);
295
296  density = 0.059*g/cm3;
297  G4Material* Radiator = new G4Material(name="Radiator",density, nel=2);
298  Radiator->AddElement(elC,1);
299  Radiator->AddElement(elH,2);
300
301  density = 0.145*g/cm3;
302  G4Material* CarbonFiber = new G4Material(name="CarbonFiber",density, nel=1);
303  CarbonFiber->AddElement(elC,1);
304
305  // Dry air (average composition)
306
307
308  density = 1.25053*mg/cm3 ;       // STP
309  G4Material* Nitrogen = new G4Material(name="N2"  , density, ncomponents=1);
310  Nitrogen->AddElement(elN, 2);
311
312  density = 1.4289*mg/cm3 ;       // STP
313  G4Material* Oxygen = new G4Material(name="O2"  , density, ncomponents=1);
314  Oxygen->AddElement(elO, 2);
315
316  density = 1.7836*mg/cm3 ;       // STP
317  G4Material* Argon = new G4Material(name="Argon"  , density, ncomponents=1);
318  Argon->AddElement(elAr, 1);
319
320  density = 1.2928*mg/cm3 ;       // STP
321  G4Material* Air = new G4Material(name="Air"  , density, ncomponents=3);
322  Air->AddMaterial( Nitrogen, fractionmass = 0.7557 ) ;
323  Air->AddMaterial( Oxygen,   fractionmass = 0.2315 ) ;
324  Air->AddMaterial( Argon,    fractionmass = 0.0128 ) ;
325
326  // Xenon as detector gas, STP
327
328  density = 5.858*mg/cm3 ;
329  a = 131.29*g/mole ;
330  G4Material* Xe  = new G4Material(name="Xenon",z=54., a, density );
331
332  // Helium as detector gas, STP
333
334  density = 0.178*mg/cm3 ;
335  a = 4.0026*g/mole ;
336  G4Material* He  = new G4Material(name="He",z=2., a, density );
337
338  // Neon as detector gas, STP
339
340  density = 0.900*mg/cm3 ;
341  a = 20.179*g/mole ;
342  G4Material* Ne  = new G4Material(name="Ne",z=10., a, density );
343
344  // Krypton as detector gas, STP
345
346  density = 3.700*mg/cm3 ;
347  a = 83.80*g/mole ;
348  G4Material* Kr  = new G4Material(name="Kr",z=36., a, density );
349
350  // Carbone dioxide, CO2 STP
351
352  density = 1.977*mg/cm3 ;
353  G4Material* CarbonDioxide = new G4Material(name="CO2", density, nel=2) ;
354  CarbonDioxide->AddElement(elC,1) ;
355  CarbonDioxide->AddElement(elO,2) ;
356
357  // Metane, STP
358
359  density = 0.7174*mg/cm3 ;
360  G4Material* metane = new G4Material(name="CH4",density,nel=2) ;
361  metane->AddElement(elC,1) ;
362  metane->AddElement(elH,4) ;
363
364  // Propane, STP
365
366  density = 2.005*mg/cm3 ;
367  G4Material* propane = new G4Material(name="C3H8",density,nel=2) ;
368  propane->AddElement(elC,3) ;
369  propane->AddElement(elH,8) ;
370
371  // iso-Butane (methylpropane), STP
372
373  density = 2.67*mg/cm3 ;
374  G4Material* isobutane = new G4Material(name="isoC4H10",density,nel=2) ;
375  isobutane->AddElement(elC,4) ;
376  isobutane->AddElement(elH,10) ;
377
378// 87.5% Xe + 7.5% CH4 + 5% C3H8, 20 C, 1 atm
379
380  density = 4.9196*mg/cm3 ;
381
382  G4Material* XeCH4C3H8 = new G4Material(name="XeCH4C3H8"  , density,
383                                                             ncomponents=3);
384  XeCH4C3H8->AddMaterial( Xe,       fractionmass = 0.971 ) ;
385  XeCH4C3H8->AddMaterial( metane,   fractionmass = 0.010 ) ;
386  XeCH4C3H8->AddMaterial( propane,  fractionmass = 0.019 ) ;
387
388  // Propane in MWPC, 2 atm, 20 C
389
390  //  density = 3.758*mg/cm3 ;
391  density = 3.736*mg/cm3 ;
392  G4Material* propaneDet = new G4Material(name="detC3H8",density,nel=2) ;
393  propaneDet->AddElement(elC,3) ;
394  propaneDet->AddElement(elH,8) ;
395
396  // 80% Ar + 20% CO2, STP
397
398  density = 1.8223*mg/cm3 ;
399  G4Material* Ar20CO2 = new G4Material(name="Ar20CO2"  , density,
400                                                             ncomponents=2);
401  Ar20CO2->AddMaterial( Argon,           fractionmass = 0.783 ) ;
402  Ar20CO2->AddMaterial( CarbonDioxide,   fractionmass = 0.217 ) ;
403
404  // 93% Ar + 7% CH4, STP
405
406  density = 1.709*mg/cm3 ;
407  G4Material* Ar7CH4 = new G4Material(name="Ar7CH4"  , density,
408                                                             ncomponents=2);
409  Ar7CH4->AddMaterial( Argon,    fractionmass = 0.971 ) ;
410  Ar7CH4->AddMaterial( metane,   fractionmass = 0.029 ) ;
411
412  // 80% Xe + 20% CO2, STP
413
414  density = 5.0818*mg/cm3 ;
415  G4Material* Xe20CO2 = new G4Material(name="Xe20CO2"  , density,
416                                                             ncomponents=2);
417  Xe20CO2->AddMaterial( Xe,              fractionmass = 0.922 ) ;
418  Xe20CO2->AddMaterial( CarbonDioxide,   fractionmass = 0.078 ) ;
419
420  // 80% Kr + 20% CO2, STP
421
422  density = 3.601*mg/cm3 ;
423  G4Material* Kr20CO2 = new G4Material(name="Kr20CO2"  , density,
424                                                             ncomponents=2);
425  Kr20CO2->AddMaterial( Kr,              fractionmass = 0.89 ) ;
426  Kr20CO2->AddMaterial( CarbonDioxide,   fractionmass = 0.11 ) ;
427
428  // 80% He + 20% CO2, STP
429
430  density = 0.5378*mg/cm3 ;
431  G4Material* He20CO2 = new G4Material(name="He20CO2"  , density,
432                                                             ncomponents=2);
433  He20CO2->AddMaterial( He,              fractionmass = 0.265 ) ;
434  He20CO2->AddMaterial( CarbonDioxide,   fractionmass = 0.735 ) ;
435
436
437
438//
439// Print the table of materials
440//
441
442// G4cout << *(G4Material::GetMaterialTable()) << G4endl;
443
444//
445////////////////////////////////////////////////////////////////////////
446//
447//
448// Checking Sandia table coefficients
449//
450  G4int numberOfMat, iMat, matIndex, nbOfElements, sanIndex, row, iSan;
451  G4double unit;
452  G4String materialName = "Air";
453  static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
454  numberOfMat = theMaterialTable->size() ;
455
456  for(iMat=0;iMat<numberOfMat;iMat++)
457  {
458    if(materialName == (*theMaterialTable)[iMat]->GetName() )
459    {
460      matIndex = (*theMaterialTable)[iMat]->GetIndex() ;
461      break ;
462    }
463  }
464 
465//
466////////////////////////////////////////////////////////////////////////
467//
468// Sandia cof according old PAI stuff
469//
470  for(iMat=0;iMat<numberOfMat;iMat++)
471  {
472     G4String matName = (*theMaterialTable)[iMat]->GetName();
473     matIndex = (*theMaterialTable)[iMat]->GetIndex();
474     nbOfElements = (*theMaterialTable)[iMat]->GetNumberOfElements();
475     density = (*theMaterialTable)[iMat]->GetDensity();
476     
477     G4cout<<matIndex<<"\t"<<matName<<G4endl<<G4endl;
478     
479     G4cout<<"Sandia cof according old PAI stuff"<<G4endl<<G4endl;
480
481     G4int* thisMaterialZ = new G4int[nbOfElements];
482     for(iSan=0;iSan<nbOfElements;iSan++)
483     {
484        thisMaterialZ[iSan] = (G4int)(*theMaterialTable)[iMat]->
485                                      GetElement(iSan)->GetZ();
486     }     
487     G4SandiaTable sandia(matIndex) ;
488
489     //  sanIndex = sandia.SandiaIntervals(thisMaterialZ,nbOfElements);   
490     //  sanIndex = sandia.SandiaMixing( thisMaterialZ ,
491     //                             (*theMaterialTable)[iMat]->GetFractionVector() ,
492     //                      nbOfElements,sanIndex) ;
493     sanIndex = sandia.GetMaxInterval() ;
494     G4cout<<"fMaxInterval = "<<sanIndex<<G4endl<<G4endl;
495
496     for(row = 0; row < sanIndex - 1 ; row++)
497     {
498       G4cout<<row+1<<"\t"<<sandia.GetPhotoAbsorpCof(row+1,0)/keV;
499       
500       unit = cm2/g;
501       for(iSan = 1; iSan < 5; iSan++)
502       {
503         unit *= keV;         
504         G4cout<<"\t"<<sandia.GetPhotoAbsorpCof(row+1,iSan)/unit;
505       }
506       G4cout<<G4endl ;
507     }
508     G4cout<<G4endl ;
509     
510//
511////////////////////////////////////////////////////////////////////////
512//
513// Sandia cof according ComputeMatSandiaMatrix()
514//
515   G4SandiaTable* sanMatrix = G4Material::GetMaterial(matName)->
516                              GetSandiaTable();
517   sanIndex = sanMatrix->GetMatNbOfIntervals();
518     
519   G4cout<<"Sandia cof according ComputeMatSandiaMatrix()"<<G4endl<<G4endl;
520
521     for (row=0; row<sanIndex; row++) {
522       G4cout<<row+1<<"\t"<<sanMatrix->GetSandiaCofForMaterial(row,0)/keV;
523               
524       unit = cm2/g;
525       for (iSan=1; iSan<5; iSan++) {
526         unit *= keV; 
527         G4cout<<"\t"<<(sanMatrix->GetSandiaCofForMaterial(row,iSan))
528                                                            /(density*unit);
529       }
530       G4cout<<G4endl;
531     }     
532     G4cout<<G4endl;     
533  }
534  return EXIT_SUCCESS;
535}
536
537//
538//
539///////////////////////  end of G4SandiaTableTest.cc  /////////////////////////
Note: See TracBrowser for help on using the repository browser.