Context Navigation

source: trunk/source/processes/electromagnetic/lowenergy/test/G4ComplexTest.cc @ 1199

Last change on this file since 1199 was 1199, checked in by garnier, 15 years ago
nvx fichiers dans CVS
File size: 24.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	// $Id: G4ComplexTest.cc,v 1.21 2006/06/29 19:43:49 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: G4ComplexTest
36	// This test provide AlongStepDoIt and PostStepDoIt
37	// tests for electromagnetic processes. The input
38	// data have to be describe in ASCII file
39	//
40	// Author: V.Ivanchenko on base of Maria Grazia Pia tests
41	//
42	// Creation date: 8 May 2001
43	//
44	// Modifications:
45	//
46	// -------------------------------------------------------------------
47
48	#include "globals.hh"
49	#include "G4ios.hh"
50	#include <fstream>
51	#include <iomanip>
52
53	#include "G4ProductionCuts.hh"
54	#include "G4ProductionCutsTable.hh"
55	#include "G4Material.hh"
56	#include "G4MaterialCutsCouple.hh"
57	#include "G4VContinuousDiscreteProcess.hh"
58	#include "G4ProcessManager.hh"
59
60	#include "G4LowEnergyIonisation.hh"
61	#include "G4LowEnergyBremsstrahlung.hh"
62	#include "G4LowEnergyCompton.hh"
63	#include "G4LowEnergyGammaConversion.hh"
64	#include "G4LowEnergyPhotoElectric.hh"
65	#include "G4LowEnergyRayleigh.hh"
66
67	#include "G4hLowEnergyIonisation.hh"
68
69	#include "G4eIonisation.hh"
70	#include "G4eBremsstrahlung.hh"
71	#include "G4ComptonScattering.hh"
72	#include "G4GammaConversion.hh"
73	#include "G4PhotoElectricEffect.hh"
74
75	#include "G4eplusAnnihilation.hh"
76
77	#include "G4MuIonisation.hh"
78	#include "G4MuBremsstrahlung.hh"
79	#include "G4MuPairProduction.hh"
80
81	#include "G4hIonisation.hh"
82
83	#include "G4MultipleScattering.hh"
84
85	#include "G4EnergyLossTables.hh"
86	#include "G4ParticleChange.hh"
87	#include "G4ParticleChange.hh"
88	#include "G4DynamicParticle.hh"
89	#include "G4AntiProton.hh"
90	#include "G4Proton.hh"
91	#include "G4Electron.hh"
92	#include "G4Positron.hh"
93	#include "G4Gamma.hh"
94
95	#include "G4Box.hh"
96	#include "G4PVPlacement.hh"
97	#include "G4VPhysicalVolume.hh"
98	#include "G4LogicalVolume.hh"
99	#include "G4RunManager.hh"
100
101	#include "G4Step.hh"
102	#include "G4GRSVolume.hh"
103
104	#include "G4UnitsTable.hh"
105
106	// New Histogramming (from AIDA and Anaphe):
107	#include <memory> // for the auto_ptr(T>
108
109	#include "AIDA/AIDA.h"
110
111
112	#include "G4Timer.hh"
113
114	int main(int argc,char** argv)
115	{
116
117	// -------------------------------------------------------------------
118	// Setup
119
120	G4int nEvt = 100;
121	G4int nPart =-1;
122	G4String nameMat = "Si";
123	G4int nProcess = 0;
124	G4bool usepaw = false;
125	G4bool postDo = true;
126	G4bool lowE = true;
127	G4int verbose = 0;
128	G4double gEnergy = 0.1*MeV;
129	G4String hFile = "";
130	G4double theStep = 1.0*micrometer;
131	G4double range = 1.0*micrometer;
132	G4double cutG = 1.0*micrometer;
133	G4double cutE = 1.0*micrometer;
134	G4Material* material = 0;
135	G4String name[6] = {"Ionisation", "Bremsstrahlung", "Compton",
136	"GammaConversion", "PhotoElectric", "Raylaigh"};
137
138
139	G4cout.setf( ios::scientific, ios::floatfield );
140
141	// -------------------------------------------------------------------
142	// Control on input
143
144	if(argc < 2) {
145	G4cout << "Input file is not specified! Exit" << G4endl;
146	exit(1);
147	}
148
149	ifstream* fin = new ifstream();
150	string fname = argv[1];
151	fin->open(fname.c_str());
152	if( !fin->is_open()) {
153	G4cout << "Input file <" << fname << "> does not exist! Exit" << G4endl;
154	exit(1);
155	}
156
157	// -------------------------------------------------------------------
158	//--------- Materials definition ---------
159
160	G4Material* m;
161	m = new G4Material("Be", 4., 9.01g/mole, 1.848g/cm3);
162	m = new G4Material("Graphite",6., 12.00g/mole, 2.265g/cm3 );
163	m = new G4Material("Al", 13., 26.98g/mole, 2.7 g/cm3);
164	m = new G4Material("Si", 14., 28.055g/mole, 2.33g/cm3);
165	m = new G4Material("LAr", 18., 39.95g/mole, 1.393g/cm3);
166	m = new G4Material("Fe", 26., 55.85g/mole, 7.87g/cm3);
167	m = new G4Material("Cu", 29., 63.55g/mole, 8.96g/cm3);
168	m = new G4Material("W", 74., 183.85g/mole, 19.30g/cm3);
169	m = new G4Material("Pb", 82., 207.19g/mole, 11.35g/cm3);
170	m = new G4Material("U", 92., 238.03g/mole, 18.95g/cm3);
171
172	G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
173	G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
174	G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
175	G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
176	G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole);
177
178	m = new G4Material("O2", 8., 16.00g/mole, 1.1g/cm3);
179
180	m = new G4Material ("Water" , 1.*g/cm3, 2);
181	m->AddElement(H,2);
182	m->AddElement(O,1);
183
184	m = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
185	m->AddElement(H,6);
186	m->AddElement(C,2);
187
188	m = new G4Material ("CsI" , 4.53*g/cm3, 2);
189	m->AddElement(Cs,1);
190	m->AddElement(I,1);
191
192	const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
193
194	G4int nMaterials = G4Material::GetNumberOfMaterials();
195
196	G4cout << "Available materials are: " << G4endl;
197	G4int mat;
198	for (mat = 0; mat < nMaterials; mat++) {
199	G4cout << mat << ") " << (*theMaterialTable)[mat]->GetName() << G4endl;
200	}
201	material = (*theMaterialTable)[0];
202
203	G4cout << "Available processes are: " << G4endl;
204	for (mat = 0; mat < 6; mat++) {
205	G4cout << mat << ") " << name[mat] << G4endl;
206	}
207
208	// Particle definitions
209
210	G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
211	G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
212	G4ParticleDefinition* positron = G4Positron::PositronDefinition();
213	G4ParticleDefinition* proton = G4Proton::ProtonDefinition();
214	G4ParticleDefinition* antiproton = G4AntiProton::AntiProtonDefinition();
215	G4ParticleDefinition* part = gamma;
216
217	// Geometry
218
219	G4double initX = 0.;
220	G4double initY = 0.;
221	G4double initZ = 1.;
222	G4double dimX = 100.0*cm;
223	G4double dimY = 100.0*cm;
224	G4double dimZ = 100.0*cm;
225
226	G4Box* sFrame = new G4Box ("Box",dimX, dimY, dimZ);
227	G4LogicalVolume* lFrame = new G4LogicalVolume(sFrame,material,"Box",0,0,0);
228	G4PVPlacement* pFrame = new G4PVPlacement(0,G4ThreeVector(),"Box",
229	lFrame,0,false,0);
230	G4RunManager* rm = new G4RunManager();
231	G4cout << "World is defined " << G4endl;
232	rm->GeometryHasBeenModified();
233	rm->DefineWorldVolume(pFrame);
234
235	// -------------------------------------------------------------------
236	// ---- Read input file
237	G4cout << "Available commands are: " << G4endl;
238	G4cout << "#events" << G4endl;
239	G4cout << "#particle" << G4endl;
240	G4cout << "#energy(MeV)" << G4endl;
241	G4cout << "#cutG(mm)" << G4endl;
242	G4cout << "#cutE(mm)" << G4endl;
243	G4cout << "#range(mm)" << G4endl;
244	G4cout << "#step(mm)" << G4endl;
245	G4cout << "#material" << G4endl;
246	G4cout << "#process" << G4endl;
247	G4cout << "#domain" << G4endl;
248	G4cout << "#test" << G4endl;
249	G4cout << "#paw" << G4endl;
250	G4cout << "#verbose" << G4endl;
251	G4cout << "#run" << G4endl;
252	G4cout << "#exit" << G4endl;
253	G4cout << pFrame << G4endl;
254
255	G4ProcessManager elecManager, positManager, *gammaManager,
256	protManager, aprotManager;
257
258	elecManager = new G4ProcessManager(electron);
259	electron->SetProcessManager(elecManager);
260
261	positManager = new G4ProcessManager(positron);
262	positron->SetProcessManager(positManager);
263
264	gammaManager = new G4ProcessManager(gamma);
265	gamma->SetProcessManager(gammaManager);
266
267	protManager = new G4ProcessManager(proton);
268	proton->SetProcessManager(protManager);
269
270	aprotManager = new G4ProcessManager(antiproton);
271	antiproton->SetProcessManager(aprotManager);
272
273	G4eIonisation* elecSTion = 0;
274	G4eBremsstrahlung* elecSTbr = 0;
275	G4LowEnergyIonisation* elecLEion = 0;
276	G4LowEnergyBremsstrahlung* elecLEbr = 0;
277	G4bool ionis = true;
278
279	string line, line1;
280	G4bool end = true;
281	for(G4int run=0; run<100; run++) {
282	do {
283	(*fin) >> line;
284	G4cout << "Next line " << line << G4endl;
285	if(line == "#events") {
286	(*fin) >> nEvt;
287	if(nEvt < 1) nEvt = 1;
288	} else if(line == "#particle") {
289	(*fin) >> nPart;
290	} else if(line == "#energy(MeV)") {
291	(*fin) >> gEnergy;
292	gEnergy *= MeV;
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	(*fin) >> line1;
311	if(line1 == "lowenergy") {lowE = true;}
312	else {lowE = false;}
313	} else if(line == "#test") {
314	(*fin) >> line1;
315	if(line1 == "PostStep") {postDo = true;}
316	else {postDo = false;}
317	} else if(line == "#paw") {
318	usepaw = true;
319	(*fin) >> hFile;
320	} else if(line == "#run") {
321	break;
322	} else if(line == "#verbose") {
323	(*fin) >> verbose;
324	} else if(line == "#exit") {
325	end = false;
326	break;
327	}
328	} while(end);
329
330	if(!end) break;
331
332	G4cout << "###### Start new run # " << run << " #####" << G4endl;
333	if(nPart == 0) {
334	part = gamma;
335	} else if(nPart == 1) {
336	part = electron;
337	} else if(nPart == 2) {
338	part = positron;
339	} else if(nPart == 3) {
340	part = proton;
341	} else if(nPart == 4) {
342	part = antiproton;
343	} else {
344	G4cout << "Particle #" << nPart
345	<< " is absent in the list of particles: Exit" << G4endl;
346	break;
347	}
348	if(nProcess < 0 \|\| nProcess > 5) {
349	G4cout << "Process #" << nProcess
350	<< " is absent in the list of processes: Exit" << G4endl;
351	break;
352	}
353
354	for (mat = 0; mat < nMaterials; mat++) {
355	material = (*theMaterialTable)[mat];
356	if(nameMat == material->GetName()) break;
357	}
358
359	G4cout << "The particle: " << part->GetParticleName() << G4endl;
360	G4cout << "The energy: " << gEnergy/MeV << " MeV" << G4endl;
361	G4cout << "The material: " << material->GetName() << G4endl;
362	G4cout << "The cut on e-:" << cutE/mm << " mm" << G4endl;
363	G4cout << "The cut on g: " << cutG/mm << " mm" << G4endl;
364	G4cout << "The step: " << theStep/mm << " mm" << G4endl;
365	if(postDo && lowE) {
366	G4cout << "Test of PostStepDoIt for " << name[nProcess]
367	<< " for lowenergy" << G4endl;
368	} else if(postDo && !lowE) {
369	G4cout << "Test of PostStepDoIt for " << name[nProcess]
370	<< " for standard" << G4endl;
371	} else if(!postDo && !lowE) {
372	G4cout << "Test of AlongStepDoIt for " << name[nProcess]
373	<< " for standard" << G4endl;
374	} else if(!postDo && lowE) {
375	G4cout << "Test of AlongStepDoIt for " << name[nProcess]
376	<< " for lowenergy" << G4endl;
377	}
378
379	// -------------------------------------------------------------------
380	// ---- HBOOK initialization
381
382	// Creating the analysis factory
383	std::auto_ptr< AIDA::IAnalysisFactory > af( AIDA_createAnalysisFactory() );
384
385	// Creating the tree factory
386	std::auto_ptr< AIDA::ITreeFactory > tf( af->createTreeFactory() );
387
388	// Creating a tree mapped to a new hbook file.
389	std::auto_ptr< AIDA::ITree > tree( tf->create( hFile,"hbook" ,false,false ) );
390	G4cout << "Tree store : " << tree->storeName() << G4endl;
391
392	// Creating a tuple factory, whose tuples will be handled by the tree
393	// std::auto_ptr< AIDA::ITupleFactory > tpf( af->createTupleFactory( *tree ) );
394
395	AIDA::IHistogram1D* hist[4];
396	//AIDA::ITuple* ntuple1 = 0;
397	//AIDA::ITuple* ntuple2 = 0;
398
399	if(usepaw) {
400
401	// ---- primary ntuple ------
402	// If using Anaphe HBOOK implementation, there is a limitation on the length of the
403	// variable names in a ntuple
404	//ntuple1 = tpf->create( "100", "Primary", "float ekin, dedx" );
405	//ntuple2 = tpf->create( "101", "Secondary", "float ekin, dedx" );
406
407
408	// Creating a histogram factory, whose histograms will be handled by the tree
409	std::auto_ptr< AIDA::IHistogramFactory > hf( af->createHistogramFactory( *tree ) );
410
411	// Creating an 1-dimensional histogram in the root directory of the tree
412
413	hist[0] = hf->createHistogram1D("11","Kinetic Energy (T/T0)", 50,0.,1.0);
414	hist[1] = hf->createHistogram1D("12","Momentum (MeV/c)", 50,0.,gEnergy*0.1/MeV);
415	hist[2] = hf->createHistogram1D("13","Number of secondaries", 20,-0.5,19.5);
416	hist[3] = hf->createHistogram1D("14","Energy deposition (MeV)", 50,0.,gEnergy*0.1/MeV);
417
418	G4cout<< "Histograms is initialised" << G4endl;
419	}
420
421	G4Timer* timer = new G4Timer();
422	timer->Start();
423	G4ProductionCutsTable* cutsTable = G4ProductionCutsTable::GetProductionCutsTable();
424
425	G4ProductionCuts* cuts = cutsTable->GetDefaultProductionCuts();
426	cuts->SetProductionCut(cutG, 0);
427	cuts->SetProductionCut(cutE, 1);
428	cuts->SetProductionCut(cutE, 2);
429	G4cout << "Cuts are defined " << G4endl;
430
431
432	cutsTable->UpdateCoupleTable();
433	(G4ProductionCutsTable::GetProductionCutsTable())->DumpCouples();
434	const G4MaterialCutsCouple* theCouple = cutsTable->GetMaterialCutsCouple(material,cuts);
435	// Processes
436
437	G4VDiscreteProcess* dProcess;
438	G4VContinuousDiscreteProcess* cdProcess;
439	dProcess = 0;
440	cdProcess = 0;
441
442	G4cout << "Start BuildPhysicsTable" << G4endl;;
443
444	if(lowE) {
445	if(ionis) {
446	elecLEion = new G4LowEnergyIonisation();
447	elecLEbr = new G4LowEnergyBremsstrahlung();
448	elecManager->AddProcess(elecLEion);
449	elecManager->AddProcess(elecLEbr);
450	elecLEion->BuildPhysicsTable(*electron);
451	elecLEbr->BuildPhysicsTable(*electron);
452	ionis = false;
453	}
454	if(nPart == 0) {
455	if(nProcess == 2) dProcess = new G4LowEnergyCompton();
456	if(nProcess == 3) dProcess = new G4LowEnergyGammaConversion();
457	if(nProcess == 4) dProcess = new G4LowEnergyPhotoElectric();
458	if(nProcess == 5) dProcess = new G4LowEnergyRayleigh();
459	if(dProcess) {
460	gammaManager->AddProcess(dProcess);
461	dProcess->BuildPhysicsTable(*gamma);
462	}
463	} else if(nPart == 1) {
464	if(nProcess == 0) cdProcess = elecLEion;
465	if(nProcess == 1) cdProcess = elecLEbr;
466	} else if(nPart == 3) {
467	if(nProcess == 0) {
468	cdProcess = new G4hLowEnergyIonisation();
469	protManager->AddProcess(cdProcess);
470	cdProcess->BuildPhysicsTable(*proton);
471	}
472	} else if(nPart == 4) {
473	if(nProcess == 0) {
474	cdProcess = new G4hLowEnergyIonisation();
475	aprotManager->AddProcess(cdProcess);
476	cdProcess->BuildPhysicsTable(*antiproton);
477	}
478	}
479
480	} else {
481	if(ionis) {
482	elecSTion = new G4eIonisation();
483	elecSTbr = new G4eBremsstrahlung();
484	elecManager->AddProcess(elecSTion);
485	elecManager->AddProcess(elecSTbr);
486	elecSTion->BuildPhysicsTable(*electron);
487	elecSTbr->BuildPhysicsTable(*electron);
488	ionis = false;
489	}
490	if(nPart == 0) {
491	if(nProcess == 2) dProcess = new G4ComptonScattering();
492	if(nProcess == 3) dProcess = new G4GammaConversion();
493	if(nProcess == 4) dProcess = new G4PhotoElectricEffect();
494	if(dProcess) {
495	gammaManager->AddProcess(dProcess);
496	dProcess->BuildPhysicsTable(*gamma);
497	}
498	} else if(nPart == 1) {
499	if(nProcess == 0) cdProcess = elecSTion;
500	if(nProcess == 1) cdProcess = elecSTbr;
501	} else if(nPart == 2) {
502	G4eIonisation* pSTion = new G4eIonisation();
503	G4eBremsstrahlung* pSTbr = new G4eBremsstrahlung();
504	if(nProcess == 0) cdProcess = pSTion;
505	if(nProcess == 1) cdProcess = pSTbr;
506	} else if(nPart == 3) {
507	if(nProcess == 0) {
508	cdProcess = new G4hIonisation();
509	protManager->AddProcess(cdProcess);
510	cdProcess->BuildPhysicsTable(*proton);
511	}
512	} else if(nPart == 4) {
513	if(nProcess == 0) {
514	cdProcess = new G4hIonisation();
515	aprotManager->AddProcess(cdProcess);
516	cdProcess->BuildPhysicsTable(*antiproton);
517	}
518	}
519	}
520
521	G4cout << "Physics tables are built" << G4endl;;
522
523	timer->Stop();
524	G4cout << " " << *timer << G4endl;
525	delete timer;
526
527
528	// Control on processes
529	if(postDo && !dProcess && !cdProcess) {
530	G4cout << "Discret Process is not found out! Exit" << G4endl;
531	break;
532	}
533	if(!postDo && !cdProcess) {
534	G4cout << "Continues Discret Process is not found out! Exit" << G4endl;
535	break;
536	}
537
538	// Create a DynamicParticle
539
540	G4ParticleMomentum gDir(initX,initY,initZ);
541	G4DynamicParticle dParticle(part,gDir,gEnergy);
542
543	// Track
544	G4ThreeVector aPosition(0.,0.,0.);
545	G4double aTime = 0. ;
546
547	G4Track* gTrack;
548	gTrack = new G4Track(&dParticle,aTime,aPosition);
549
550	// Step
551
552	G4Step* step;
553	step = new G4Step();
554	step->SetTrack(gTrack);
555
556	G4StepPoint aPoint, bPoint;
557	aPoint = new G4StepPoint();
558	aPoint->SetPosition(aPosition);
559	aPoint->SetMaterial(material);
560	aPoint->SetMaterialCutsCouple(theCouple);
561	G4double safety = 10000.*cm;
562	aPoint->SetSafety(safety);
563	step->SetPreStepPoint(aPoint);
564
565	bPoint = aPoint;
566	G4ThreeVector bPosition(0.,0.,theStep);
567	bPoint->SetPosition(bPosition);
568	step->SetPostStepPoint(bPoint);
569	step->SetStepLength(theStep);
570
571	// --------- Test the DoIt
572	G4int nElectrons = 0;
573	G4int nPositrons = 0;
574	G4int nPhotons = 0;
575	G4double rmax = 0.0;
576	G4double de = 0.0;
577	G4double de2 = 0.0;
578
579	G4cout << "dProcess= " << dProcess << " cdProcess= " << cdProcess;
580	if(dProcess) G4cout << " name= " << dProcess->GetProcessName() << G4endl;
581	if(cdProcess) G4cout << " name= " << cdProcess->GetProcessName() << G4endl;
582
583
584	timer = new G4Timer();
585	timer->Start();
586
587	for (G4int iter=0; iter<nEvt; iter++) {
588
589	gTrack->SetStep(step);
590
591	if(verbose) {
592	G4cout << "Iteration = " << iter
593	<< " - Step Length = "
594	<< step->GetStepLength()/mm << " mm "
595	<< G4endl;
596	}
597
598	G4VParticleChange* dummy = 0;
599	if(postDo) {
600	if(dProcess) dummy = dProcess->PostStepDoIt(gTrack, step);
601	if(cdProcess) {
602	dummy = cdProcess->PostStepDoIt(gTrack, step);
603	}
604	} else {
605	dummy = cdProcess->AlongStepDoIt(gTrack, step);
606	}
607	G4ParticleChange* particleChange = (G4ParticleChange*) dummy;
608
609	// Primary physical quantities
610
611	G4double energyChange = particleChange->GetEnergyChange();
612	G4double deltaE = gEnergy - energyChange ;
613	G4double dedx = deltaE / (step->GetStepLength());
614
615	G4ThreeVector change = particleChange->CalcMomentum(energyChange,
616	*(particleChange->GetMomentumChange()),
617	part->GetPDGMass());
618
619	G4double pxChange = change.x();
620	G4double pyChange = change.y();
621	G4double pzChange = change.z();
622
623	if(verbose) {
624	G4cout << "---- Primary after the step ---- " << G4endl;
625	G4cout << "---- Energy: " << energyChange/MeV << " MeV, "
626	<< "(px,py,pz): ("
627	<< pxChange/MeV << ","
628	<< pyChange/MeV << ","
629	<< pzChange/MeV << ") MeV"
630	<< G4endl;
631	G4cout << "---- Energy loss (dE) = " << deltaE/keV << " keV;"
632	<< "Stopping power (dE/dx)=" << dedx*mm/keV << " keV/mm"
633	<< "; rmax(mm)= " << rmax << G4endl;
634	}
635
636	// Primary
637
638	G4int nsec = particleChange->GetNumberOfSecondaries();
639	/*
640	if(ntuple1) {
641	ntuple1->column("epri", gEnergy/MeV);
642	ntuple1->column("efin", energyChange/MeV);
643	ntuple1->column("dedx", dedx*mm/MeV);
644	ntuple1->column("nsec", nsec);
645	ntuple1->column("nele", nElectrons);
646	ntuple1->column("npho", nPhotons);
647	ntuple1->dumpData();
648	}
649	*/
650	de += deltaE;
651	de2 += deltaE*deltaE;
652
653	// Secondaries physical quantities
654
655	if(usepaw) {
656	hist[2]->fill((float)nsec, (float)1.0);
657	hist[3]->fill((float)(particleChange->GetLocalEnergyDeposit()/MeV), (float)1.0);
658	}
659
660	for (G4int i = 0; i<nsec; i++) {
661	// The following two items should be filled per event, not
662	// per secondary; filled here just for convenience, to avoid
663	// complicated logic to dump ntuple when there are no secondaries
664
665	G4Track* finalParticle = particleChange->GetSecondary(i) ;
666
667	G4double e = finalParticle->GetTotalEnergy();
668	G4double eKin = finalParticle->GetKineticEnergy();
669	G4double px = (finalParticle->GetMomentum()).x();
670	G4double py = (finalParticle->GetMomentum()).y();
671	G4double pz = (finalParticle->GetMomentum()).z();
672	//G4double theta= (finalParticle->GetMomentum()).theta();
673	G4double p = std::sqrt(pxpx + pypy + pz*pz);
674
675	if (eKin > gEnergy) {
676	G4cout << "WARNING: eFinal > eInit in event #" << iter << G4endl;
677	}
678
679	G4String partName = finalParticle->GetDefinition()->GetParticleName();
680	if(verbose) {
681	G4cout << "==== Final "
682	<< partName << " "
683	<< "E= " << e/MeV << " MeV, "
684	<< "eKin: " << eKin/MeV << " MeV, "
685	<< "(px,py,pz): ("
686	<< px/MeV << ","
687	<< py/MeV << ","
688	<< pz/MeV << ") MeV,"
689	<< " p= " << p << " MeV"
690	<< G4endl;
691	}
692
693	if(usepaw) {
694	hist[0]->fill((float)(eKin/gEnergy), (float)1.0);
695	hist[1]->fill((float)(p/MeV), (float)1.0);
696	}
697
698	G4int partType = 0;
699	if (partName == "e-") {
700	partType = 1;
701	nElectrons++;
702
703	} else if (partName == "e+") {
704	partType = 2;
705	nPositrons++;
706
707	} else if (partName == "gamma") {
708	partType = 0;
709	nPhotons++;
710	}
711
712	// Fill the secondaries ntuple
713	/*
714	if(ntuple2) {
715	ntuple2->column("eprimary",gEnergy);
716	ntuple2->column("px", px);
717	ntuple2->column("py", py);
718	ntuple2->column("pz", pz);
719	ntuple2->column("p", p);
720	ntuple2->column("e", e);
721	ntuple2->column("theta", theta);
722	ntuple2->column("ekin", eKin);
723	ntuple2->column("type", partType);
724	ntuple2->dumpData();
725	}
726	*/
727	delete particleChange->GetSecondary(i);
728	}
729
730	particleChange->Clear();
731
732	}
733	G4cout << "###### Statistics:" << G4endl;
734	G4cout << "Average number of secondary electrons= "
735	<< (G4double)nElectrons/(G4double)nEvt << G4endl;
736	G4cout << "Average number of secondary positrons= "
737	<< (G4double)nPositrons/(G4double)nEvt << G4endl;
738	G4cout << "Average number of secondary photons= "
739	<< (G4double)nPhotons/(G4double)nEvt << G4endl;
740	G4double x = de/(G4double)nEvt;
741	G4double y = de2/(G4double)nEvt - x*x;
742	if(0.0 < y) y = std::sqrt(y);
743	G4cout << "Average energy deposition(MeV)= "
744	<< x/MeV << " +- " << y/MeV << G4endl;
745
746	timer->Stop();
747	G4cout << " " << *timer << G4endl;
748	delete timer;
749
750	if(usepaw) {
751	tree->commit();
752	std::cout << "Closing the tree..." << std::endl;
753	tree->close();
754	G4cout << "# hbook is writed" << G4endl;
755	}
756
757	G4cout << "###### End of run # " << run << " ######" << G4endl;
758
759	} while(end);
760	G4cout << "###### End of test #####" << G4endl;
761	}
762
763
764
765
766
767
768
769
770
771
772
773

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

Download in other formats: