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

Last change on this file since 1212 was 1199, checked in by garnier, 16 years ago

nvx fichiers dans CVS
File size: 17.1 KB
RevLine 
[1199]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: materials-V09-02-18 $
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.