Context Navigation

G4eIonisationTest.cc

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

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

Context Navigation

source: trunk/source/processes/electromagnetic/lowenergy/test/G4eIonisationTest.cc

Download in other formats: