Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

G4LowEnergyTest.cc @ 1350

Last change on this file since 1350 was 1350, checked in by garnier, 13 years ago
update to last version 4.9.4
File size: 21.5 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: G4LowEnergyTest.cc,v 1.9 2006/06/29 19:44:10 gunter Exp $
28	// GEANT4 tag $Name: geant4-09-04-ref-00 $
29	//
30	// KaonMinusAtRestTest.cc
31	// -------------------------------------------------------------------
32	// GEANT 4 class file --- Copyright CERN 1998
33	// CERN Geneva Switzerland
34	//
35	//
36	// File name: G4LowEnergyTest
37	//
38	// Author: Christian Voelcker (from M. Maire)
39	//
40	// Creation date: ?
41	//
42	// Modifications:
43	//
44	// Alessandra Forti march 1999 added ntuples
45	// Alessandra Forti march 1999 added processes
46	// Alessandra Forti march 1999 added cut in energy
47	// -------------------------------------------------------------------
48
49	#include "G4ios.hh"
50	#include <fstream>
51	#include <iomanip>
52
53	#include "G4Material.hh"
54
55	#include "G4ProcessManager.hh"
56	#include "G4LowEnergyPhotoElectric.hh"
57	#include "G4LowEnergyCompton.hh"
58	#include "G4LowEnergyRayleigh.hh"
59	#include "G4LowEnergyGammaConversion.hh"
60	#include "G4LowEnergyBremsstrahlung.hh"
61	#include "G4LowEnergyIonisation.hh"
62	#include "G4eIonisation.hh"
63	#include "G4hIonisation.hh"
64
65	#include "G4VParticleChange.hh"
66	#include "G4ParticleChange.hh"
67	#include "G4DynamicParticle.hh"
68	#include "G4Electron.hh"
69	#include "G4Positron.hh"
70	#include "G4Gamma.hh"
71
72	#include "G4Box.hh"
73	#include "G4PVPlacement.hh"
74
75	#include "G4Step.hh"
76	#include "G4GRSVolume.hh"
77
78	#include "CLHEP/Hist/TupleManager.h"
79	#include "CLHEP/Hist/HBookFile.h"
80	#include "CLHEP/Hist/Histogram.h"
81	#include "CLHEP/Hist/Tuple.h"
82
83	HepTupleManager* hbookManager;
84
85	main()
86	{
87
88	// Setup
89
90	G4int niter=1e4;
91	G4int imat=2;
92	G4int verboseLevel=1;
93	G4int processID=6;
94
95	G4cout << "How many interactions? [10], Which material? [3], which Verbose Level? [1]" << G4endl;
96	G4cin >> niter >> imat >> verboseLevel;
97
98	G4cout<<"which process?"<<G4endl<<std::setw(60)<<"[1] = G4LowEnergyPhotoElectric, [2] = G4LowEnergyCompton"<<G4endl;
99	G4cout<<std::setw(60)<<"[3] = G4LowEnergyRayleigh, [4] = G4LowEnergyGammaconversion"<<G4endl;
100	G4cout<<std::setw(60)<<"[5] = G4LowEnergyBremstrahlung"<<"[6] = G4LowEnergyIonisation"<<G4endl;
101
102	G4cin >> processID;
103
104	G4double InitEnergy = 1e-3, InitX = 0., InitY = 0., InitZ = 1.;
105	G4cout<<"Enter the initial particle energy E and its direction"<<G4endl;
106	G4cin >> InitEnergy >> InitX >> InitY >> InitZ;
107
108	G4double gammaCut = 1e-3, electronCut = 1e-3;
109	G4cout<<"Set photons 1e-3 mm and electrons cuts 1e-3 mm"<<G4endl;
110	G4cin>>gammaCut>>electronCut;
111
112	//-------- write results onto a file --------
113
114	// std::ofstream outFile1( "lowenergypri.out", std::ios::out);
115	// std::ofstream outFile2( "lowenergysec.out", std::ios::out);
116	// std::ofstream outFile3( "lowenergymfp.out", std::ios::out);
117
118	// outFile1.setf( std::ios::scientific, std::ios::floatfield);
119	// outFile2.setf( std::ios::scientific, std::ios::floatfield);
120	// outFile3.setf( std::ios::scientific, std::ios::floatfield);
121
122	// outFile1.setf(std::ios::left);
123	// outFile2.setf(std::ios::left);
124	// outFile3.setf(std::ios::left);
125
126	G4cout.setf( std::ios::scientific, std::ios::floatfield );
127	// -------------------------------------------------------------------
128
129	// ALE ---- HBOOK initialization
130	if(processID == 1){
131
132	hbookManager = new HBookFile("photoelectric.hbook", 58);
133	assert (hbookManager != 0);
134	}
135	else if(processID == 2){
136
137	hbookManager = new HBookFile("compton.hbook", 58);
138	assert (hbookManager != 0);
139	}
140	else if(processID == 3){
141
142	hbookManager = new HBookFile("rayleigh.hbook", 58);
143	assert (hbookManager != 0);
144	}
145	else if(processID == 4){
146
147	hbookManager = new HBookFile("gammaconv.hbook", 58);
148	assert (hbookManager != 0);
149	}
150	else if(processID == 5){
151
152	hbookManager = new HBookFile("bremstrahlung.hbook", 58);
153	assert (hbookManager != 0);
154	}
155	else if(processID == 6){
156
157	hbookManager = new HBookFile("ionisation.hbook", 58);
158	assert (hbookManager != 0);
159	}
160	// ALE ---- Book a histogram and ntuples
161	G4cout<<"Hbook file name: "<<((HBookFile*) hbookManager)->filename()<<G4endl;
162
163	// ---- primary ntuple ------
164	HepTuple* ntuple1 = hbookManager->ntuple("Primary Ntuple");
165	assert (ntuple1 != 0);
166
167	// ---- secondaries ntuple ------
168	HepTuple* ntuple2 = hbookManager->ntuple("Secondaries Ntuple");
169	assert (ntuple2 != 0);
170
171	// ---- mfp ntuple ------
172	HepTuple* ntuple3 = hbookManager->ntuple("MeanFreePath Ntuple");
173	assert (ntuple3 != 0);
174
175	// ---- secondaries histos ----
176	HepHistogram* hEKin;
177	hEKin = hbookManager->histogram("Kinetic Energy", 100,0.,200.);
178	assert (hEKin != 0);
179
180	HepHistogram* hP;
181	hP = hbookManager->histogram("Momentum", 100,0.,1000.);
182	assert (hP != 0);
183
184	HepHistogram* hNSec;
185	hNSec = hbookManager->histogram("Number of secondaries", 40,0.,40.);
186	assert (hNSec != 0);
187
188	HepHistogram* hDebug;
189	hDebug = hbookManager->histogram("Debug", 100,0.,200.);
190	assert (hDebug != 0);
191
192
193	//--------- Materials definition ---------
194
195	G4Material* Be = new G4Material("Beryllium", 4., 9.01g/mole, 1.848g/cm3);
196	G4Material* Graphite = new G4Material("Graphite",6., 12.00g/mole, 2.265g/cm3 );
197	G4Material* Al = new G4Material("Aluminium", 13., 26.98g/mole, 2.7 g/cm3);
198	G4Material* LAr = new G4Material("LArgon", 18., 39.95g/mole, 1.393g/cm3);
199	G4Material* Fe = new G4Material("Iron", 26., 55.85g/mole, 7.87g/cm3);
200	G4Material* Cu = new G4Material("Copper", 29., 63.55g/mole, 8.96g/cm3);
201	G4Material* W = new G4Material("Tungsten", 74., 183.85g/mole, 19.30g/cm3);
202	G4Material* Pb = new G4Material("Lead", 82., 207.19g/mole, 11.35g/cm3);
203	G4Material* U = new G4Material("Uranium", 92., 238.03g/mole, 18.95g/cm3);
204
205	G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
206	G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
207	G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
208	G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
209	G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole);
210
211	G4Material* maO = new G4Material("Oxygen", 8., 16.00g/mole, 1.1g/cm3);
212
213	G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2);
214	water->AddElement(H,2);
215	water->AddElement(O,1);
216
217	G4Material* ethane = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
218	ethane->AddElement(H,6);
219	ethane->AddElement(C,2);
220
221	G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2);
222	csi->AddElement(Cs,1);
223	csi->AddElement(I,1);
224
225	static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
226
227	G4cout<<"The material is: "<<(*theMaterialTable)(imat)->GetName()<<G4endl;
228
229	// Geometry definitions
230	G4Box* theFrame = new G4Box ("Frame",92mm, 92mm, 92*mm);
231
232	G4LogicalVolume* LogicalFrame = new G4LogicalVolume(theFrame,
233	(*theMaterialTable)(imat),
234	"LFrame", 0, 0, 0);
235
236	G4PVPlacement* PhysicalFrame = new G4PVPlacement(0,G4ThreeVector(),
237	"PFrame",LogicalFrame,0,false,0);
238
239	// the center-of-mass of the cube should be located at the origin!
240
241	//--------- Particle definition ---------
242
243	G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
244	G4ParticleDefinition* positron = G4Positron::PositronDefinition();
245	G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
246	G4ParticleDefinition* proton = G4Proton::ProtonDefinition();
247
248	//--------- Processes definition ---------
249
250	G4ProcessManager* theGammaProcessManager = new G4ProcessManager(gamma);
251	gamma->SetProcessManager(theGammaProcessManager);
252
253	G4ProcessManager* theElectronProcessManager = new G4ProcessManager(electron);
254	electron->SetProcessManager(theElectronProcessManager);
255
256	G4ProcessManager* thePositronProcessManager = new G4ProcessManager(positron);
257	positron->SetProcessManager(thePositronProcessManager);
258
259	G4ProcessManager* theProtonProcessManager = new G4ProcessManager(proton);
260	proton->SetProcessManager(theProtonProcessManager);
261
262	G4LowEnergyPhotoElectric PhotoElectricProcess;
263	G4LowEnergyCompton ComptonProcess;
264	G4LowEnergyRayleigh RayleighProcess;
265	G4LowEnergyGammaConversion GammaConvProcess;
266
267	theGammaProcessManager->AddDiscreteProcess(&PhotoElectricProcess);
268	theGammaProcessManager->AddDiscreteProcess(&ComptonProcess);
269	theGammaProcessManager->AddDiscreteProcess(&RayleighProcess);
270	theGammaProcessManager->AddDiscreteProcess(&GammaConvProcess);
271
272	G4LowEnergyIonisation IonisationProcess;
273	theElectronProcessManager->AddProcess(&IonisationProcess);
274	theElectronProcessManager->SetProcessOrdering(&IonisationProcess,idxAlongStep,1);
275	theElectronProcessManager->SetProcessOrdering(&IonisationProcess,idxPostStep,1);
276
277	G4LowEnergyBremsstrahlung BremstrahlungProcess;
278	theElectronProcessManager->AddProcess(&BremstrahlungProcess);
279	theElectronProcessManager->SetProcessOrdering(&BremstrahlungProcess,idxAlongStep,1);
280	theElectronProcessManager->SetProcessOrdering(&BremstrahlungProcess,idxPostStep,1);
281
282	G4eIonisation IonisationPlusProcess;
283	thePositronProcessManager->AddProcess(&IonisationPlusProcess);
284	thePositronProcessManager->SetProcessOrdering(&IonisationPlusProcess,idxAlongStep,1);
285	thePositronProcessManager->SetProcessOrdering(&IonisationPlusProcess,idxPostStep,1);
286
287	G4hIonisation hIonisationProcess;
288	theProtonProcessManager->AddProcess(&hIonisationProcess);
289	theProtonProcessManager->SetProcessOrdering(&hIonisationProcess,idxAlongStep,1);
290	theProtonProcessManager->SetProcessOrdering(&hIonisationProcess,idxPostStep,1);
291
292	G4ForceCondition* condition;
293
294	// ------- set cut and Build CrossSection Tables -------
295	//
296	G4Gamma::SetEnergyRange(2.5e-4MeV,1e5MeV);
297	G4Electron::SetEnergyRange(2.5e-4MeV,1e5MeV);
298	G4Positron::SetEnergyRange(2.5e-4MeV,1e5MeV);
299
300	gamma->SetCuts(1e-6*mm);
301	electron->SetCuts(1e-6*mm);
302	positron->SetCuts(1e-6*mm);
303
304	G4cout<<"the cut in energy for gamma in: "<<(*theMaterialTable)(imat)->GetName()
305	<<" is: "<<G4Gamma::GetCutsInEnergy()[imat]<<G4endl;
306	G4cout<<"the cut in energy for e- in: "<<(*theMaterialTable)(imat)->GetName()
307	<<" is: "<<G4Electron::GetCutsInEnergy()[imat]<<G4endl;
308
309	// -------- create 1 Dynamic Particle ----
310
311	G4double gammaEnergy = InitEnergy*MeV;
312
313	G4ParticleMomentum gammaDirection(InitX,InitY,InitZ);
314
315
316	G4DynamicParticle photon(G4Gamma::Gamma(),gammaDirection,gammaEnergy);
317	G4DynamicParticle elecT(G4Electron::Electron(),gammaDirection,gammaEnergy);
318
319	//--------- track definition (for this test ONLY!)------------
320
321	G4ThreeVector aPosition(0.,0.,0.);
322	G4double aTime = 0. ;
323
324	G4Track* ptrack;
325	if(processID != 5 && processID != 6)
326	ptrack = new G4Track(&photon,aTime,aPosition) ;
327	else
328	ptrack = new G4Track(&elecT,aTime,aPosition) ;
329
330	G4Track& aTrack = (*ptrack) ;
331
332	// do I really need this?
333
334	G4GRSVolume* touche = new G4GRSVolume(PhysicalFrame, NULL, aPosition);
335	ptrack->SetTouchable(touche);
336
337	// -------- create 1 Step (for this test only)----
338
339	G4Step* Step = new G4Step();
340	G4Step& aStep = (*Step);
341	Step->SetTrack(ptrack);
342
343	// --------- check applicability
344	G4ParticleDefinition* PhotonDefinition = photon.GetDefinition();
345	G4ParticleDefinition* ElectronDefinition = elecT.GetDefinition();
346
347	if(!PhotoElectricProcess.IsApplicable(*PhotonDefinition) \|\|
348	!ComptonProcess.IsApplicable(*PhotonDefinition) \|\|
349	!RayleighProcess.IsApplicable(*PhotonDefinition) \|\|
350	!GammaConvProcess.IsApplicable(*PhotonDefinition) \|\|
351	!BremstrahlungProcess.IsApplicable(*ElectronDefinition)\|\|
352	!IonisationProcess.IsApplicable(*ElectronDefinition)) {
353
354	G4cout
355	<< PhotonDefinition->GetParticleName()
356	<< " is not a Photon! or "
357	<< ElectronDefinition->GetParticleName()
358	<<" is not an Electron"<<G4endl;
359	G4Exception("FAIL: * exit program *\n");
360	// return ;
361	}
362
363	// PhotoElectricProcess.SetVerboseLevel(verboseLevel);
364
365	// Initialize the physics tables for ALL processes
366	// IonisationProcess.MinusNbOfProcesses();
367	IonisationProcess.BuildPhysicsTable(*electron);
368	BremstrahlungProcess.BuildPhysicsTable(*electron);
369
370	IonisationPlusProcess.MinusNbOfProcesses();
371	IonisationPlusProcess.BuildPhysicsTable(*positron);
372
373	hIonisationProcess.BuildPhysicsTable(*proton);
374
375	if(processID == 1){
376
377	PhotoElectricProcess.BuildPhysicsTable(*gamma);
378	}
379
380	else if(processID == 2) {
381
382	ComptonProcess.BuildPhysicsTable(*gamma);
383	}
384
385	else if(processID == 3) {
386
387	RayleighProcess.BuildPhysicsTable(*gamma);
388	}
389
390	else if(processID == 4) {
391
392	GammaConvProcess.BuildPhysicsTable(*gamma);
393	}
394	else if(processID == 5) {
395
396	cout<<"**** BR BuildPhysicsTable:Table already Built *****"<<G4endl;
397	}
398	else if(processID == 6) {
399
400	cout<<"**** IO BuildPhysicsTable:Table already Built *****"<<G4endl;
401	}
402
403	else {
404
405	G4Exception("No such process!\n");
406	}
407
408	// Mean FreePath Test
409	G4Material* apttoMaterial ;
410	G4String MaterialName ;
411	///***********************************************************************
412	// UNCOMMENT LATER OR NEVER I HAVE IT IN 6 diff files
413	//***********************************************************************
414	G4double minArg = 100eV,maxArg = 100GeV, argStp;
415	const G4int pntNum = 300;
416	G4double Tkin[pntNum+1];
417	G4double meanFreePath ;
418
419	argStp = (std::log10(maxArg)-std::log10(minArg))/pntNum;
420
421	for(G4int d = 0; d < pntNum+1; d++){
422
423	Tkin[d] = std::pow(10,(std::log10(minArg) + d*argStp));
424	}
425
426	for ( G4int J = 0 ; J < G4Material::GetNumberOfMaterials() ; J++ ){
427
428	apttoMaterial = (*theMaterialTable)[ J ] ;
429	MaterialName = apttoMaterial->GetName() ;
430
431	LogicalFrame->SetMaterial(apttoMaterial);
432
433	for (G4int i=0 ; i<pntNum; i++){
434
435	if(processID != 5 && processID != 6)
436	photon.SetKineticEnergy(Tkin[i]);
437	else
438	elecT.SetKineticEnergy(Tkin[i]);
439
440	if(processID == 1){
441
442	meanFreePath = PhotoElectricProcess.GetMeanFreePath(aTrack, 1., condition);
443	}
444
445	else if(processID == 2) {
446
447	meanFreePath = ComptonProcess.GetMeanFreePath(aTrack, 1., condition);
448	}
449
450	else if(processID == 3) {
451
452	meanFreePath = RayleighProcess.GetMeanFreePath(aTrack, 1., condition);
453	}
454
455	else if(processID == 4) {
456
457	meanFreePath = GammaConvProcess.GetMeanFreePath(aTrack, 1., condition);
458	}
459
460	else if(processID == 5) {
461
462	meanFreePath = BremstrahlungProcess.GetMeanFreePath(aTrack, 1., condition);
463	}
464
465	else if(processID == 6) {
466
467	meanFreePath = IonisationProcess.GetMeanFreePath(aTrack, 1., condition);
468	}
469
470	else {
471
472	G4Exception("No such process!\n");
473	}
474
475	ntuple3->column("matind",J);
476	ntuple3->column("kinen",Tkin[i]);
477	ntuple3->column("mfp",meanFreePath);
478	ntuple3->dumpData();
479
480	// outFile3<<std::setw(4)<<J<<std::setw(14)<<Tkin[i]<<std::setw(14)<<meanFreePath<<G4endl;
481	}
482	}// for loop on materials
483	//END OF COMMENT */
484	// --------- Test the DoIt for the Photon Absorption
485
486	apttoMaterial = (*theMaterialTable)(imat) ;
487
488	LogicalFrame->SetMaterial(apttoMaterial);
489	if(processID != 5 && processID != 6){
490	photon.SetKineticEnergy(InitEnergy*MeV);
491	photon.SetMomentumDirection(InitX, InitY, InitZ);
492	}
493	else{
494	elecT.SetKineticEnergy(InitEnergy*MeV);
495	elecT.SetMomentumDirection(InitX, InitY, InitZ);
496	}
497	// PostStepDoIt calls
498	G4int iteration = 0;
499
500	G4VParticleChange* adummy;
501	G4Track* aFinalParticle;
502	G4String aParticleName;
503
504	do {
505
506	if(processID == 1){
507
508	adummy = PhotoElectricProcess.PostStepDoIt(aTrack, aStep);
509
510	}
511
512	else if(processID == 2) {
513
514	adummy = ComptonProcess.PostStepDoIt(aTrack, aStep);
515	}
516
517	else if(processID == 3) {
518
519	adummy = RayleighProcess.PostStepDoIt(aTrack, aStep);
520	}
521
522	else if(processID == 4) {
523
524	adummy = GammaConvProcess.PostStepDoIt(aTrack, aStep);
525	}
526
527	else if(processID == 5) {
528
529	adummy = BremstrahlungProcess.PostStepDoIt(aTrack, aStep);
530	}
531
532	else if(processID == 6) {
533
534	adummy = IonisationProcess.PostStepDoIt(aTrack, aStep);
535	}
536
537	else {
538
539	G4Exception("No such process!\n");
540	}
541
542	G4ParticleChange* aParticleChange = (G4ParticleChange*) adummy;
543
544	// ------------ book primary physical quantities -------------
545	G4double pEnChange = 0, pxChange = 0, pyChange = 0, pzChange = 0, PChange = 0;
546
547	pEnChange = aParticleChange->GetEnergyChange();
548	pxChange = aParticleChange->GetMomentumChange()->x();
549	pyChange = aParticleChange->GetMomentumChange()->y();
550	pzChange = aParticleChange->GetMomentumChange()->z();
551	PChange = std::sqrt(pxChangepxChange+pyChangepyChange+pzChange*pzChange);
552
553	// ---- secondaries histos ----
554	G4cout<<"E and p of the primary particle: "<<pEnChange<<" "<<pxChange<<" "
555	<<pyChange<<" "<<pzChange<<G4endl;
556
557	ntuple1->column("ench", pEnChange);
558	ntuple1->column("pxch", pxChange);
559	ntuple1->column("pych", pyChange);
560	ntuple1->column("pzch", pzChange);
561	ntuple1->column("pch", PChange);
562
563	// outFile1<<std::setw(13)<<pEnChange<<std::setw(13)<<pxChange<<std::setw(13)
564	// <<pyChange<<std::setw(13)<<pzChange<<G4endl;
565
566	ntuple1->dumpData();
567
568	// ------------ book secondaries physical quantities ---------
569
570	G4double e = 0, eKin = 0, Px = 0, Py = 0, Pz = 0, P = 0, ShID = 0;
571
572	// hNSec->accumulate(aParticleChange->GetNumberOfSecondaries());
573	// hDebug->accumulate(aParticleChange->GetLocalEnergyDeposit());
574
575	for (G4int i = 0; i < (aParticleChange->GetNumberOfSecondaries()); i++) {
576
577	// The following two items should be filled per event, not
578	// per secondary; filled here just for convenience, to avoid
579	// complicated logic to dump ntuple when there are no secondaries
580
581	ntuple2->column("nsec",aParticleChange->GetNumberOfSecondaries());
582	ntuple2->column("deposit",aParticleChange->GetLocalEnergyDeposit());
583
584	aFinalParticle = aParticleChange->GetSecondary(i) ;
585
586	e = aFinalParticle->GetTotalEnergy();
587	eKin = aFinalParticle->GetKineticEnergy();
588	Px = (aFinalParticle->GetMomentum()).x();
589	Py = (aFinalParticle->GetMomentum()).y();
590	Pz = (aFinalParticle->GetMomentum()).z();
591	P = std::sqrt(PxPx+PyPy+Pz*Pz);
592	if(processID == 1){
593
594	ShID = PhotoElectricProcess.GetTransitionShell(i);
595	}
596	else if(processID == 6){
597
598	ShID = IonisationProcess.GetTransitionShell(i);
599	}
600
601	aParticleName = aFinalParticle->GetDefinition()->GetParticleName();
602	G4cout<<aParticleName<<": "
603	<<" "<<e<<" "<<eKin<<" "<<Px<<" "<<Py<<" "<<Pz<<" ***"<<G4endl;
604	hEKin->accumulate(eKin);
605	hP->accumulate(std::sqrt(PxPx+PyPy+Pz*Pz));
606
607	G4int ptype;
608	if(aParticleName == "gamma") ptype = 0;
609	else if(aParticleName == "e-") ptype = -1;
610	else if(aParticleName == "e+") ptype = 1;
611
612	// Fill the secondaries ntuple
613	ntuple2->column("px", Px);
614	ntuple2->column("py", Py);
615	ntuple2->column("pz", Pz);
616	ntuple2->column("p", P);
617	ntuple2->column("e", e);
618	ntuple2->column("ekin", eKin);
619	ntuple2->column("ptype", ptype);
620	ntuple2->column("sh", ShID);
621
622	ntuple2->dumpData();
623
624	// Print secondaries on a file
625	// outFile2<<std::setw(3)<<aParticleChange->GetNumberOfSecondaries()<<std::setw(14)
626	// <<aParticleChange->GetLocalEnergyDeposit()<<std::setw(14)<<e<<std::setw(14)
627	// <<eKin<<std::setw(14)<<std::sqrt(PxPx+PyPy+Pz*Pz)<<std::setw(14)<<Px<<std::setw(14)
628	// <<Py<<std::setw(14)<<Pz<<std::setw(3)<<ptype<<G4endl;
629
630	delete aParticleChange->GetSecondary(i);
631	}
632
633	aParticleChange->Clear();
634	iteration++;
635	G4cout << "******* Iteration = " << iteration << G4endl;
636
637	} while (iteration < niter) ;
638	cout<<"Iteration number: "<<G4endl;
639	hbookManager->write();
640	delete hbookManager;
641
642	// delete materials and elements
643	delete Be;
644	delete Graphite;
645	delete Al;
646	delete LAr;
647	delete Fe;
648	delete Cu;
649	delete W;
650	delete Pb;
651	delete U;
652	delete H;
653	delete maO;
654	delete C;
655	delete Cs;
656	delete I;
657	delete O;
658	delete water;
659	delete ethane;
660	delete csi;
661	delete Step;
662	delete touche;
663	// outFile1.close();
664	// outFile2.close();
665	// outFile3.close();
666
667	cout<<"END OF THE MAIN PROGRAM"<<G4endl;
668	}
669
670
671
672
673
674
675
676
677
678
679
680
681

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

Context Navigation

source: trunk/source/processes/electromagnetic/lowenergy/test/G4LowEnergyTest.cc @ 1350

Download in other formats: