Context Navigation

testLowEnergyPAI.cc

Last change on this file was 1199, checked in by garnier, 15 years ago
nvx fichiers dans CVS
File size: 14.6 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	//
28	//
29	//
30	// Test routine based on G4PAIdNdxTest for low energy extension of PAI model
31	//
32	// History:
33	//
34	// 18.09.09, V. Grichine implementation based on G4PAIdNdxTest
35
36	#include "G4ios.hh"
37	#include <fstream>
38	#include <cmath>
39	#include "globals.hh"
40	#include "Randomize.hh"
41
42	#include "G4Isotope.hh"
43	#include "G4Element.hh"
44	#include "G4Material.hh"
45	#include "G4MaterialTable.hh"
46	#include "G4SandiaTable.hh"
47
48	// #include "G4PAIonisation.hh"
49	#include "G4PAIxSection.hh"
50
51
52
53	G4double GetXscSimple(G4Element* el, G4double Tkin)
54	{
55	G4int shell, NbOfShells, Nshell;
56	G4double tmpxsc, lnt, t, B, xsc = 0.;
57	G4double minust, plust, tplus, pert;
58
59	G4double ryd = 13.6*eV; // Rydberg number
60	G4double ryd2 = ryd*ryd;
61
62	G4double cof = 4piBohr_radiusBohr_radiusryd2;
63
64	NbOfShells = el->GetNbOfAtomicShells();
65
66
67	for ( shell = 0; shell < NbOfShells; shell++ )
68	{
69	B = el->GetAtomicShell(shell);
70	Nshell = el->GetNbOfShellElectrons(shell);
71
72	t = Tkin/B;
73	lnt = std::log(t);
74	tplus = t + 1.;
75	pert = 1./t;
76	minust = 1. - pert;
77	plust = 1. + pert;
78
79	tmpxsc = 0.5minustplust*lnt + minust - lnt/tplus;
80	tmpxsc *= Nshell/B/B/(tplus+1);
81	xsc += tmpxsc;
82	}
83	xsc *= cof;
84	return xsc;
85	}
86
87
88
89
90
91
92
93
94	int main()
95	{
96	std::ofstream outFile("PAIdNdx.out", std::ios::out ) ;
97	outFile.setf( std::ios::scientific, std::ios::floatfield );
98
99	std::ofstream fileOut("PAICrenkovPlasmon.out", std::ios::out ) ;
100	fileOut.setf( std::ios::scientific, std::ios::floatfield );
101
102	// std::ifstream fileRead("exp.dat", std::ios::out ) ;
103	// fileRead.setf( std::ios::scientific, std::ios::floatfield );
104
105	std::ofstream fileWrite("exp.dat", std::ios::out ) ;
106	fileWrite.setf( std::ios::scientific, std::ios::floatfield );
107
108	std::ofstream fileWrite1("mprrpai.dat", std::ios::out ) ;
109	fileWrite1.setf( std::ios::scientific, std::ios::floatfield );
110
111	// Create materials
112
113
114	G4int iz , n, nel, ncomponents ;
115	G4double a, z, ez, density , temperature, pressure, fractionmass ;
116	G4State state ;
117	G4String name, symbol ;
118
119	// G4Element* elH = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
120
121	a = 14.01*g/mole;
122	G4Element* elN = new G4Element(name="Nitrogen", symbol="N", ez=7., a);
123
124	a = 16.00*g/mole;
125	// G4Element* elO = new G4Element(name="Oxigen", symbol="O", ez=8., a);
126
127	a = 12.01*g/mole;
128	G4Element* elC = new G4Element(name="Carbon",symbol="C", ez=6., a);
129
130	a = 55.85*g/mole;
131	G4Element* elFe = new G4Element(name="Iron",symbol="Fe", ez=26., a);
132
133	a = 16.00*g/mole;
134	G4Element* elO = new G4Element(name="Oxygen",symbol="O", ez=8., a);
135
136	a = 1.01*g/mole;
137	G4Isotope* ih1 = new G4Isotope("Hydrogen",iz=1,n=1,a);
138
139	a = 2.01*g/mole;
140	G4Isotope* ih2 = new G4Isotope("Deuterium",iz=1,n=2,a);
141
142	G4Element* elH = new G4Element(name="Hydrogen",symbol="H",2);
143	elH->AddIsotope(ih1,.999);
144	elH->AddIsotope(ih2,.001);
145
146	a = 39.948*g/mole;
147	G4Element* elAr = new G4Element(name="Argon", symbol="Ar", z=18., a);
148
149	a = 131.29*g/mole;
150	G4Element* elXe = new G4Element(name="Xenon", symbol="Xe", z=54., a);
151
152	a = 19.00*g/mole;
153	G4Element* elF = new G4Element(name="Fluorine", symbol="F", z=9., a);
154
155	a = 69.723*g/mole;
156	G4Element* elGa = new G4Element(name="Ga", symbol="Ga", z=31., a);
157
158	a = 74.9216*g/mole;
159	G4Element* elAs = new G4Element(name="As", symbol="As", z=33., a);
160
161
162	// G4Isotope::DumpInfo();
163	// G4Element::DumpInfo();
164	// G4Material::DumpInfo();
165
166	/* ***************************************************************
167
168	a = 9.012*g/mole;
169	density = 1.848*g/cm3;
170	G4Material* Be = new G4Material(name="Beryllium", z=4. , a, density);
171
172	density = 1.390*g/cm3;
173	a = 39.95*g/mole;
174	G4Material* lAr = new G4Material(name="liquidArgon", z=18., a, density);
175
176	density = 19.32*g/cm3;
177	a =196.97*g/mole;
178	G4Material* Au = new G4Material(name="Gold" , z=79., a, density);
179
180	// Carbon dioxide
181
182	density = 1.977*mg/cm3;
183	G4Material* CO2 = new G4Material(name="CO2", density, nel=2,
184	kStateGas,273.15kelvin,1.atmosphere);
185	CO2->AddElement(elC,1);
186	CO2->AddElement(elO,2);
187
188	density = 1.290*mg/cm3; // old air from elements
189	G4Material* air = new G4Material(name="air" , density, ncomponents=2);
190	Air->AddElement(elN, fractionmass=0.7);
191	Air->AddElement(elO, fractionmass=0.3);
192
193
194	density = 1.25053*mg/cm3 ; // STP
195	a = 14.01*g/mole ; // get atomic weight !!!
196	// a = 28.016*g/mole;
197	G4Material* newN2 = new G4Material(name="newN2", z= 7.,a,density) ;
198
199	density = 1.25053*mg/cm3 ; // STP
200	G4Material* anotherN2 = new G4Material(name="anotherN2", density,ncomponents=2);
201	anotherN2->AddElement(elN, 1);
202	anotherN2->AddElement(elN, 1);
203
204	density = 1.000*g/cm3;
205	G4Material* H2O = new G4Material(name="Water", density, ncomponents=2);
206	H2O->AddElement(elH, natoms=2);
207	H2O->AddElement(elO, natoms=1);
208
209
210
211	a = 26.98159*g/mole;
212	density = 2.7*g/cm3;
213	G4Material* Al = new G4Material(name="Aluminium", z=13., a, density);
214
215	// Silicon as detector material
216
217	density = 2.330*g/cm3;
218	a = 28.0855*g/mole;
219	G4Material* Si = new G4Material(name="Silicon", z=14., a, density);
220
221
222	density = 7.870*g/cm3;
223	a = 55.85*g/mole;
224	G4Material* Fe = new G4Material(name="Iron" , z=26., a, density);
225
226	density = 8.960*g/cm3;
227	a = 63.55*g/mole;
228	G4Material* Cu = new G4Material(name="Copper" , z=29., a, density);
229
230	density = 11.35*g/cm3;
231	a = 207.19*g/mole;
232	G4Material* Pb = new G4Material(name="Lead" , z=82., a, density);
233
234	// Polypropelene
235
236	G4Material* CH2 = new G4Material ("Polypropelene" , 0.91*g/cm3, 2);
237	CH2->AddElement(elH,2);
238	CH2->AddElement(elC,1);
239
240	// Kapton (polyimide)
241
242	density = 1.39*g/cm3;
243	G4Material* Kapton = new G4Material(name="Kapton", density, nel=3);
244	Kapton->AddElement(elO,2);
245	Kapton->AddElement(elC,5);
246	Kapton->AddElement(elH,4);
247
248
249	// Germanium as detector material
250
251	density = 5.323*g/cm3;
252	a = 72.59*g/mole;
253	G4Material* Ge = new G4Material(name="Ge", z=32., a, density);
254
255	// GaAs detectors
256
257	density = 5.32*g/cm3;
258	G4Material* GaAs = new G4Material(name="GaAs",density, nel=2);
259	GaAs->AddElement(elGa,1);
260	GaAs->AddElement(elAs,1);
261
262	// Diamond detectors
263
264	density = 3.5*g/cm3;
265	G4Material* Diamond = new G4Material(name="Diamond",density, nel=1);
266	Diamond->AddElement(elC,1);
267
268	G4double TRT_Xe_density = 5.485*mg/cm3;
269	G4Material* TRT_Xe = new G4Material(name="TRT_Xe", TRT_Xe_density, nel=1,
270	kStateGas,293.15kelvin,1.atmosphere);
271	TRT_Xe->AddElement(elXe,1);
272
273	G4double TRT_CO2_density = 1.842*mg/cm3;
274	G4Material* TRT_CO2 = new G4Material(name="TRT_CO2", TRT_CO2_density, nel=2,
275	kStateGas,293.15kelvin,1.atmosphere);
276	TRT_CO2->AddElement(elC,1);
277	TRT_CO2->AddElement(elO,2);
278
279	G4double TRT_CF4_density = 3.9*mg/cm3;
280	G4Material* TRT_CF4 = new G4Material(name="TRT_CF4", TRT_CF4_density, nel=2,
281	kStateGas,293.15kelvin,1.atmosphere);
282	TRT_CF4->AddElement(elC,1);
283	TRT_CF4->AddElement(elF,4);
284
285	// ATLAS TRT straw tube gas mixture (20 C, 1 atm)
286
287	G4double XeCO2CF4_density = 4.76*mg/cm3;
288	G4Material* XeCO2CF4 = new G4Material(name="XeCO2CF4", XeCO2CF4_density,
289	ncomponents=3,
290	kStateGas,293.15kelvin,1.atmosphere);
291	XeCO2CF4->AddMaterial(TRT_Xe,0.807);
292	XeCO2CF4->AddMaterial(TRT_CO2,0.039);
293	XeCO2CF4->AddMaterial(TRT_CF4,0.154);
294
295	// TRT_CH2
296
297	density = 0.935*g/cm3;
298	G4Material* TRT_CH2 = new G4Material(name="TRT_CH2",density, nel=2);
299	TRT_CH2->AddElement(elC,1);
300	TRT_CH2->AddElement(elH,2);
301
302	// Radiator
303
304	density = 0.059*g/cm3;
305	G4Material* Radiator = new G4Material(name="Radiator",density, nel=2);
306	Radiator->AddElement(elC,1);
307	Radiator->AddElement(elH,2);
308
309	// Carbon Fiber
310
311	density = 0.145*g/cm3;
312	G4Material* CarbonFiber = new G4Material(name="CarbonFiber",density, nel=1);
313	CarbonFiber->AddElement(elC,1);
314
315
316	// Dry air (average composition)
317
318
319	density = 1.25053*mg/cm3 ; // STP
320	G4Material* Nitrogen = new G4Material(name="N2" , density, ncomponents=1);
321	Nitrogen->AddElement(elN, 2);
322
323	density = 1.4289*mg/cm3 ; // STP
324	G4Material* Oxygen = new G4Material(name="O2" , density, ncomponents=1);
325	Oxygen->AddElement(elO, 2);
326
327	density = 1.7836*mg/cm3 ; // STP
328	G4Material* Argon = new G4Material(name="Argon" , density, ncomponents=1);
329	Argon->AddElement(elAr, 1);
330
331	density = 1.2928*mg/cm3 ; // STP
332	G4Material* Air = new G4Material(name="Air" , density, ncomponents=3);
333	Air->AddMaterial( Nitrogen, fractionmass = 0.7557 ) ;
334	Air->AddMaterial( Oxygen, fractionmass = 0.2315 ) ;
335	Air->AddMaterial( Argon, fractionmass = 0.0128 ) ;
336
337	// Xenon as detector gas, STP
338
339	density = 5.858*mg/cm3 ;
340	a = 131.29*g/mole ;
341	G4Material* Xe = new G4Material(name="Xenon",z=54., a, density );
342
343	// Helium as detector gas, STP
344
345	density = 0.178*mg/cm3 ;
346	a = 4.0026*g/mole ;
347	G4Material* He = new G4Material(name="He",z=2., a, density );
348
349
350	// Krypton as detector gas, STP
351
352	density = 3.700*mg/cm3 ;
353	a = 83.80*g/mole ;
354	G4Material* Kr = new G4Material(name="Kr",z=36., a, density );
355
356	****************************************************** */
357
358	// Neon as detector gas, STP
359
360	density = 0.900*mg/cm3 ;
361	a = 20.179*g/mole ;
362	G4Material* Ne = new G4Material(name="Ne",z=10., a, density );
363
364	// Carbone dioxide, CO2 STP
365
366	density = 1.977*mg/cm3 ;
367	G4Material* CarbonDioxide = new G4Material(name="CO2", density, nel=2) ;
368	CarbonDioxide->AddElement(elC,1) ;
369	CarbonDioxide->AddElement(elO,2) ;
370
371	/* *****************************************************
372
373	// Metane, STP
374
375	density = 0.7174*mg/cm3 ;
376	G4Material* metane = new G4Material(name="CH4",density,nel=2) ;
377	metane->AddElement(elC,1) ;
378	metane->AddElement(elH,4) ;
379
380	// Propane, STP
381
382	density = 2.005*mg/cm3 ;
383	G4Material* propane = new G4Material(name="C3H8",density,nel=2) ;
384	propane->AddElement(elC,3) ;
385	propane->AddElement(elH,8) ;
386
387	// iso-Butane (methylpropane), STP
388
389	density = 2.67*mg/cm3 ;
390	G4Material* isobutane = new G4Material(name="isoC4H10",density,nel=2) ;
391	isobutane->AddElement(elC,4) ;
392	isobutane->AddElement(elH,10) ;
393
394	// 87.5% Xe + 7.5% CH4 + 5% C3H8, 20 C, 1 atm
395
396	density = 4.9196*mg/cm3 ;
397
398	G4Material* XeCH4C3H8 = new G4Material(name="XeCH4C3H8" , density,
399	ncomponents=3);
400	XeCH4C3H8->AddMaterial( Xe, fractionmass = 0.971 ) ;
401	XeCH4C3H8->AddMaterial( metane, fractionmass = 0.010 ) ;
402	XeCH4C3H8->AddMaterial( propane, fractionmass = 0.019 ) ;
403
404	// Propane in MWPC, 2 atm, 20 C
405
406	// density = 3.758*mg/cm3 ;
407	density = 3.736*mg/cm3 ;
408	G4Material* propaneDet = new G4Material(name="detC3H8",density,nel=2) ;
409	propaneDet->AddElement(elC,3) ;
410	propaneDet->AddElement(elH,8) ;
411
412	************************************************** */
413
414	// 90% Ne + 10% CO2, STP
415
416	density = 1.0077*mg/cm3 ;
417	G4Material* Ne10CO2 = new G4Material(name="Ne10CO2" , density,
418
419	ncomponents=2);
420	Ne10CO2->AddMaterial( Ne, fractionmass = 0.8038 ) ;
421	Ne10CO2->AddMaterial( CarbonDioxide, fractionmass = 0.1962 ) ;
422
423	/* *****************************************************
424
425	// 80% Ar + 20% CO2, STP
426
427	density = 1.8223*mg/cm3 ;
428	G4Material* Ar20CO2 = new G4Material(name="Ar20CO2" , density,
429	ncomponents=2);
430	Ar20CO2->AddMaterial( Argon, fractionmass = 0.783 ) ;
431	Ar20CO2->AddMaterial( CarbonDioxide, fractionmass = 0.217 ) ;
432
433	// 93% Ar + 7% CH4, STP
434
435	density = 1.709*mg/cm3 ;
436	G4Material* Ar7CH4 = new G4Material(name="Ar7CH4" , density,
437	ncomponents=2);
438	Ar7CH4->AddMaterial( Argon, fractionmass = 0.971 ) ;
439	Ar7CH4->AddMaterial( metane, fractionmass = 0.029 ) ;
440
441	// 80% Xe + 20% CO2, STP
442
443	density = 5.0818*mg/cm3 ;
444	G4Material* Xe20CO2 = new G4Material(name="Xe20CO2" , density,
445	ncomponents=2);
446	Xe20CO2->AddMaterial( Xe, fractionmass = 0.922 ) ;
447	Xe20CO2->AddMaterial( CarbonDioxide, fractionmass = 0.078 ) ;
448
449	// 80% Kr + 20% CO2, STP
450
451	density = 3.601*mg/cm3 ;
452	G4Material* Kr20CO2 = new G4Material(name="Kr20CO2" , density,
453	ncomponents=2);
454	Kr20CO2->AddMaterial( Kr, fractionmass = 0.89 ) ;
455	Kr20CO2->AddMaterial( CarbonDioxide, fractionmass = 0.11 ) ;
456
457	// 80% He + 20% CO2, STP
458
459	density = 0.5378*mg/cm3 ;
460	G4Material* He20CO2 = new G4Material(name="He20CO2" , density,
461	ncomponents=2);
462	He20CO2->AddMaterial( He, fractionmass = 0.265 ) ;
463	He20CO2->AddMaterial( CarbonDioxide, fractionmass = 0.735 ) ;
464
465	*/ //////////////////////
466
467
468	// G4cout << *(G4Material::GetMaterialTable()) << G4endl;
469
470
471	// G4int i, j, k, numOfMaterials, iSan, nbOfElements, sanIndex, row ;
472	// G4double maxEnergyTransfer, kineticEnergy ;
473	// G4double tau, gamma, bg2, beta2, rateMass, Tmax, Tmin, Tkin ;
474
475	const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable() ;
476
477
478	return EXIT_SUCCESS;
479
480	}
481
482
483
484
485
486
487
488
489
490

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: trunk/source/processes/electromagnetic/standard/test/testLowEnergyPAI.cc

Download in other formats: