Context Navigation

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

Last change on this file since 1228 was 1199, checked in by garnier, 15 years ago
nvx fichiers dans CVS
File size: 27.2 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	// GEANT 4 class file --- Copyright CERN 1998
30	// CERN Geneva Switzerland
31	//
32	//
33	// File name: G4StoppingPowerTest
34	// This test provide cross sections
35	// tests for electromagnetic processes. The input
36	// data have to be describe in ASCII file
37	//
38	// Author: V.Ivanchenko
39	//
40	// Creation date: 23 May 2001
41	//
42	// Modifications:
43	//
44	// -------------------------------------------------------------------
45
46	#include "globals.hh"
47	#include "G4ios.hh"
48	#include <fstream>
49	#include <iomanip>
50
51	#include "G4Material.hh"
52	#include "G4VContinuousDiscreteProcess.hh"
53	#include "G4ProcessManager.hh"
54	#include "G4VParticleChange.hh"
55
56	#include "G4LowEnergyIonisation.hh"
57	#include "G4LowEnergyBremsstrahlung.hh"
58	#include "G4LowEnergyCompton.hh"
59	#include "G4LowEnergyGammaConversion.hh"
60	#include "G4LowEnergyPhotoElectric.hh"
61	#include "G4LowEnergyRayleigh.hh"
62	#include "G4hLowEnergyIonisation.hh"
63
64	#include "G4VeEnergyLoss.hh"
65	#include "G4eIonisation.hh"
66	#include "G4eBremsstrahlung.hh"
67	#include "G4ComptonScattering.hh"
68	#include "G4GammaConversion.hh"
69	#include "G4PhotoElectricEffect.hh"
70
71	#include "G4eplusAnnihilation.hh"
72
73	#include "G4MuIonisation.hh"
74	#include "G4MuBremsstrahlung.hh"
75	#include "G4MuPairProduction.hh"
76
77	#include "G4hIonisation.hh"
78	#include "G4MultipleScattering.hh"
79	#include "G4EnergyLossTables.hh"
80	#include "G4ParticleChange.hh"
81	#include "G4ParticleChange.hh"
82	#include "G4DynamicParticle.hh"
83	#include "G4AntiProton.hh"
84	#include "G4Proton.hh"
85	#include "G4Electron.hh"
86	#include "G4Positron.hh"
87	#include "G4Gamma.hh"
88	#include "G4ForceCondition.hh"
89	#include "G4Box.hh"
90	#include "G4PVPlacement.hh"
91	#include "G4Step.hh"
92	#include "G4GRSVolume.hh"
93
94	// New Histogramming (from AIDA and Anaphe):
95	#include <memory> // for the auto_ptr(T>
96
97
98	#include "AIDA/AIDA.h"
99
100	#include "hTest/include/G4IonC12.hh"
101	#include "hTest/include/G4IonAr40.hh"
102
103	#include "G4Timer.hh"
104
105	int main(int argc,char** argv)
106	{
107	// HepTupleManager* hbookManager;
108
109	// -------------------------------------------------------------------
110	// Setup
111
112	G4int nPart =-1;
113	G4String nameMat = "Si";
114	G4int nProcess = 3;
115	G4bool usepaw = false;
116	G4bool fluct = false;
117	G4bool lowE = true;
118	G4int verbose = 0;
119	G4double emin = 0.01*MeV;
120	G4double emax = 100.0*MeV;
121	G4int nstatf = 10;
122	G4double xstatf = 1.0/(G4double)nstatf;
123	G4int nbin = 1000;
124	G4String hFile = "hbook.paw";
125	G4double theStep = 0.01*micrometer;
126	G4double range = 1.0*micrometer;
127	G4double cutG = 10.0*mm;
128	G4double cutE = 10.0*mm;
129	G4Material* material = 0;
130	G4String name[3] = {"Ionisation", "Bremsstrahlung",
131	"Ionisation+Bremsstrahlung"};
132	G4bool setBarkasOff = false;
133	G4bool setNuclearOff= true;
134
135	G4cout.setf( ios::scientific, ios::floatfield );
136
137
138	// -------------------------------------------------------------------
139	// Control on input
140
141	if(argc < 2) {
142	G4cout << "Input file is not specified! Exit" << G4endl;
143	exit(1);
144	}
145
146	ifstream* fin = new ifstream();
147	string fname = argv[1];
148	fin->open(fname.c_str());
149	if( !fin->is_open()) {
150	G4cout << "Input file <" << fname << "> does not exist! Exit" << G4endl;
151	exit(1);
152	}
153
154	// -------------------------------------------------------------------
155	//--------- Materials definition ---------
156
157	G4Material* ma[16];
158	ma[0] = new G4Material("Be", 4., 9.01g/mole, 1.848g/cm3);
159	ma[1] = new G4Material("Graphite",6., 12.00g/mole, 2.265g/cm3 );
160	ma[1]->SetChemicalFormula("Graphite");
161	ma[2] = new G4Material("Al", 13., 26.98g/mole, 2.7 g/cm3);
162	ma[3] = new G4Material("Si", 14., 28.055g/mole, 2.33g/cm3);
163
164	ma[4] = new G4Material("LAr", 18., 39.95g/mole, 1.393g/cm3);
165	ma[5] = new G4Material("Fe", 26., 55.85g/mole, 7.87g/cm3);
166	ma[6] = new G4Material("Cu", 29., 63.55g/mole, 8.96g/cm3);
167	ma[7] = new G4Material("W", 74., 183.85g/mole, 19.30g/cm3);
168	ma[8] = new G4Material("Pb",82., 207.19g/mole, 11.35g/cm3);
169	ma[9] = new G4Material("U", 92., 238.03g/mole, 18.95g/cm3);
170
171	G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
172	G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
173	G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
174	G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
175	G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole);
176
177	ma[10] = new G4Material("O2", 8., 16.00g/mole, 1.1g/cm3);
178	ma[10]->SetChemicalFormula("O_2");
179
180	ma[11] = new G4Material ("Water" , 1.*g/cm3, 2);
181	ma[11]->AddElement(H,2);
182	ma[11]->AddElement(O,1);
183	ma[11]->SetChemicalFormula("H_2O");
184
185	ma[12] = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
186	ma[12]->AddElement(H,6);
187	ma[12]->AddElement(C,2);
188	ma[12]->SetChemicalFormula("C_2H_6");
189
190	ma[13] = new G4Material ("CsI" , 4.53*g/cm3, 2);
191	ma[13]->AddElement(Cs,1);
192	ma[13]->AddElement(I,1);
193	ma[13]->SetChemicalFormula("CsI");
194
195	ma[14] = new G4Material("H2", 1., 1.00794g/mole, 1.g/cm3);
196	ma[14]->SetChemicalFormula("H_2");
197	ma[15] = new G4Material("Au", 79., 196.96655g/mole, 19.32g/cm3);
198
199
200	static const G4MaterialTable* theMaterialTable =
201	G4Material::GetMaterialTable();
202
203	G4int nMaterials = G4Material::GetNumberOfMaterials();
204	G4cout << "Available materials are: " << G4endl;
205	G4int mat;
206	for (mat = 0; mat < nMaterials; mat++) {
207	G4cout << mat << ") " << ma[mat]->GetName() << G4endl;
208	}
209
210	G4cout << "Available processes are: " << G4endl;
211	for (mat = 0; mat < 2; mat++) {
212	G4cout << mat << ") " << name[mat] << G4endl;
213	}
214
215	// Particle definitions
216
217	G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
218	G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
219	G4ParticleDefinition* positron = G4Positron::PositronDefinition();
220	G4ParticleDefinition* proton = G4Proton::ProtonDefinition();
221	G4ParticleDefinition* antiproton = G4AntiProton::AntiProtonDefinition();
222	G4ParticleDefinition* c12 = G4IonC12::IonC12Definition();
223	G4ParticleDefinition* ar40 = G4IonAr40::IonAr40Definition();
224	G4ParticleDefinition* part = electron;
225
226	G4hLowEnergyIonisation* hionle = 0;
227	G4hLowEnergyIonisation* ionle = 0;
228	G4LowEnergyIonisation* eionle = 0;
229	G4LowEnergyBremsstrahlung* ebrle = 0;
230	G4hIonisation* ionst = 0;
231	G4hIonisation* hionst = 0;
232	G4eIonisation* eionst = 0;
233	G4eBremsstrahlung* ebrst = 0;
234
235	G4cout << "Process is initialized" << G4endl;;
236
237	// Geometry
238
239	G4double initX = 0.;
240	G4double initY = 0.;
241	G4double initZ = 1.;
242	G4double dimX = 100.0*cm;
243	G4double dimY = 100.0*cm;
244	G4double dimZ = 100.0*cm;
245
246	G4Box* sFrame = new G4Box ("Box",dimX, dimY, dimZ);
247	G4LogicalVolume* lFrame = new G4LogicalVolume(sFrame,material,"Box",0,0,0);
248	G4PVPlacement* pFrame = new G4PVPlacement(0,G4ThreeVector(),"Box",
249	lFrame,0,false,0);
250
251	// -------------------------------------------------------------------
252	// ---- Read input file
253	G4cout << "Available commands are: " << G4endl;
254	G4cout << "#events" << G4endl;
255	G4cout << "#particle" << G4endl;
256	G4cout << "#emin(MeV)" << G4endl;
257	G4cout << "#emax(MeV)" << G4endl;
258	G4cout << "#nbin" << G4endl;
259	G4cout << "#cutG(mm)" << G4endl;
260	G4cout << "#cutE(mm)" << G4endl;
261	G4cout << "#range(mm)" << G4endl;
262	G4cout << "#step(mm)" << G4endl;
263	G4cout << "#material" << G4endl;
264	G4cout << "#process" << G4endl;
265	G4cout << "#domain" << G4endl;
266	G4cout << "#paw" << G4endl;
267	G4cout << "#verbose" << G4endl;
268	G4cout << "#run" << G4endl;
269	G4cout << "#exit" << G4endl;
270	G4cout << "#barkas" << G4endl;
271	G4cout << "#nuclear" << G4endl;
272	G4cout << pFrame << G4endl;
273
274	G4ProcessManager gmanager, elecManager, *positManager,
275	protManager, aprotManager, *ionManager;
276
277	string line, line1;
278	G4bool end = true;
279	for(G4int run=0; run<100; run++) {
280	do {
281	(*fin) >> line;
282	G4cout << "Next line " << line << G4endl;
283	if(line == "#particle") {
284	(*fin) >> nPart;
285	} else if(line == "#emin(MeV)") {
286	(*fin) >> emin;
287	emin *= MeV;
288	} else if(line == "#emax(MeV)") {
289	(*fin) >> emax;
290	emax *= MeV;
291	} else if(line == "#nbin") {
292	(*fin) >> nbin;
293	} else if(line == "#cutG(mm)") {
294	(*fin) >> cutG;
295	cutG *= mm;
296	} else if(line == "#cutE(mm)") {
297	(*fin) >> cutE;
298	cutE *= mm;
299	} else if(line == "#range(mm)") {
300	(*fin) >> range;
301	range *= mm;
302	} else if(line == "#step(mm)") {
303	(*fin) >> theStep;
304	theStep *= mm;
305	} else if(line == "#material") {
306	(*fin) >> nameMat;
307	} else if(line == "#process") {
308	(*fin) >> nProcess;
309	} else if(line == "#domain") {
310	line1 = "";
311	(*fin) >> line1;
312	if(line1 == "lowenergy") {lowE = true;}
313	else {lowE = false;}
314	} else if(line == "#paw") {
315	usepaw = true;
316	(*fin) >> hFile;
317	} else if(line == "#run") {
318	break;
319	} else if(line == "#verbose") {
320	(*fin) >> verbose;
321	} else if(line == "#fluct") {
322	fluct = true;
323	} else if(line == "#exit") {
324	end = false;
325	break;
326	} else if(line == "#barkas") {
327	line1 = "";
328	(*fin) >> line1;
329	G4cout << line1 << G4endl;
330	if(line1 == "off") setBarkasOff = true;
331	if(line1 == "on") setBarkasOff = false;
332	} else if(line == "#nuclear") {
333	line1 = "";
334	(*fin) >> line1;
335	G4cout << line1 << G4endl;
336	if(line1 == "off") setNuclearOff = true;
337	if(line1 == "on") setNuclearOff = false;
338	}
339	} while(end);
340
341	if(!end) break;
342
343	G4cout << "###### Start new run # " << run << " #####" << G4endl;
344	if(nPart == 0) {
345	part = gamma;
346	} else if(nPart == 1) {
347	part = electron;
348	} else if(nPart == 2) {
349	part = positron;
350	} else if(nPart == 3) {
351	part = proton;
352	} else if(nPart == 4) {
353	part = antiproton;
354	} else if(nPart == 5) {
355	part = c12;
356	} else if(nPart == 6) {
357	part = ar40;
358	} else {
359	G4cout << "Particle #" << nPart
360	<< " is absent in the list of particles: Exit" << G4endl;
361	end = false;
362	break;
363	}
364	if(nProcess < 0 \|\| nProcess > 2) {
365	G4cout << "Process #" << nProcess
366	<< " is absent in the list of processes: Exit" << G4endl;
367	end = false;
368	break;
369	}
370
371
372	for (mat = 0; mat < nMaterials; mat++) {
373	material = ma[mat];
374	if(nameMat == material->GetName()) break;
375	}
376
377	G4cout << "The particle: " << part->GetParticleName() << G4endl;
378	G4cout << "The material: " << material->GetName() << G4endl;
379	G4cout << "The cut on e-:" << cutE/mm << " mm" << G4endl;
380	G4cout << "The cut on g: " << cutG/mm << " mm" << G4endl;
381	G4cout << "The step: " << theStep/mm << " mm" << G4endl;
382
383	// -------------------------------------------------------------------
384	// ---- HBOOK initialization
385
386	G4double emin10 = std::log10(emin/MeV);
387	G4double emax10 = std::log10(emax/MeV);
388	G4double bin = (emax10 - emin10) / (G4double)(nbin-1);
389
390	// Creating the analysis factory
391	std::auto_ptr< AIDA::IAnalysisFactory > af( AIDA_createAnalysisFactory() );
392
393	// Creating the tree factory
394	std::auto_ptr< AIDA::ITreeFactory > tf( af->createTreeFactory() );
395
396	// Creating a tree mapped to a new hbook file.
397	std::auto_ptr< AIDA::ITree > tree( tf->create( hFile,"hbook",false,false ) );
398	std::cout << "Tree store : " << tree->storeName() << std::endl;
399
400	// Creating a tuple factory, whose tuples will be handled by the tree
401	std::auto_ptr< AIDA::ITupleFactory > tpf( af->createTupleFactory( *tree ) );
402
403	AIDA::IHistogram1D* hist[4];
404	AIDA::ITuple* ntuple1 = 0;
405
406	if(usepaw) {
407
408	// ---- primary ntuple ------
409	// If using Anaphe HBOOK implementation, there is a limitation on the length of the
410	// variable names in a ntuple
411	ntuple1 = tpf->create( "100", "tuple", "float ekin, dedx" );
412
413
414	// Creating a histogram factory, whose histograms will be handled by the tree
415	std::auto_ptr< AIDA::IHistogramFactory > hf( af->createHistogramFactory( *tree ) );
416
417	// Creating an 1-dimensional histogram in the root directory of the tree
418
419	hist[0] = hf->createHistogram1D("11","Stopping power (MeVcm*2/g)",
420	nbin,emin10,emax10);
421	hist[1] = hf->createHistogram1D("12","Stopping power (MeV/mm)",
422	nbin,emin10,emax10);
423	hist[2] = hf->createHistogram1D("13","Step limit (mm)",
424	nbin,emin10,emax10);
425	hist[3] = hf->createHistogram1D("14","Number of secondaries",
426	nbin,emin10,emax10);
427
428
429	G4cout<< "Histograms is initialised" << G4endl;
430	}
431
432	gamma->SetCuts(cutG);
433	electron->SetCuts(cutE);
434	// positron->SetCuts(cutE);
435
436	// Processes - all new
437	G4bool success = false;
438
439	gmanager = new G4ProcessManager(gamma);
440	gmanager->AddDiscreteProcess(new G4LowEnergyPhotoElectric());
441	gmanager->AddDiscreteProcess(new G4LowEnergyCompton());
442	gmanager->AddDiscreteProcess(new G4LowEnergyGammaConversion());
443
444	if(lowE) {
445
446	elecManager = new G4ProcessManager(electron);
447	electron->SetProcessManager(elecManager);
448	eionle = new G4LowEnergyIonisation();
449	if(!fluct) eionle->SetEnlossFluc(false);
450	ebrle = new G4LowEnergyBremsstrahlung();
451	if(!fluct) ebrle->SetEnlossFluc(false);
452	elecManager->AddProcess(eionle);
453	elecManager->AddProcess(ebrle);
454	eionle->BuildPhysicsTable(*electron);
455	ebrle->BuildPhysicsTable(*electron);
456	if(nPart == 1) {
457	if(nProcess == 2) {
458	success = true;
459	}
460
461	} else if (nPart == 3) {
462	protManager = new G4ProcessManager(proton);
463	proton->SetProcessManager(protManager);
464	hionle = new G4hLowEnergyIonisation();
465	if(!fluct) hionle->SetEnlossFluc(false);
466	protManager->AddProcess(hionle);
467	if(setNuclearOff) hionle->SetNuclearStoppingOff();
468	if(!setNuclearOff) hionle->SetNuclearStoppingOn();
469	if(setBarkasOff) hionle->SetBarkasOff();
470	if(!setBarkasOff) hionle->SetBarkasOn();
471	hionle->SetVerboseLevel(verbose);
472	hionle->BuildPhysicsTable(*proton);
473	success = true;
474
475	} else if (nPart == 4) {
476	aprotManager = new G4ProcessManager(antiproton);
477	antiproton->SetProcessManager(aprotManager);
478	hionle = new G4hLowEnergyIonisation();
479	if(!fluct) hionle->SetEnlossFluc(false);
480	aprotManager->AddProcess(hionle);
481	if(setNuclearOff) hionle->SetNuclearStoppingOff();
482	if(!setNuclearOff) hionle->SetNuclearStoppingOn();
483	if(setBarkasOff) hionle->SetBarkasOff();
484	if(!setBarkasOff) hionle->SetBarkasOn();
485	hionle->SetVerboseLevel(verbose);
486	hionle->BuildPhysicsTable(*antiproton);
487	success = true;
488
489	} else if (nPart == 5 \|\| nPart == 6) {
490	protManager = new G4ProcessManager(proton);
491	proton->SetProcessManager(protManager);
492	hionle = new G4hLowEnergyIonisation();
493	if(!fluct) hionle->SetEnlossFluc(false);
494	protManager->AddProcess(hionle);
495	if(setNuclearOff) hionle->SetNuclearStoppingOff();
496	if(!setNuclearOff) hionle->SetNuclearStoppingOn();
497	if(setBarkasOff) hionle->SetBarkasOff();
498	if(!setBarkasOff) hionle->SetBarkasOn();
499	hionle->SetVerboseLevel(verbose);
500	hionle->BuildPhysicsTable(*proton);
501	ionManager = new G4ProcessManager(part);
502	part->SetProcessManager(ionManager);
503	ionle = new G4hLowEnergyIonisation();
504	if(!fluct) ionle->SetEnlossFluc(false);
505	if(setNuclearOff) ionle->SetNuclearStoppingOff();
506	if(!setNuclearOff) ionle->SetNuclearStoppingOn();
507	if(setBarkasOff) ionle->SetBarkasOff();
508	if(!setBarkasOff) ionle->SetBarkasOn();
509	ionManager->AddProcess(ionle);
510	ionle->SetVerboseLevel(verbose);
511	ionle->BuildPhysicsTable(*part);
512	success = true;
513	}
514
515	} else {
516
517	elecManager = new G4ProcessManager(electron);
518	electron->SetProcessManager(elecManager);
519	eionst = new G4eIonisation();
520	if(!fluct) eionst->SetEnlossFluc(false);
521	elecManager->AddProcess(eionst);
522	ebrst = new G4eBremsstrahlung();
523	if(!fluct) ebrst->SetEnlossFluc(false);
524	elecManager->AddProcess(ebrst);
525	eionst->BuildPhysicsTable(*electron);
526	ebrst->BuildPhysicsTable(*electron);
527	if(nPart == 1) {
528	if(nProcess == 2) {
529	success = true;
530	}
531	} else if(nPart == 2) {
532	positManager = new G4ProcessManager(positron);
533	positron->SetProcessManager(positManager);
534	if(nProcess == 0) {
535	eionst = new G4eIonisation();
536	if(!fluct) eionst->SetEnlossFluc(false);
537	positManager->AddProcess(eionst);
538	eionst->BuildPhysicsTable(*positron);
539	success = true;
540	} else if(nProcess == 1) {
541	ebrst = new G4eBremsstrahlung();
542	if(!fluct) ebrst->SetEnlossFluc(false);
543	positManager->AddProcess(ebrst);
544	ebrst->BuildPhysicsTable(*positron);
545	success = true;
546	} else if(nProcess == 2) {
547	eionst = new G4eIonisation();
548	ebrst = new G4eBremsstrahlung();
549	if(!fluct) eionst->SetEnlossFluc(false);
550	if(!fluct) ebrst->SetEnlossFluc(false);
551	positManager->AddProcess(eionst);
552	positManager->AddProcess(ebrst);
553	eionst->BuildPhysicsTable(*positron);
554	ebrst->BuildPhysicsTable(*positron);
555	success = true;
556	}
557
558	} else if (nPart == 3) {
559	protManager = new G4ProcessManager(proton);
560	proton->SetProcessManager(protManager);
561	hionst = new G4hIonisation();
562	if(!fluct) hionst->SetEnlossFluc(false);
563	protManager->AddProcess(hionst);
564	// hionst->SetVerboseLevel(verbose);
565	hionst->BuildPhysicsTable(*proton);
566	success = true;
567
568	} else if (nPart == 4) {
569	aprotManager = new G4ProcessManager(antiproton);
570	antiproton->SetProcessManager(aprotManager);
571	hionst = new G4hIonisation();
572	if(!fluct) hionst->SetEnlossFluc(false);
573	aprotManager->AddProcess(hionst);
574	hionst->SetVerboseLevel(verbose);
575	hionst->BuildPhysicsTable(*antiproton);
576	success = true;
577
578	} else if (nPart == 5 \|\| nPart == 6) {
579	protManager = new G4ProcessManager(proton);
580	proton->SetProcessManager(protManager);
581	hionst = new G4hIonisation();
582	if(!fluct) hionst->SetEnlossFluc(false);
583	protManager->AddProcess(hionst);
584	// hionst->SetVerboseLevel(verbose);
585	hionst->BuildPhysicsTable(*proton);
586	ionManager = new G4ProcessManager(part);
587	part->SetProcessManager(ionManager);
588	ionst = new G4hIonisation();
589	if(!fluct) ionst->SetEnlossFluc(false);
590	ionManager->AddProcess(ionst);
591	ionst->SetVerboseLevel(verbose);
592	ionst->BuildPhysicsTable(*part);
593	success = true;
594	}
595	}
596
597	if(success) G4cout << "Physics tables are built" << G4endl;
598	else G4cout << "Physics tables are not built!!!" << G4endl;
599
600	G4cout << "gCut(MeV)= " << gamma->GetEnergyThreshold(material)/MeV << G4endl;
601	G4cout << "eCut(MeV)= " << electron->GetEnergyThreshold(material)/MeV << G4endl;
602
603
604	// Create a DynamicParticle
605
606	G4ParticleMomentum gDir(initX,initY,initZ);
607	G4double gEnergy = emax;
608	G4DynamicParticle dParticle(part,gDir,gEnergy);
609
610	// Track
611	G4ThreeVector aPosition(0.,0.,0.);
612	G4double aTime = 0. ;
613
614	G4Track* gTrack;
615	gTrack = new G4Track(&dParticle,aTime,aPosition);
616
617	// Step
618
619	G4Step* step;
620	step = new G4Step();
621	step->SetTrack(gTrack);
622
623	G4StepPoint aPoint, bPoint;
624	aPoint = new G4StepPoint();
625	aPoint->SetPosition(aPosition);
626	aPoint->SetMaterial(material);
627	G4double safety = 10000.*cm;
628	aPoint->SetSafety(safety);
629	step->SetPreStepPoint(aPoint);
630
631	bPoint = aPoint;
632	G4ThreeVector bPosition(0.,0.,theStep);
633	bPoint->SetPosition(bPosition);
634	step->SetPostStepPoint(bPoint);
635	step->SetStepLength(theStep);
636
637	if(!fluct) {
638	nstatf = 1;
639	xstatf = 1.0;
640	}
641
642	G4Timer* timer = new G4Timer();
643	timer->Start();
644
645	G4double le = emin10 - bin;
646
647	for (G4int iter=0; iter<nbin; iter++) {
648
649	le += bin;
650	G4double e = std::pow(10.0,le) * MeV;
651	gTrack->SetStep(step);
652	gTrack->SetKineticEnergy(e);
653
654	for (G4int jj=0; jj<nstatf; jj++) {
655
656	G4double x = 0.0;
657	G4VParticleChange* aChange = 0;
658
659	if(lowE) {
660
661	if(nPart == 1) {
662	if(nProcess == 0) {
663	x = eionle->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
664	aChange = eionle->AlongStepDoIt(gTrack,step);
665
666	} else if(nProcess == 1) {
667	x = ebrle->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
668	aChange = ebrle->AlongStepDoIt(gTrack,step);
669	} else if(nProcess == 2) {
670	x = eionle->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
671	aChange = eionle->AlongStepDoIt(gTrack,step);
672	}
673
674	} else if (nPart == 3 ) {
675	x = hionle->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
676	aChange = hionle->AlongStepDoIt(gTrack,step);
677
678	} else if (nPart == 4) {
679	x = hionle->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
680	aChange = hionle->AlongStepDoIt(gTrack,step);
681
682	} else if (nPart == 5 \|\| nPart == 6) {
683	x = ionle->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
684	aChange = ionle->AlongStepDoIt(gTrack,step);
685
686	}
687
688	} else {
689
690
691	if(nPart == 1) {
692	if(nProcess == 0) {
693	x = eionst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
694	aChange = eionst->AlongStepDoIt(gTrack,step);
695
696	} else if(nProcess == 1) {
697	x = ebrst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
698	aChange = ebrst->AlongStepDoIt(gTrack,step);
699	} else if(nProcess == 2) {
700	x = eionst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
701	aChange = eionst->AlongStepDoIt(gTrack,step);
702	}
703
704	} else if(nPart == 2) {
705	if(nProcess == 0) {
706	x = eionst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
707	aChange = eionst->AlongStepDoIt(gTrack,step);
708
709	} else if(nProcess == 1) {
710	x = ebrst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
711	aChange = ebrst->AlongStepDoIt(gTrack,step);
712	} else if(nProcess == 2) {
713	x = eionst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
714	aChange = eionst->AlongStepDoIt(gTrack,step);
715	}
716
717	} else if (nPart == 3) {
718	x = hionst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
719	aChange = hionst->AlongStepDoIt(gTrack,step);
720
721	} else if (nPart == 4) {
722	x = hionst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
723	aChange = hionst->AlongStepDoIt(gTrack,step);
724
725	} else if (nPart == 5 \|\| nPart == 6) {
726	x = ionst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
727	aChange = ionst->AlongStepDoIt(gTrack,step);
728	}
729	}
730
731	G4double delx = theStep;
732	G4double de = aChange->GetLocalEnergyDeposit();
733	G4int n = aChange->GetNumberOfSecondaries();
734
735
736	//G4cout << " de(MeV) = " << de/MeV << " n= " << n << G4endl;
737
738	if(n > 0) {
739	for(G4int i=0; i<n; i++) {
740	de += (aChange->GetSecondary(i))->GetKineticEnergy();
741	if(verbose) {
742	G4cout << "add "
743	<< ((aChange->GetSecondary(i))->GetKineticEnergy())/eV
744	<< " eV" << G4endl;
745	}
746	}
747	}
748	G4double st = de/(delx*(material->GetDensity()));
749	st = gram/(cmcm*MeV);
750	G4double s = demm/(delxMeV);
751
752	G4cout << "E(MeV)= " << e/MeV << " dE/dx(MeV/mm)= " << s << G4endl;
753
754	if(verbose) {
755	G4cout << "Iteration = " << iter
756	<< " E = " << e/MeV << " MeV; StepLimit= "
757	<< x/mm << " mm; de= "
758	<< de/eV << " eV; dE/dx= "
759	<< st << " MeV*cm^2/g" << G4endl;
760	}
761
762	if(x > 1000.0meter) x = 1000.0meter;
763
764	if (usepaw) {
765	float st10 = -5.0;
766	if(st > 1.e-5) st10 = (float)log10(st);
767	if(verbose>1) {
768	G4cout << " de(MeV) = " << de/MeV
769	<< G4endl;
770	G4cout << " n1= " << ntuple1->findColumn("ekin")
771	<< " n2= " << ntuple1->findColumn("dedx")
772	<< G4endl;
773	}
774	ntuple1->fill( ntuple1->findColumn("ekin"), (float)le);
775	ntuple1->fill( ntuple1->findColumn("dedx"), st10);
776	ntuple1->addRow();
777	// G4cout << "ntuple is filled " << G4endl;
778
779	hist[0]->fill(le,st*xstatf);
780	hist[1]->fill(le,s*xstatf);
781	hist[2]->fill(le,x*xstatf/mm);
782	hist[3]->fill(le,xstatf*(G4double)n);
783	}
784	}
785	}
786
787	timer->Stop();
788	G4cout << " " << *timer << G4endl;
789	delete timer;
790
791	// Committing the transaction with the tree
792	if(usepaw) {
793	std::cout << "Committing..." << std::endl;
794	tree->commit();
795	std::cout << "Closing the tree..." << std::endl;
796	tree->close();
797	}
798	G4cout << "###### End of run # " << run << " ######" << G4endl;
799
800
801	} while(end);
802	G4cout << "###### End of test #####" << G4endl;
803	}
804
805	#include "hTest/src/G4IonC12.cc"
806	#include "hTest/src/G4IonAr40.cc"

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

Download in other formats: