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

Last change on this file since 1350 was 1350, checked in by garnier, 15 years ago

update to last version 4.9.4
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-04-ref-00 $
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
80HepTupleManager* hbookManager;
81
82int 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.01*g/mole, 1.848*g/cm3);
132 G4Material* Graphite = new G4Material("Graphite",6., 12.00*g/mole, 2.265*g/cm3 );
133 G4Material* Al = new G4Material("Aluminium", 13., 26.98*g/mole, 2.7 *g/cm3);
134 G4Material* Si = new G4Material("Silicon", 14., 28.055*g/mole, 2.33*g/cm3);
135 G4Material* LAr = new G4Material("LArgon", 18., 39.95*g/mole, 1.393*g/cm3);
136 G4Material* Fe = new G4Material("Iron", 26., 55.85*g/mole, 7.87*g/cm3);
137 G4Material* Cu = new G4Material("Copper", 29., 63.55*g/mole, 8.96*g/cm3);
138 G4Material* W = new G4Material("Tungsten", 74., 183.85*g/mole, 19.30*g/cm3);
139 G4Material* Pb = new G4Material("Lead", 82., 207.19*g/mole, 11.35*g/cm3);
140 G4Material* U = new G4Material("Uranium", 92., 238.03*g/mole, 18.95*g/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.00*g/mole, 1.1*g/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(pxChange*pxChange + pyChange*pyChange + 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(px*px+py*py+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.