Context Navigation

source: trunk/source/processes/electromagnetic/lowenergy/test/G4hTestStoppingPower.cc @ 1228

Last change on this file since 1228 was 1199, checked in by garnier, 15 years ago
nvx fichiers dans CVS
File size: 33.7 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: G4hTestStoppingPower.cc,v 1.18 2006/06/29 19:44:46 gunter Exp $
28	// GEANT4 tag $Name: geant4-09-03-cand-01 $
29	//
30	// -------------------------------------------------------------------
31	// GEANT 4 class file --- Copyright CERN 1998
32	// CERN Geneva Switzerland
33	//
34	//
35	// File name: G4hTestStoppingPower.cc
36	//
37	// Author: Vladimir Ivanchenko
38	//
39	// Creation date: 28 July 2000
40	//
41	// Modifications:
42	//
43	// -------------------------------------------------------------------
44
45	#include "globals.hh"
46	#include "G4ios.hh"
47	#include "G4UnitsTable.hh"
48	#include <fstream>
49	#include <iomanip>
50
51	#include "G4Material.hh"
52	#include "G4Element.hh"
53	#include "G4ProcessManager.hh"
54	#include "G4hLowEnergyIonisation.hh"
55	#include "G4hBetheBlochModel.hh"
56	#include "G4hParametrisedLossModel.hh"
57	#include "G4hNuclearStoppingModel.hh"
58	#include "G4QAOLowEnergyLoss.hh"
59	#include "G4VParticleChange.hh"
60	#include "G4ParticleChange.hh"
61	#include "G4DynamicParticle.hh"
62
63	#include "G4hZiegler1977p.hh"
64	#include "G4hZiegler1985p.hh"
65	#include "G4hZiegler1977He.hh"
66	#include "G4hICRU49p.hh"
67	#include "G4hICRU49He.hh"
68	#include "G4hSRIM2000p.hh"
69	#include "G4hSRIM2003p.hh"
70
71	#include "G4hZiegler1977Nuclear.hh"
72	#include "G4hZiegler1985Nuclear.hh"
73	//#include "G4hMollereNuclear.hh" // exist no more
74
75
76	#include "G4Box.hh"
77	#include "G4PVPlacement.hh"
78
79	#include "G4Step.hh"
80	#include "G4GRSVolume.hh"
81
82	#include "G4ParticleTypes.hh"
83	#include "G4ParticleTable.hh"
84	#include "G4ParticleDefinition.hh"
85	#include "G4ParticleWithCuts.hh"
86	#include "G4Ions.hh"
87
88	// New Histogramming (from AIDA and Anaphe):
89	#include <memory> // for the auto_ptr(T>
90
91	#include "AIDA/AIDA.h"
92
93
94	#include "hTest/include/G4IonC12.hh"
95	#include "hTest/include/G4IonAr40.hh"
96
97	#include "G4Timer.hh"
98
99	int main()
100	{
101	// ---- HBOOK initialization
102
103	G4String hFile = "htest.hbook";
104
105	//--------- Materials definition ---------
106
107	G4Material* Be = new G4Material("Beryllium", 4., 9.01g/mole, 1.848g/cm3);
108	G4Material* Graphite = new G4Material("Graphite",6.,12.0g/mole,2.265g/cm3);
109	G4Material* Al = new G4Material("Aluminium", 13., 26.98g/mole, 2.7 g/cm3);
110	G4Material* Si = new G4Material("Silicon", 14., 28.055g/mole, 2.33g/cm3);
111	G4Material* LAr = new G4Material("LArgon", 18., 39.95g/mole, 1.393g/cm3);
112	G4Material* Ti = new G4Material("Titan", 22., 47.867g/mole, 4.54g/cm3);
113	G4Material* Fe = new G4Material("Iron", 26., 55.85g/mole, 7.87g/cm3);
114	G4Material* Cu = new G4Material("Copper", 29., 63.55g/mole, 8.96g/cm3);
115	G4Material* W = new G4Material("Tungsten",74., 183.85g/mole, 19.30g/cm3);
116	G4Material* Pb = new G4Material("Lead", 82., 207.19g/mole, 11.35g/cm3);
117	G4Material* U = new G4Material("Uranium", 92., 238.03g/mole, 18.95g/cm3);
118
119	G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
120	G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
121	G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
122	G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
123	G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole);
124
125	// G4Material* water = new G4Material ("Water" ,"H_2O", 1.*g/cm3, 2);
126	G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2);
127	water->AddElement(H,2);
128	water->AddElement(O,1);
129
130	// G4Material* ethane = new G4Material ("Ethane" ,"C_2H_6", 0.4241*g/cm3, 2);
131	G4Material* ethane = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
132	ethane->AddElement(H,6);
133	ethane->AddElement(C,2);
134
135	// G4Material* csi = new G4Material ("CsI" , "CsI", 4.53*g/cm3, 2);
136	G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2);
137	csi->AddElement(Cs,1);
138	csi->AddElement(I,1);
139
140	G4Material* el[93];
141	G4double z;
142	G4int j;
143
144	for ( j=1; j<93; j++)
145	{
146	z = G4double(j);
147	el[j] = new G4Material ("atom", z , 2.0zg/mole, 10.0 *g/cm3);
148	}
149
150	const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
151
152	// create table
153	G4int numOfMaterials = G4Material::GetNumberOfMaterials();
154
155	G4double dimx = 1mm, dimy = 1mm, dimz = 1*mm;
156	G4int imat = 0;
157	G4cout << "Available materials are: " << G4endl;
158	for (imat = 0; imat < numOfMaterials; imat++) {
159	G4cout << imat << ") " << (*theMaterialTable)[imat]->GetName() << G4endl;
160	}
161
162
163	// Geometry definitions
164	G4Box* theFrame = new G4Box ("Frame",dimx, dimy, dimz);
165
166	G4LogicalVolume* LogicalFrame = new G4LogicalVolume(theFrame,
167	(*theMaterialTable)[0],
168	"LFrame", 0, 0, 0);
169
170	G4PVPlacement* PhysicalFrame = new G4PVPlacement(0,G4ThreeVector(),
171	"PFrame",LogicalFrame,0,false,0);
172
173	G4cout << PhysicalFrame << G4endl;
174
175	//--------- Particle definition ---------
176
177	G4ParticleDefinition* gamma = G4Gamma::Gamma();
178	G4ParticleDefinition* electron = G4Electron::Electron();
179	G4ParticleDefinition* proton = G4Proton::Proton();
180	G4ParticleDefinition* antiproton =G4AntiProton::AntiProton();
181	G4ParticleDefinition* deuteron = G4Deuteron::Deuteron();
182	G4ParticleDefinition* alpha = G4Alpha::Alpha();
183
184	G4Ions* iC12 = new G4Ions::G4Ions(
185	"IonC12", 11.14945GeV, 0.0MeV, +6.0*eplus,
186	0, +1, 0,
187	0, 0, 0,
188	"static_nucleus", 0, +12, 0,
189	true, -1.0, NULL);
190
191	G4Ions* iAr40 = new G4Ions::G4Ions(
192	"IonAr40", 37.291GeV, 0.0MeV, +18.0*eplus,
193	0, +1, 0,
194	0, 0, 0,
195	"static_nucleus", 0, +40, 0,
196	true, -1.0, NULL);
197
198	G4Ions* iFe56 = new G4Ions::G4Ions(
199	"IonFe56", 52.0308GeV, 0.0MeV, +26.0*eplus,
200	0, +1, 0,
201	0, 0, 0,
202	"static_nucleus", 0, +56, 0,
203	true, -1.0, NULL);
204
205	G4ParticleDefinition* ionC12 = iC12->IonsDefinition();
206	G4ParticleDefinition* ionAr40 = iAr40->IonsDefinition();
207	G4ParticleDefinition* ionFe56 = iFe56->IonsDefinition();
208
209	G4ParticleDefinition* part[7];
210	part[0] = proton;
211	part[1] = antiproton;
212	part[2] = deuteron;
213	part[3] = alpha;
214	part[4] = ionC12;
215	part[5] = ionAr40;
216	part[6] = ionFe56;
217
218	G4double ecut = 1000.0*mm;
219	G4double pcut = 0.0001*mm;
220	gamma->SetCuts(ecut);
221	electron->SetCuts(ecut);
222	G4cout << "Cuts are following: cutElectron = " << ecut
223	<< " mm; cutProton = " << pcut << " mm" << G4endl;
224
225	//--------- Ionisation processes definition and build physics table --
226
227	// Define models for parametrisation of electronic energy losses
228	// G4VLowEnergyModel* theBetheBlochModel =
229	// new G4hBetheBlochModel("Bethe-Bloch") ;
230	// theProtonModel = new G4hParametrisedLossModel(theProtonTable) ;
231	// theAntiProtonModel = new G4QAOLowEnergyLoss(theAntiProtonTable) ;
232
233	G4hNuclearStoppingModel* theNuclearStoppingModel =
234	// new G4hNuclearStoppingModel("Ziegler1977") ;
235	// new G4hNuclearStoppingModel("Ziegler1985") ;
236	new G4hNuclearStoppingModel("ICRU_R49") ;
237
238	G4VLowEnergyModel* theIonEffChargeModel =
239	new G4hIonEffChargeSquare("Ziegler1988") ;
240
241
242	G4hLowEnergyIonisation* hIon[8];
243	G4ProcessManager* theProcessManager[8];
244	G4int i;
245
246	G4cout << "Define processes!" << G4endl;
247
248	for( i=0; i<7; i++) {
249	G4cout << "Ionisation process for particle " << i << " " << part[i]->GetParticleName() << G4endl;
250	theProcessManager[i] = new G4ProcessManager(part[i]);
251	part[i]->SetProcessManager(theProcessManager[i]);
252	hIon[i] = new G4hLowEnergyIonisation();
253	hIon[i]->SetEnlossFluc(false) ;
254
255	// hIon[i]->SetBarkasOff();
256	// hIon[i]->SetNuclearStoppingOff();
257	// hIon[i]->SetStoppingPowerTableName("ICRU_R49p");
258
259	theProcessManager[i]->AddProcess(hIon[i]);
260	hIon[i]->BuildPhysicsTable(*part[i]);
261	}
262
263
264
265	//----------- Histogram -------------------------------------
266
267	G4cout << "Fill Hbook!" << G4endl;
268
269	// Creating the analysis factory
270	std::auto_ptr< AIDA::IAnalysisFactory > af( AIDA_createAnalysisFactory() );
271
272	// Creating the tree factory
273	std::auto_ptr< AIDA::ITreeFactory > tf( af->createTreeFactory() );
274
275	// Creating a tree mapped to a new hbook file.
276	std::auto_ptr< AIDA::ITree > tree( tf->create( hFile,"hbook", false,false) );
277	std::cout << "Tree store : " << tree->storeName() << std::endl;
278
279	// Creating a histogram factory, whose histograms will be handled by the tree
280	std::auto_ptr< AIDA::IHistogramFactory > hf( af->createHistogramFactory( *tree ) );
281
282	// G4Material* material ;
283
284	G4double minE = 1.0eV, maxE = 10000.0MeV, s;
285	const G4int num = 200;
286	G4double tkin = 0.0;
287	s = (std::log10(maxE)-std::log10(minE))/num;
288
289	AIDA::IHistogram1D* h[71] ;
290
291	// Test on Stopping Powers for all elements
292
293	h[1] = hf->createHistogram1D("1","p 40 keV (keV*cm2/10^15!atoms) Ziegler1977p"
294	,92,0.5,92.5) ;
295	h[2] = hf->createHistogram1D("2","p 100 keV (keV*cm2/10^15!atoms) Ziegler1977p"
296	,92,0.5,92.5) ;
297	h[3] = hf->createHistogram1D("3","p 400 keV (keV*cm2/10^15!atoms) Ziegler1977p"
298	,92,0.5,92.5) ;
299
300	h[4] = hf->createHistogram1D("4","p 1 MeV (keV*cm2/10^15!atoms) Ziegler1977p"
301	,92,0.5,92.5) ;
302	h[5] = hf->createHistogram1D("5","p 4 MeV (keV*cm2/10^15!atoms) Ziegler1977p"
303	,92,0.5,92.5) ;
304
305	h[6] = hf->createHistogram1D("6","p 40 keV (keV*cm2/10^15!atoms) ICRU_49p"
306	,92,0.5,92.5) ;
307	h[7] = hf->createHistogram1D("7","p 100 keV (keV*cm2/10^15!atoms) ICRU_49p"
308	,92,0.5,92.5) ;
309	h[8] = hf->createHistogram1D("8","p 400 keV (keV*cm2/10^15!atoms) ICRU_49p"
310	,92,0.5,92.5) ;
311	h[9] = hf->createHistogram1D("9","p 1 MeV (keV*cm2/10^15!atoms) ICRU_49p"
312	,92,0.5,92.5) ;
313	h[10] = hf->createHistogram1D("10","p 4 MeV (keV*cm2/10^15!atoms) ICRU_49p"
314	,92,0.5,92.5) ;
315
316	h[11] = hf->createHistogram1D("11","p 40 keV (keV*cm2/10^15!atoms) Ziegler1977He"
317	,92,0.5,92.5) ;
318	h[12] = hf->createHistogram1D("12","p 100 keV (keV*cm2/10^15!atoms) Ziegler1977He"
319	,92,0.5,92.5) ;
320	h[13] = hf->createHistogram1D("13","p 400 keV (keV*cm2/10^15!atoms) Ziegler1977He"
321	,92,0.5,92.5) ;
322	h[14] = hf->createHistogram1D("14","p 1 MeV (keV*cm2/10^15!atoms) Ziegler1977He"
323	,92,0.5,92.5) ;
324	h[15] = hf->createHistogram1D("15","p 4 MeV (keV*cm2/10^15!atoms) Ziegler1977He"
325	,92,0.5,92.5) ;
326
327	h[16] = hf->createHistogram1D("16","He 10 keV (keV*cm2/10^15!atoms) Ziegler1977He"
328	,92,0.5,92.5) ;
329	h[17] = hf->createHistogram1D("17","He 40 keV (keV*cm2/10^15!atoms) Ziegler1977He"
330	,92,0.5,92.5) ;
331	h[18] = hf->createHistogram1D("18","He 100 keV (keV*cm2/10^15!atoms) Ziegler1977He"
332	,92,0.5,92.5) ;
333
334	h[19] = hf->createHistogram1D("19","He 400 keV (keV*cm2/10^15!atoms) Ziegler1977He"
335	,92,0.5,92.5) ;
336	h[20] = hf->createHistogram1D("20","He 1 MeV (keV*cm2/10^15!atoms) Ziegler1977He"
337	,92,0.5,92.5) ;
338	h[21] = hf->createHistogram1D("21","He 4 MeV (keV*cm2/10^15!atoms) Ziegler1977He"
339	,92,0.5,92.5) ;
340	h[22] = hf->createHistogram1D("22","He 10 keV (keV*cm2/10^15!atoms) ICRU49He"
341	,92,0.5,92.5) ;
342	h[23] = hf->createHistogram1D("23","He 40 keV (keV*cm2/10^15!atoms) ICRU49He"
343	,92,0.5,92.5) ;
344	h[24] = hf->createHistogram1D("24","He 100 keV (keV*cm2/10^15!atoms) ICRU49He"
345	,92,0.5,92.5) ;
346	h[25] = hf->createHistogram1D("25","He 400 keV (keV*cm2/10^15!atoms) ICRU49He"
347	,92,0.5,92.5) ;
348	h[26] = hf->createHistogram1D("26","He 1 MeV (keV*cm2/10^15!atoms) ICRU49He"
349	,92,0.5,92.5) ;
350	h[27] = hf->createHistogram1D("27","He 4 MeV (keV*cm2/10^15!atoms) ICRU49He"
351	,92,0.5,92.5) ;
352
353	h[28]= hf->createHistogram1D("28","p in C (MeV/mm) ICRU49p"
354	,num,std::log10(minE),std::log10(maxE)) ;
355	h[29]= hf->createHistogram1D("29","p in Al (MeV/mm) ICRU49p"
356	,num,std::log10(minE),std::log10(maxE)) ;
357	h[30]= hf->createHistogram1D("30","p in Si (MeV/mm) ICRU49p"
358	,num,std::log10(minE),std::log10(maxE)) ;
359	h[31]= hf->createHistogram1D("31","p in Cu (MeV/mm) ICRU49p"
360	,num,std::log10(minE),std::log10(maxE)) ;
361	h[32]= hf->createHistogram1D("32","p in Fe (MeV/mm) ICRU49p"
362	,num,std::log10(minE),std::log10(maxE)) ;
363	h[33]= hf->createHistogram1D("33","p in Pb (MeV/mm) ICRU49p"
364	,num,std::log10(minE),std::log10(maxE)) ;
365	h[34]= hf->createHistogram1D("34","p in C2H6 (MeV/mm) ICRU49p"
366	,num,std::log10(minE),std::log10(maxE)) ;
367	h[35]= hf->createHistogram1D("35","p in H2O (MeV/mm) ICRU49p"
368	,num,std::log10(minE),std::log10(maxE)) ;
369	h[36]= hf->createHistogram1D("36","p in lAr (MeV/mm) ICRU49p"
370	,num,std::log10(minE),std::log10(maxE)) ;
371	h[37]= hf->createHistogram1D("37","p in CsI (MeV/mm) ICRU49p"
372	,num,std::log10(minE),std::log10(maxE)) ;
373
374
375	h[38]= hf->createHistogram1D("38","p in C (MeV/mm)Ziegler1985p"
376	,num,std::log10(minE),std::log10(maxE)) ;
377	h[39]= hf->createHistogram1D("39","p in Al (MeV/mm)Ziegler1985p"
378	,num,std::log10(minE),std::log10(maxE)) ;
379	h[40]= hf->createHistogram1D("40","p in Si (MeV/mm)Ziegler1985p"
380	,num,std::log10(minE),std::log10(maxE)) ;
381	h[41]= hf->createHistogram1D("41","p in Cu (MeV/mm)Ziegler1985p"
382	,num,std::log10(minE),std::log10(maxE)) ;
383	h[42]= hf->createHistogram1D("42","p in Fe (MeV/mm)Ziegler1985p"
384	,num,std::log10(minE),std::log10(maxE)) ;
385	h[43]= hf->createHistogram1D("43","p in Pb (MeV/mm)Ziegler1985p"
386	,num,std::log10(minE),std::log10(maxE)) ;
387	h[44]= hf->createHistogram1D("44","p in C2H6 (MeV/mm)Ziegler1985p"
388	,num,std::log10(minE),std::log10(maxE)) ;
389	h[45]= hf->createHistogram1D("45","p in H2O (MeV/mm)Ziegler1985p"
390	,num,std::log10(minE),std::log10(maxE)) ;
391	h[46]= hf->createHistogram1D("46","p in lAr (MeV/mm)Ziegler1985p"
392	,num,std::log10(minE),std::log10(maxE)) ;
393	h[47]= hf->createHistogram1D("47","p in CsI (MeV/mm)Ziegler1985p"
394	,num,std::log10(minE),std::log10(maxE)) ;
395
396
397	h[48]= hf->createHistogram1D("48","He effective charge for Cu"
398	,num,std::log10(minE),std::log10(maxE)) ;
399	h[49]= hf->createHistogram1D("49","C12 effective charge in Cu"
400	,num,std::log10(minE),std::log10(maxE)) ;
401
402	h[50]= hf->createHistogram1D("50","He in Al (MeV/(mg/cm2)) ICRU49p"
403	,num,std::log10(minE),std::log10(maxE)) ;
404	h[51]= hf->createHistogram1D("51","C12 in Al (MeV/(mg/cm2)) ICRU49p"
405	,num,std::log10(minE),std::log10(maxE)) ;
406	h[52]= hf->createHistogram1D("52","Ar40 in Al (MeV/(mg/cm2)) ICRU49p"
407	,num,std::log10(minE),std::log10(maxE)) ;
408
409	// Table with the data
410	h[53] = hf->createHistogram1D("53","Data p 40 keV (keV*cm2/10^15!atoms) Ziegler1977p"
411	,92,0.5,92.5) ;
412	h[54] = hf->createHistogram1D("54","Data He 40 keV (keV*cm2/10^15!atoms) Ziegler1977He"
413	,92,0.5,92.5) ;
414
415	h[55] = hf->createHistogram1D("55","p 40 keV (keV*cm2/10^15!atoms) Ziegler1985p"
416	,92,0.5,92.5) ;
417	h[56] = hf->createHistogram1D("56","p 100 keV (keV*cm2/10^15!atoms) Ziegler1985p"
418	,92,0.5,92.5) ;
419	h[57] = hf->createHistogram1D("57","p 400 keV (keV*cm2/10^15!atoms) Ziegler1985p"
420	,92,0.5,92.5) ;
421	h[58] = hf->createHistogram1D("58","p 1 MeV (keV*cm2/10^15!atoms) Ziegler1985p"
422	,92,0.5,92.5) ;
423	h[59] = hf->createHistogram1D("59","p 4 MeV (keV*cm2/10^15!atoms) Ziegler1985p"
424	,92,0.5,92.5) ;
425	// Histo for Antiproton
426
427	h[60]= hf->createHistogram1D("60","pbar in C (MeV/mm) QAOLoss"
428	,num,std::log10(minE),std::log10(maxE)) ;
429	h[61]= hf->createHistogram1D("61","pbar in Al (MeV/mm) QAOLoss"
430	,num,std::log10(minE),std::log10(maxE)) ;
431	h[62]= hf->createHistogram1D("62","pbar in Si (MeV/mm) QAOLoss"
432	,num,std::log10(minE),std::log10(maxE)) ;
433	h[63]= hf->createHistogram1D("63","pbar in Cu (MeV/mm) QAOLoss"
434	,num,std::log10(minE),std::log10(maxE)) ;
435	h[64]= hf->createHistogram1D("64","pbar in Fe (MeV/mm) QAOLoss"
436	,num,std::log10(minE),std::log10(maxE)) ;
437	h[65]= hf->createHistogram1D("65","pbar in Pb (MeV/mm) QAOLoss"
438	,num,std::log10(minE),std::log10(maxE)) ;
439	h[66]= hf->createHistogram1D("66","pbar in C2H6 (MeV/mm) QAOLoss"
440	,num,std::log10(minE),std::log10(maxE)) ;
441	h[67]= hf->createHistogram1D("67","pbar in H2O (MeV/mm) QAOLoss"
442	,num,std::log10(minE),std::log10(maxE)) ;
443	h[68]= hf->createHistogram1D("68","pbar in lAr (MeV/mm) QAOLoss"
444	,num,std::log10(minE),std::log10(maxE)) ;
445	h[69]= hf->createHistogram1D("69","pbar in CsI (MeV/mm) QAOLoss"
446	,num,std::log10(minE),std::log10(maxE)) ;
447
448	h[70] = hf->createHistogram1D("70","p 6.5 MeV (keV*cm2/10^15!atoms) Ziegler77p"
449	,92,0.5,92.5) ;
450
451
452	// G4VhElectronicStoppingPower* Z77p = new G4hZiegler1977p() ;
453	// G4VhElectronicStoppingPower* Z85p = new G4hZiegler1985p() ;
454	G4VhElectronicStoppingPower* Z77p = new G4hSRIM2000p() ;
455	G4VhElectronicStoppingPower* Z85p = new G4hSRIM2003p() ;
456	G4VhElectronicStoppingPower* Z77He = new G4hZiegler1977He() ;
457	G4VhElectronicStoppingPower* I49p = new G4hICRU49p() ;
458	G4VhElectronicStoppingPower* I49He = new G4hICRU49He() ;
459
460	G4double de;
461
462	// Ziegler1977p
463	G4double tau = 40.0*keV ;
464	G4double q2;
465	G4double rateMass=4.0026/1.007276;
466
467	for ( j=1; j<93; j++)
468	{
469	z = G4double(j);
470	q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin);
471	de = Z77p->ElectronicStoppingPower(z, tau) ;
472	h[1]->fill(z,de) ;
473	de = Z85p->ElectronicStoppingPower(z, tau) ;
474	h[55]->fill(z,de) ;
475	de = I49p->ElectronicStoppingPower(z, tau) ;
476	h[6]->fill(z,de) ;
477	de = Z77He->ElectronicStoppingPower(z, tau) ;
478	h[11]->fill(z,de) ;
479	de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ;
480	h[17]->fill(z,de) ;
481	de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ;
482	h[23]->fill(z,de) ;
483	}
484	tau = 100.0*keV ;
485	for ( j=1; j<93; j++)
486	{
487	z = G4double(j);
488	de = Z77p->ElectronicStoppingPower(z, tau) ;
489	h[2]->fill(z,de) ;
490	de = Z85p->ElectronicStoppingPower(z, tau) ;
491	h[56]->fill(z,de) ;
492	de = I49p->ElectronicStoppingPower(z, tau) ;
493	h[7]->fill(z,de) ;
494	de = Z77He->ElectronicStoppingPower(z, tau) ;
495	h[12]->fill(z,de) ;
496	q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin);
497	de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ;
498	h[18]->fill(z,de) ;
499	de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ;
500	h[24]->fill(z,de) ;
501	}
502
503	tau = 400.0*keV ;
504	for ( j=1; j<93; j++)
505	{
506	z = G4double(j);
507	de = Z77p->ElectronicStoppingPower(z, tau) ;
508	h[3]->fill(z,de) ;
509	de = Z85p->ElectronicStoppingPower(z, tau) ;
510	h[57]->fill(z,de) ;
511	de = I49p->ElectronicStoppingPower(z, tau) ;
512	h[8]->fill(z,de) ;
513	de = Z77He->ElectronicStoppingPower(z, tau) ;
514	h[13]->fill(z,de) ;
515	q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin);
516	de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ;
517	h[19]->fill(z,de) ;
518	de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ;
519	h[25]->fill(z,de) ;
520	}
521	tau = 1000.0*keV ;
522	for ( j=1; j<93; j++)
523	{
524	z = G4double(j);
525	de = Z77p->ElectronicStoppingPower(z, tau) ;
526	h[4]->fill(z,de) ;
527	de = Z85p->ElectronicStoppingPower(z, tau) ;
528	h[58]->fill(z,de) ;
529	de = I49p->ElectronicStoppingPower(z, tau) ;
530	h[9]->fill(z,de) ;
531	de = Z77He->ElectronicStoppingPower(z, tau) ;
532	h[14]->fill(z,de) ;
533	q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin);
534	de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ;
535	h[20]->fill(z,de) ;
536	de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ;
537	h[26]->fill(z,de) ;
538	}
539	tau = 4000*keV ;
540	for ( j=1; j<93; j++)
541	{
542	z = G4double(j);
543	de = Z77p->ElectronicStoppingPower(z, tau) ;
544	h[5]->fill(z,de) ;
545	de = Z85p->ElectronicStoppingPower(z, tau) ;
546	h[59]->fill(z,de) ;
547	de = I49p->ElectronicStoppingPower(z, tau) ;
548	h[10]->fill(z,de) ;
549	de = Z77He->ElectronicStoppingPower(z, tau) ;
550	h[15]->fill(z,de) ;
551	q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin);
552	de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ;
553	h[21]->fill(z,de) ;
554	de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ;
555	h[27]->fill(z,de) ;
556	}
557	tau = 6.5*MeV ;
558	for ( j=1; j<93; j++)
559	{
560	z = G4double(j);
561	de = Z77p->ElectronicStoppingPower(z, tau) ;
562	h[70]->fill(z,de) ;
563	}
564	tau = 10.0*keV ;
565	for ( j=1; j<93; j++)
566	{
567	z = G4double(j);
568	q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin);
569	de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ;
570	h[16]->fill(z,de) ;
571	de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ;
572	h[22]->fill(z,de) ;
573	}
574
575	for ( j=1; j<93; j++)
576	{
577	static G4double p40[92] = {
578	6.22, 6.55, 7.61, 10.4, 12.9, 13.9, 15.9, 14.6, 11.7, 11.1,
579	14.2, 20.6, 20.7, 22.3, 18.1, 19.3, 27.5, 30.5, 27.9, 29.8,
580	28.3, 26.7, 25.1, 22.5, 19.8, 20.1, 18.0, 20.1, 20.4, 23.4,
581	27.5, 29.0, 29.3, 31.0, 31.1, 34.8, 31.5, 35.0, 35.5, 37.3,
582	38.3, 35.9, 38.0, 34.4, 33.4, 29.8, 30.9, 32.7, 34.9, 36.0,
583	40.9, 39.1, 43.1, 45.8, 40.8, 44.1, 44.9, 42.0, 41.0, 40.0,
584	39.0, 38.0, 37.1, 38.2, 35.3, 31.5, 29.9, 29.1, 28.3, 27.5,
585	28.1, 28.9, 27.3, 26.3, 29.9, 29.1, 28.4, 25.8, 28.0, 24.9,
586	27.3, 30.7, 34.3, 35.4, 35.6, 35.5, 39.8, 42.9, 43.7, 44.1,
587	42.4, 41.8} ;
588	h[53]->fill(double(j),p40[j-1]) ;
589
590	static G4double he40[92] = {
591	11.8, 16.7, 21.9, 24.1, 34.7, 36.1, 45.5, 46.1, 45.8, 44.9,
592	46.9, 51.7, 54.9, 63.2, 63.8, 67.4, 79.8, 80.8, 78.0, 83.2,
593	85.4, 84.2, 90.6, 85.7, 84.1, 86.4, 82.0, 77.7, 74.4, 75.1,
594	75.6, 80.7, 84.9, 87.0, 88.6, 103.0, 103.0, 113.0, 116.0, 123.0,
595	120.0, 115.0, 118.0, 115.0, 113.0, 112.0, 111.0, 110.0, 116.0, 117.0,
596	119.0, 123.0, 122.0, 139.0, 138.0, 143.0, 152.0, 141.0, 139.0, 137.0,
597	135.0, 134.0, 130.0, 137.0, 129.0, 128.0, 124.0, 125.0, 120.0, 118.0,
598	119.0, 120.0, 121.0, 122.0, 125.0, 127.0, 128.0, 120.0, 125.0, 128.0,
599	136.0, 138.0, 144.0, 148.0, 151.0, 162.0, 166.0, 177.0, 179.0, 183.0,
600	177.0, 176.0 } ;
601	h[54]->fill(double(j),he40[j-1]) ;
602	}
603
604	G4cout << "Stopping Power histograms are filled!" << G4endl;
605
606	// dedx
607	for (j = 0 ; j < num-1 ; j++) {
608	tkin = std::pow(10.0,(std::log10(minE) + (G4double(j)+0.5)*s));
609	de = hIon[0]->ComputeDEDX(part[0],Graphite,tkin) ;
610	h[28]->fill(std::log10(tkin),de) ;
611	de = hIon[0]->ComputeDEDX(part[0],Al,tkin) ;
612	h[29]->fill(std::log10(tkin),de) ;
613	de = hIon[0]->ComputeDEDX(part[0],Si,tkin) ;
614	h[30]->fill(std::log10(tkin),de) ;
615	de = hIon[0]->ComputeDEDX(part[0],Cu,tkin) ;
616	h[31]->fill(std::log10(tkin),de) ;
617	de = hIon[0]->ComputeDEDX(part[0],Fe,tkin) ;
618	h[32]->fill(std::log10(tkin),de) ;
619	de = hIon[0]->ComputeDEDX(part[0],Pb,tkin) ;
620	h[33]->fill(std::log10(tkin),de) ;
621	de = hIon[0]->ComputeDEDX(part[0],ethane,tkin) ;
622	// G4cout << "ethane: E = " << tkin << "; dedx = " << de << G4endl;
623	h[34]->fill(std::log10(tkin),de) ;
624	de = hIon[0]->ComputeDEDX(part[0],water,tkin) ;
625	// de += theNuclearStoppingModel->TheValue(part[1],water,tkin) ;
626	// G4cout << "water : E = " << tkin << "; dedx = " << de << G4endl;
627	h[35]->fill(std::log10(tkin),de) ;
628	de = hIon[0]->ComputeDEDX(part[0],LAr,tkin) ;
629	h[36]->fill(std::log10(tkin),de) ;
630	de = hIon[0]->ComputeDEDX(part[0],csi,tkin) ;
631	h[37]->fill(std::log10(tkin),de) ;
632	}
633
634	G4cout << "Proton's dEdx histograms are filled!" << G4endl;
635
636	for (j = 0 ; j < num-1 ; j++) {
637	tkin = std::pow(10.0,(std::log10(minE) + (G4double(j)+0.5)*s));
638	de = hIon[0]->ComputeDEDX(part[1],Graphite,tkin) ;
639	h[60]->fill(std::log10(tkin),de) ;
640	de = hIon[0]->ComputeDEDX(part[1],Al,tkin) ;
641	h[61]->fill(std::log10(tkin),de) ;
642	de = hIon[0]->ComputeDEDX(part[1],Si,tkin) ;
643	h[62]->fill(std::log10(tkin),de) ;
644	de = hIon[0]->ComputeDEDX(part[1],Cu,tkin) ;
645	h[63]->fill(std::log10(tkin),de) ;
646	de = hIon[0]->ComputeDEDX(part[1],Fe,tkin) ;
647	h[64]->fill(std::log10(tkin),de) ;
648	de = hIon[0]->ComputeDEDX(part[1],Pb,tkin) ;
649	h[65]->fill(std::log10(tkin),de) ;
650	de = hIon[0]->ComputeDEDX(part[1],ethane,tkin) ;
651	// G4cout << "ethane: E = " << tkin << "; dedx = " << de << G4endl;
652	h[66]->fill(std::log10(tkin),de) ;
653	de = hIon[0]->ComputeDEDX(part[1],water,tkin) ;
654	// G4cout << "water : E = " << tkin << "; dedx = " << de << G4endl;
655	h[67]->fill(std::log10(tkin),de) ;
656	de = hIon[0]->ComputeDEDX(part[1],LAr,tkin) ;
657	h[68]->fill(std::log10(tkin),de) ;
658	de = hIon[0]->ComputeDEDX(part[1],csi,tkin) ;
659	h[69]->fill(std::log10(tkin),de) ;
660	}
661
662	G4cout << "AntiProton's dEdx histograms are filled!" << G4endl;
663
664	G4double mProt = part[0]->GetPDGMass()*1.007276;
665	G4double fact = cm/(2700.0*MeV) ; // to MeV/mg/cm^2
666
667	for (j = 0 ; j < num-1 ; j++) {
668	tkin = std::pow(10.0,(std::log10(minE) + (G4double(j)+0.5)*s));
669	de = theIonEffChargeModel->TheValue(part[3],Cu,tkin) ;
670	// G4cout << "E = " << tkin << "; dedx = " << de << G4endl;
671	h[48]->fill(std::log10(tkin),de) ;
672	de = theIonEffChargeModel->TheValue(part[4],Cu,tkin) ;
673	h[49]->fill(std::log10(tkin),de) ;
674	G4double tRed = tkin * (part[3]->GetPDGMass())/mProt ;
675	de = hIon[3]->ComputeDEDX(part[3],Al,tRed) ;
676	de += theNuclearStoppingModel->TheValue(part[3],Al,tRed) ;
677	h[50]->fill(std::log10(tkin),de*fact) ;
678	tRed = tkin * (part[4]->GetPDGMass())/mProt ;
679	de = hIon[4]->ComputeDEDX(part[4],Al,tRed) ;
680	//de += theNuclearStoppingModel->TheValue(part[4],Al,tRed) ;
681	h[51]->fill(std::log10(tkin),de*fact) ;
682	tRed = tkin * (part[5]->GetPDGMass())/mProt ;
683	de = hIon[5]->ComputeDEDX(part[5],Al,tRed) ;
684	//de += theNuclearStoppingModel->TheValue(part[5],Al,tRed) ;
685	h[52]->fill(std::log10(tkin),de*fact) ;
686	}
687
688	G4cout << "Ions's dEdx histograms are filled!" << G4endl;
689
690	theProcessManager[7] = new G4ProcessManager(part[0]);
691	part[0]->SetProcessManager(theProcessManager[7]);
692	hIon[7] = new G4hLowEnergyIonisation();
693	hIon[7]->SetElectronicStoppingPowerModel(part[0],"Ziegler1985p");
694	hIon[7]->SetEnlossFluc(false) ;
695	//hIon[7]->SetNuclearStoppingOff();
696	// hIon[7]->SetBarkasOff();
697	theProcessManager[7]->AddProcess(hIon[7]);
698	hIon[7]->BuildPhysicsTable(*part[0]);
699
700	G4cout << "Ziegler's dEdx histograms will be filled" << G4endl;
701
702	for (j = 0 ; j < num-1 ; j++) {
703	tkin = std::pow(10.0,(std::log10(minE) + (G4double(j)+0.5)*s));
704	de = hIon[7]->ComputeDEDX(part[0],Graphite,tkin) ;
705	h[38]->fill(std::log10(tkin),de) ;
706	de = hIon[7]->ComputeDEDX(part[0],Al,tkin) ;
707	h[39]->fill(std::log10(tkin),de) ;
708	de = hIon[7]->ComputeDEDX(part[0],Si,tkin) ;
709	h[40]->fill(std::log10(tkin),de) ;
710	de = hIon[7]->ComputeDEDX(part[0],Cu,tkin) ;
711	h[41]->fill(std::log10(tkin),de) ;
712	de = hIon[7]->ComputeDEDX(part[0],Fe,tkin) ;
713	h[42]->fill(std::log10(tkin),de) ;
714	de = hIon[7]->ComputeDEDX(part[0],Pb,tkin) ;
715	h[43]->fill(std::log10(tkin),de) ;
716	de = hIon[7]->ComputeDEDX(part[0],ethane,tkin) ;
717	// G4cout << "ethane: E = " << tkin << "; dedx = " << de << G4endl;
718	h[44]->fill(std::log10(tkin),de) ;
719	de = hIon[7]->ComputeDEDX(part[0],water,tkin) ;
720	// G4cout << "water: E = " << tkin << "; dedx = " << de << G4endl;
721	h[45]->fill(std::log10(tkin),de) ;
722	de = hIon[7]->ComputeDEDX(part[0],LAr,tkin) ;
723	h[46]->fill(std::log10(tkin),de) ;
724	de = hIon[7]->ComputeDEDX(part[0],csi,tkin) ;
725	h[47]->fill(std::log10(tkin),de) ;
726	}
727
728	//---------------------- Fill Ntuple ------------------------
729	/*
730	G4cout << "Fill Ntuple!" << G4endl;
731
732	// ---- primary ntuple ------
733	HepTuple* ntuple = hbookManager->ntuple("dEdx ntuple");
734	assert (ntuple != 0);
735
736	for( i = 0; i < num; i++) {
737	tkin = std::pow(10,(std::log10(minE) + i*s));
738
739	for ( G4int j = 0 ; j < numOfMaterials; j++ ) {
740
741	material = (*theMaterialTable)[j] ;
742	// get elements in the actual material,
743	const G4ElementVector* theElementVector = material->GetElementVector() ;
744	const G4double* theAtomicNumDensityVector =
745	material->GetAtomicNumDensityVector() ;
746	const G4int NumberOfElements = material->GetNumberOfElements() ;
747
748	// loop for the elements in the material
749	// to find out average values Z, vF, lF
750	G4double z = 0.0, norm = 0.0 ;
751
752	if( 1 == NumberOfElements ) {
753	z = material->GetZ() ;
754
755	} else {
756	for (G4int iel=0; iel<NumberOfElements; iel++) {
757	const G4Element* element = (*theElementVector)(iel) ;
758	G4double z2 = element->GetZ() ;
759	const G4double weight = theAtomicNumDensityVector[iel] ;
760	norm += weight ;
761	z += z2 * weight ;
762	}
763	z /= norm ;
764	}
765
766	for (G4int k=0 ; k<7; k++) {
767	G4double dedx = hIon[k]->ComputeDEDX(part[k],material,tkin) ;
768	G4double q = theIonEffChargeModel->TheValue(part[k],material,tkin);
769	ntuple->column("part",k);
770	ntuple->column("mat",j);
771	ntuple->column("ie",i);
772	ntuple->column("zeff",z);
773	ntuple->column("tkin",tkin/MeV);
774	ntuple->column("mass",(part[k]->GetPDGMass())/MeV);
775	ntuple->column("char",(part[k]->GetPDGCharge())/eplus);
776	ntuple->column("cha2",q);
777	ntuple->column("dedx",dedx*mm/MeV);
778	ntuple->dumpData();
779	}
780	}
781	}
782	*/
783	//----------- End of work -------------------------------------
784
785	std::cout << "Committing..." << std::endl;
786	tree->commit();
787	std::cout << "Closing the tree..." << std::endl;
788	tree->close();
789
790	G4cout << "Ntuple and Hbook are saved" << G4endl;
791
792	// delete materials and elements
793	delete Be;
794	delete Graphite;
795	delete Al;
796	delete Si;
797	delete LAr;
798	delete Fe;
799	delete Cu;
800	delete W;
801	delete Pb;
802	delete U;
803	delete H;
804	delete C;
805	delete Cs;
806	delete I;
807	delete Ti;
808	delete O;
809	delete water;
810	delete ethane;
811	delete csi;
812	G4cout << "Materials are deleted" << G4endl;
813	delete theIonEffChargeModel;
814	delete Z77p;
815	delete Z77He;
816	delete I49p;
817	delete I49He;
818
819	cout<<"END OF THE MAIN PROGRAM"<<G4endl;
820	}

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

Download in other formats: