Context Navigation

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

Last change on this file since 1199 was 1199, checked in by garnier, 15 years ago
nvx fichiers dans CVS
File size: 15.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: G4PhotoelectricTest.cc,v 1.11 2006/06/29 19:44:30 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: G4PhotoelectricTest
36	//
37	// Author: Maria Grazia Pia
38	//
39	// Creation date: 1 May 2001
40	//
41	// Modifications:
42	//
43	// -------------------------------------------------------------------
44
45	#include "globals.hh"
46	#include "G4ios.hh"
47	#include <fstream>
48	#include <iomanip>
49
50	#include "G4Material.hh"
51	#include "G4VContinuousDiscreteProcess.hh"
52	#include "G4ProcessManager.hh"
53	#include "G4LowEnergyBremsstrahlung.hh"
54	#include "G4eBremsstrahlung.hh"
55	#include "G4LowEnergyPhotoElectric.hh"
56	#include "G4PhotoElectricEffect.hh"
57	#include "G4LowEnergyIonisation.hh"
58	#include "G4eIonisation.hh"
59	#include "G4VeLowEnergyLoss.hh"
60	#include "G4EnergyLossTables.hh"
61	#include "G4VParticleChange.hh"
62	#include "G4ParticleChange.hh"
63	#include "G4DynamicParticle.hh"
64	#include "G4Electron.hh"
65	#include "G4Positron.hh"
66	#include "G4Gamma.hh"
67
68	#include "G4Box.hh"
69	#include "G4PVPlacement.hh"
70
71	#include "G4Step.hh"
72	#include "G4GRSVolume.hh"
73
74	#include "G4UnitsTable.hh"
75	#include "CLHEP/Hist/TupleManager.h"
76	#include "CLHEP/Hist/HBookFile.h"
77	#include "CLHEP/Hist/Histogram.h"
78	#include "CLHEP/Hist/Tuple.h"
79
80	HepTupleManager* hbookManager;
81
82	int main()
83	{
84
85	// Setup
86
87	G4int nIterations = 100000;
88	G4int materialId = 3;
89
90	G4cout.setf( ios::scientific, ios::floatfield );
91
92	// -------------------------------------------------------------------
93
94	// ---- HBOOK initialization
95
96
97	hbookManager = new HBookFile("photoeltest.hbook", 58);
98	assert (hbookManager != 0);
99
100	// ---- Book a histogram and ntuples
101	G4cout<<"Hbook file name: "<<((HBookFile*) hbookManager)->filename()<<endl;
102
103	// ---- primary ntuple ------
104	HepTuple* ntuple1 = hbookManager->ntuple("Primary Ntuple");
105	assert (ntuple1 != 0);
106
107	// ---- secondary ntuple ------
108	HepTuple* ntuple2 = hbookManager->ntuple("Secondary Ntuple");
109	assert (ntuple2 != 0);
110
111	// ---- secondaries histos ----
112	HepHistogram* hEKin;
113	hEKin = hbookManager->histogram("Kinetic Energy", 100,0.,200.);
114	assert (hEKin != 0);
115
116	HepHistogram* hP;
117	hP = hbookManager->histogram("Momentum", 100,0.,1000.);
118	assert (hP != 0);
119
120	HepHistogram* hNSec;
121	hNSec = hbookManager->histogram("Number of secondaries", 40,0.,40.);
122	assert (hNSec != 0);
123
124	HepHistogram* hDebug;
125	hDebug = hbookManager->histogram("Debug", 100,0.,200.);
126	assert (hDebug != 0);
127
128
129	//--------- Materials definition ---------
130
131	G4Material* Be = new G4Material("Beryllium", 4., 9.01g/mole, 1.848g/cm3);
132	G4Material* Graphite = new G4Material("Graphite",6., 12.00g/mole, 2.265g/cm3 );
133	G4Material* Al = new G4Material("Aluminium", 13., 26.98g/mole, 2.7 g/cm3);
134	G4Material* Si = new G4Material("Silicon", 14., 28.055g/mole, 2.33g/cm3);
135	G4Material* LAr = new G4Material("LArgon", 18., 39.95g/mole, 1.393g/cm3);
136	G4Material* Fe = new G4Material("Iron", 26., 55.85g/mole, 7.87g/cm3);
137	G4Material* Cu = new G4Material("Copper", 29., 63.55g/mole, 8.96g/cm3);
138	G4Material* W = new G4Material("Tungsten", 74., 183.85g/mole, 19.30g/cm3);
139	G4Material* Pb = new G4Material("Lead", 82., 207.19g/mole, 11.35g/cm3);
140	G4Material* U = new G4Material("Uranium", 92., 238.03g/mole, 18.95g/cm3);
141
142	G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
143	G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
144	G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
145	G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
146	G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole);
147
148	G4Material* maO = new G4Material("Oxygen", 8., 16.00g/mole, 1.1g/cm3);
149
150	G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2);
151	water->AddElement(H,2);
152	water->AddElement(O,1);
153
154	G4Material* ethane = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
155	ethane->AddElement(H,6);
156	ethane->AddElement(C,2);
157
158	G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2);
159	csi->AddElement(Cs,1);
160	csi->AddElement(I,1);
161
162
163	// Interactive set-up
164
165	G4cout << "How many interactions? " << G4endl;
166	G4cin >> nIterations;
167
168	if (nIterations <= 0) G4Exception("Wrong input");
169
170	G4double initEnergy = 1*MeV;
171	G4double initX = 0.;
172	G4double initY = 0.;
173	G4double initZ = 1.;
174
175	G4cout << "Enter the initial particle energy E (MeV)" << G4endl;
176	G4cin >> initEnergy ;
177
178	initEnergy = initEnergy * MeV;
179
180	if (initEnergy <= 0.) G4Exception("Wrong input");
181
182	static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
183
184	G4int nMaterials = G4Material::GetNumberOfMaterials();
185
186	G4cout << "Available materials are: " << G4endl;
187	for (G4int mat = 0; mat < nMaterials; mat++)
188	{
189	G4cout << mat << ") "
190	<< (*theMaterialTable)[mat]->GetName()
191	<< G4endl;
192	}
193
194	G4cout << "Which material? " << G4endl;
195	G4cin >> materialId;
196
197	G4Material* material = (*theMaterialTable)[materialId] ;
198
199	G4cout << "The selected material is: "
200	<< material->GetName()
201	<< G4endl;
202
203	G4double dimX = 1*mm;
204	G4double dimY = 1*mm;
205	G4double dimZ = 1*mm;
206
207	// Geometry
208
209	G4Box* theFrame = new G4Box ("Frame",dimX, dimY, dimZ);
210
211	G4LogicalVolume* logicalFrame = new G4LogicalVolume(theFrame,
212	(*theMaterialTable)[materialId],
213	"LFrame", 0, 0, 0);
214	logicalFrame->SetMaterial(material);
215
216	G4PVPlacement* physicalFrame = new G4PVPlacement(0,G4ThreeVector(),
217	"PFrame",logicalFrame,0,false,0);
218
219	// Particle definitions
220
221	G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
222	G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
223	G4ParticleDefinition* positron = G4Positron::PositronDefinition();
224
225	gamma->SetCuts(1e-3*mm);
226	electron->SetCuts(1e-3*mm);
227	positron->SetCuts(1e-3*mm);
228
229	// Processes
230
231	G4int processType;
232	G4cout << "LowEnergy [1] or Standard [2] Bremsstrahlung?" << G4endl;
233	cin >> processType;
234	if ( !(processType == 1 \|\| processType == 2))
235	{
236	G4Exception("Wrong input");
237	}
238
239	G4VContinuousDiscreteProcess* bremProcess;
240	G4VContinuousDiscreteProcess* ioniProcess;
241
242	if (processType == 1)
243	{
244	bremProcess = new G4LowEnergyBremsstrahlung;
245	ioniProcess = new G4LowEnergyIonisation;
246	}
247	else
248	{
249	bremProcess = new G4eBremsstrahlung;
250	ioniProcess = new G4eIonisation;
251	}
252
253	G4ProcessManager* eProcessManager = new G4ProcessManager(electron);
254	electron->SetProcessManager(eProcessManager);
255	eProcessManager->AddProcess(bremProcess);
256
257	G4ProcessManager* positronProcessManager = new G4ProcessManager(positron);
258	positron->SetProcessManager(positronProcessManager);
259	positronProcessManager->AddProcess(bremProcess);
260
261	// Initialize the physics tables
262	bremProcess->BuildPhysicsTable(*electron);
263	ioniProcess->BuildPhysicsTable(*electron);
264
265	// Photon process
266	G4VDiscreteProcess* photonProcess;
267	if (processType == 1)
268	{
269	photonProcess = new G4LowEnergyPhotoElectric;
270	}
271	else
272	{
273	photonProcess = new G4PhotoElectricEffect;
274	}
275
276	G4ProcessManager* gProcessManager = new G4ProcessManager(gamma);
277	gamma->SetProcessManager(gProcessManager);
278	gProcessManager->AddProcess(photonProcess);
279	photonProcess->BuildPhysicsTable(*gamma);
280
281	// Create a DynamicParticle
282
283	G4double gEnergy = initEnergy*MeV;
284	G4ParticleMomentum gDirection(initX,initY,initZ);
285	G4DynamicParticle dynamicPhoton(G4Gamma::Gamma(),gDirection,gEnergy);
286
287	// Track
288
289	G4ThreeVector aPosition(0.,0.,0.);
290	G4ThreeVector newPosition(0.,0.,1.*mm);
291	G4double aTime = 0. ;
292
293	G4Track* gTrack = new G4Track(&dynamicPhoton,aTime,aPosition);
294
295	// do I really need this?
296
297	G4GRSVolume* touche = new G4GRSVolume(physicalFrame, 0, aPosition);
298	gTrack->SetTouchable(touche);
299
300	// Step
301
302	G4Step* step = new G4Step();
303	step->SetTrack(gTrack);
304
305	G4StepPoint* aPoint = new G4StepPoint();
306	aPoint->SetPosition(aPosition);
307	aPoint->SetMaterial(material);
308	G4double safety = 10000.*cm;
309	aPoint->SetSafety(safety);
310	step->SetPreStepPoint(aPoint);
311	G4StepPoint* newPoint = new G4StepPoint();
312	newPoint->SetPosition(newPosition);
313	newPoint->SetMaterial(material);
314	step->SetPostStepPoint(newPoint);
315
316	// Check applicability
317
318	if (! (photonProcess->IsApplicable(*gamma)))
319	{
320	G4Exception("Not Applicable");
321	}
322
323	// --------- Test the DoIt
324
325	G4cout << "DoIt in material " << material->GetName() << G4endl;
326
327	for (G4int iter=0; iter<nIterations; iter++)
328	{
329	step->SetStepLength(1*micrometer);
330
331	gTrack->SetStep(step);
332
333	G4cout << "Iteration = " << iter
334	<< " - Step Length = "
335	<< step->GetStepLength()/mm << " mm "
336	<< G4endl;
337
338	G4VParticleChange* dummy;
339	dummy = photonProcess->PostStepDoIt(gTrack, step);
340
341	G4ParticleChange* particleChange = (G4ParticleChange*) dummy;
342
343	// Primary physical quantities
344
345	G4double energyChange = particleChange->GetEnergyChange();
346	G4double dedx = initEnergy - energyChange ;
347	G4double dedxNow = dedx / (step->GetStepLength());
348
349	G4ThreeVector eChange = particleChange->CalcMomentum(energyChange,
350	(*particleChange->GetMomentumChange()),
351	particleChange->GetMassChange());
352
353	G4double pxChange = eChange.x();
354	G4double pyChange = eChange.y();
355	G4double pzChange = eChange.z();
356	G4double pChange = std::sqrt(pxChangepxChange + pyChangepyChange + pzChange*pzChange);
357
358	G4double xChange = particleChange->GetPositionChange()->x();
359	G4double yChange = particleChange->GetPositionChange()->y();
360	G4double zChange = particleChange->GetPositionChange()->z();
361	G4double thetaChange = particleChange->GetPositionChange()->theta();
362
363	G4cout << "---- Primary after the step ---- " << G4endl;
364
365	// G4cout << "Position (x,y,z) = "
366	// << xChange << " "
367	// << yChange << " "
368	// << zChange << " "
369	// << G4endl;
370
371	G4cout << "---- Energy: " << energyChange/MeV << " MeV, "
372	<< "(px,py,pz): ("
373	<< pxChange/MeV << ","
374	<< pyChange/MeV << ","
375	<< pzChange/MeV << ") MeV"
376	<< G4endl;
377
378	G4cout << "---- Energy loss (dE) = " << dedx/keV << " keV" << G4endl;
379	// G4cout << "Stopping power (dE/dx)=" << dedxNow << G4endl;
380
381	// Secondaries
382
383	ntuple1->column("eprimary", initEnergy);
384	ntuple1->column("energyf", energyChange);
385	ntuple1->column("de", dedx);
386	ntuple1->column("dedx", dedxNow);
387	ntuple1->column("pxch", xChange);
388	ntuple1->column("pych", pyChange);
389	ntuple1->column("pzch", pzChange);
390	ntuple1->column("pch", zChange);
391	ntuple1->column("thetach", thetaChange);
392
393	// Secondaries physical quantities
394
395	hNSec->accumulate(particleChange->GetNumberOfSecondaries());
396	hDebug->accumulate(particleChange->GetLocalEnergyDeposit());
397
398	G4int nElectrons = 0;
399	G4int nPositrons = 0;
400	G4int nPhotons = 0;
401
402	for (G4int i = 0; i < (particleChange->GetNumberOfSecondaries()); i++)
403	{
404	// The following two items should be filled per event, not
405	// per secondary; filled here just for convenience, to avoid
406	// complicated logic to dump ntuple when there are no secondaries
407
408	G4Track* finalParticle = particleChange->GetSecondary(i) ;
409
410	G4double e = finalParticle->GetTotalEnergy();
411	G4double eKin = finalParticle->GetKineticEnergy();
412	G4double px = (finalParticle->GetMomentum()).x();
413	G4double py = (finalParticle->GetMomentum()).y();
414	G4double pz = (finalParticle->GetMomentum()).z();
415	G4double theta = (finalParticle->GetMomentum()).theta();
416	G4double p = std::sqrt(pxpx+pypy+pz*pz);
417
418	if (e > initEnergy)
419	{
420	G4cout << "WARNING: eFinal > eInit " << G4endl;
421	// << e
422	// << " > " initEnergy
423
424	}
425
426	G4String particleName = finalParticle->GetDefinition()->GetParticleName();
427	G4cout << "==== Final "
428	<< particleName << " "
429	<< "energy: " << e/MeV << " MeV, "
430	<< "eKin: " << eKin/MeV << " MeV, "
431	<< "(px,py,pz): ("
432	<< px/MeV << ","
433	<< py/MeV << ","
434	<< pz/MeV << ") MeV "
435	<< G4endl;
436
437	hEKin->accumulate(eKin);
438	hP->accumulate(p);
439
440	G4int partType = 0;
441	if (particleName == "e-")
442	{
443	partType = 1;
444	nElectrons++;
445	}
446	else if (particleName == "e+")
447	{
448	partType = 2;
449	nPositrons++;
450	}
451	else if (particleName == "gamma")
452	{
453	partType = 3;
454	nPhotons++;
455	}
456	// Fill the secondaries ntuple
457	ntuple2->column("eprimary",initEnergy);
458	ntuple2->column("px", px);
459	ntuple2->column("py", py);
460	ntuple2->column("pz", pz);
461	ntuple2->column("p", p);
462	ntuple2->column("e", e);
463	ntuple2->column("theta", theta);
464	ntuple2->column("ekin", eKin);
465	ntuple2->column("type", partType);
466
467	ntuple2->dumpData();
468
469	delete particleChange->GetSecondary(i);
470	}
471
472	ntuple1->column("nelectrons",nElectrons);
473	ntuple1->column("npositrons",nPositrons);
474	ntuple1->column("nphotons",nPhotons);
475	ntuple1->dumpData();
476
477	particleChange->Clear();
478
479	}
480
481
482	cout << "Iteration number: " << G4endl;
483	hbookManager->write();
484	delete hbookManager;
485
486	// delete materials and elements
487	// delete Be;
488	// delete Graphite;
489	// delete Al;
490	// delete Si;
491	// delete LAr;
492	// delete Fe;
493	// delete Cu;
494	// delete W;
495	// delete Pb;
496	// delete U;
497	// delete H;
498	// delete maO;
499	// delete C;
500	// delete Cs;
501	// delete I;
502	// delete O;
503	// delete water;
504	// delete ethane;
505	// delete csi;
506	// delete step;
507	// delete touche;
508	// delete Be;
509	// delete Graphite;
510	// delete Al;
511	// delete Si;
512	// delete LAr;
513	// delete Fe;
514	// delete Cu;
515	// delete W;
516	// delete Pb;
517	// delete U;
518	// delete H;
519	// delete maO;
520	// delete C;
521	// delete Cs;
522	// delete I;
523	// delete O;
524	// delete water;
525	// delete ethane;
526	// delete csi;
527	delete step;
528	// delete touche;
529
530	cout << "END OF THE MAIN PROGRAM" << G4endl;
531	}
532
533
534
535
536
537
538
539
540
541
542
543

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

Download in other formats: