source: trunk/source/processes/electromagnetic/lowenergy/test/G4LowEnergyGammaConversionTest.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: 19.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: G4LowEnergyGammaConversionTest.cc,v 1.8 2006/06/29 19:44:03 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: G4LowEnergyGammaConversionTest.cc
36//
37// Author: Francesco Longo
38//
39// Creation date: 04 january 2001
40//
41// Modifications:
42//
43// -------------------------------------------------------------------
44
45#include "globals.hh"
46#include "G4ios.hh"
47#include <fstream>
48#include <iomanip>
49
50#include "G4ParticleDefinition.hh"
51#include "G4ParticleTypes.hh"
52#include "G4ParticleTable.hh"
53#include "G4Material.hh"
54#include "G4MaterialTable.hh"
55#include "G4VDiscreteProcess.hh"
56#include "G4VProcess.hh"
57#include "G4ProcessManager.hh"
58
59#include "G4LowEnergyGammaConversion.hh"
60#include "G4GammaConversion.hh"
61
62#include "G4EnergyLossTables.hh"
63#include "G4VParticleChange.hh"
64#include "G4ParticleChange.hh"
65#include "G4DynamicParticle.hh"
66
67#include "G4LowEnergyBremsstrahlung.hh"
68#include "G4LowEnergyIonisation.hh"
69#include "G4eIonisation.hh"
70#include "G4MultipleScattering.hh"
71#include "G4eIonisation.hh"
72#include "G4eBremsstrahlung.hh"
73#include "G4eplusAnnihilation.hh"
74
75//#include "G4ComptonScattering.hh"
76//#include "G4PhotoElectricEffect.hh"
77
78#include "G4Electron.hh"
79#include "G4Positron.hh"
80#include "G4Gamma.hh"
81
82#include "G4GRSVolume.hh"
83#include "G4Box.hh"
84#include "G4PVPlacement.hh"
85#include "G4Step.hh"
86
87#include "G4UnitsTable.hh"
88
89#include "CLHEP/Hist/TupleManager.h"
90#include "CLHEP/Hist/HBookFile.h"
91#include "CLHEP/Hist/Histogram.h"
92#include "CLHEP/Hist/Tuple.h"
93
94
95HepTupleManager* hbookManager;
96
97G4int main()
98{
99
100 // Setup
101
102 G4int nIterations = 100000;
103 G4int materialId = 3;
104
105 G4cout.setf( ios::scientific, ios::floatfield );
106
107 // -------------------------------------------------------------------
108
109 // ---- HBOOK initialization
110
111
112 hbookManager = new HBookFile("gammatest.hbook", 58);
113 assert (hbookManager != 0);
114
115 // ---- Book a histogram and ntuples
116 G4cout<<"Hbook file name: "<<((HBookFile*) hbookManager)->filename()<<endl;
117
118 // ---- primary ntuple ------
119 HepTuple* ntuple1 = hbookManager->ntuple("Primary Ntuple");
120 assert (ntuple1 != 0);
121
122 // ---- secondary ntuple ------
123 HepTuple* ntuple2 = hbookManager->ntuple("Secondary Ntuple");
124 assert (ntuple2 != 0);
125
126 // ---- table ntuple ------
127 HepTuple* ntuple3 = hbookManager->ntuple("Mean Free Path Ntuple");
128 assert (ntuple3 != 0);
129
130 // ---- secondaries histos ----
131 HepHistogram* hEKin;
132 hEKin = hbookManager->histogram("Kinetic Energy", 100,0.,200.);
133 assert (hEKin != 0);
134
135 HepHistogram* hP;
136 hP = hbookManager->histogram("Momentum", 100,0.,1000.);
137 assert (hP != 0);
138
139 HepHistogram* hNSec;
140 hNSec = hbookManager->histogram("Number of secondaries", 40,0.,40.);
141 assert (hNSec != 0);
142
143 HepHistogram* hDebug;
144 hDebug = hbookManager->histogram("Debug", 100,0.,200.);
145 assert (hDebug != 0);
146
147
148 //--------- Materials definition ---------
149
150 G4Material* Si = new G4Material("Silicon", 14., 28.055*g/mole, 2.33*g/cm3);
151 G4Material* Fe = new G4Material("Iron", 26., 55.85*g/mole, 7.87*g/cm3);
152 G4Material* Cu = new G4Material("Copper", 29., 63.55*g/mole, 8.96*g/cm3);
153 G4Material* W = new G4Material("Tungsten", 74., 183.85*g/mole, 19.30*g/cm3);
154 G4Material* Pb = new G4Material("Lead", 82., 207.19*g/mole, 11.35*g/cm3);
155 G4Material* U = new G4Material("Uranium", 92., 238.03*g/mole, 18.95*g/cm3);
156
157 G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
158 G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
159 G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
160 G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
161 G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole);
162
163 G4Material* maO = new G4Material("Oxygen", 8., 16.00*g/mole, 1.1*g/cm3);
164
165 G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2);
166 water->AddElement(H,2);
167 water->AddElement(O,1);
168
169 G4Material* ethane = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
170 ethane->AddElement(H,6);
171 ethane->AddElement(C,2);
172
173 G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2);
174 csi->AddElement(Cs,1);
175 csi->AddElement(I,1);
176
177
178 // Interactive set-up
179
180 G4cout << "How many interactions? " << G4endl;
181 G4cin >> nIterations;
182
183 if (nIterations <= 0) G4Exception("Wrong input");
184
185 G4double initEnergy = 1*MeV;
186 G4double initX = 0.;
187 G4double initY = 0.;
188 G4double initZ = 1.;
189
190 G4cout << "Enter the initial particle energy E (MeV)" << G4endl;
191 G4cin >> initEnergy ;
192
193 initEnergy = initEnergy * MeV;
194
195 if (initEnergy <= 0.) G4Exception("Wrong input");
196
197 static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
198
199 G4int nMaterials = G4Material::GetNumberOfMaterials();
200
201 G4cout << "Available materials are: " << G4endl;
202 for (G4int mat = 0; mat < nMaterials; mat++)
203 {
204 G4cout << mat << ") "
205 << (*theMaterialTable)[mat]->GetName()
206 << G4endl;
207 }
208
209 G4cout << "Which material? " << G4endl;
210 G4cin >> materialId;
211
212 G4Material* material = (*theMaterialTable)[materialId] ;
213
214 G4cout << "The selected material is: "
215 << material->GetName()
216 << G4endl;
217
218 G4double dimX = 1*mm;
219 G4double dimY = 1*mm;
220 G4double dimZ = 1*mm;
221
222 // Geometry
223
224 G4Box* theFrame = new G4Box ("Frame",dimX, dimY, dimZ);
225
226 G4LogicalVolume* logicalFrame = new G4LogicalVolume(theFrame,
227 (*theMaterialTable)[materialId],
228 "LFrame", 0, 0, 0);
229 logicalFrame->SetMaterial(material);
230
231 G4PVPlacement* physicalFrame = new G4PVPlacement(0,G4ThreeVector(),
232 "PFrame",logicalFrame,0,false,0);
233
234 // Particle definitions
235
236 G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
237 G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
238 G4ParticleDefinition* positron = G4Positron::PositronDefinition();
239
240 gamma->SetCuts(1*micrometer);
241 electron->SetCuts(1*micrometer);
242 positron->SetCuts(1*micrometer);
243
244 G4Gamma::SetEnergyRange(2.5e-4*MeV,1e5*MeV);
245 G4Electron::SetEnergyRange(2.5e-4*MeV,1e5*MeV);
246 G4Positron::SetEnergyRange(2.5e-4*MeV,1e5*MeV);
247
248 G4cout<<"the cut in energy for gamma in: "<<
249 (*theMaterialTable)[materialId]->GetName()
250 <<" is: "<<G4Gamma::GetCutsInEnergy()[materialId]<<G4endl;
251 G4cout<<"the cut in energy for e- in: "<<
252 (*theMaterialTable)[materialId]->GetName()
253 <<" is: "<<G4Electron::GetCutsInEnergy()[materialId]<<G4endl;
254
255 // Processes
256
257
258 G4int processType;
259 G4cout << "LowEnergy [1] or Standard [2] Gamma Conversion?" << G4endl;
260 G4cin >> processType;
261 if ( !(processType == 1 || processType == 2))
262 {
263 G4Exception("Wrong input");
264 }
265
266 G4VDiscreteProcess* gammaProcess;
267
268 if (processType == 1)
269 {
270 gammaProcess = new G4LowEnergyGammaConversion();
271 }
272 else
273 {
274 gammaProcess = new G4GammaConversion();
275 }
276
277
278 G4VProcess* theeminusMultipleScattering = new G4MultipleScattering();
279 G4VProcess* theeminusIonisation = new G4eIonisation();
280 G4VProcess* theeminusBremsstrahlung = new G4eBremsstrahlung();
281 G4VProcess* theeplusMultipleScattering = new G4MultipleScattering();
282 G4VProcess* theeplusIonisation = new G4eIonisation();
283 G4VProcess* theeplusBremsstrahlung = new G4eBremsstrahlung();
284 G4VProcess* theeplusAnnihilation = new G4eplusAnnihilation();
285
286 //----------------
287 // process manager
288 //----------------
289
290 // gamma
291
292 G4ProcessManager* gProcessManager = new G4ProcessManager(gamma);
293 gamma->SetProcessManager(gProcessManager);
294 gProcessManager->AddDiscreteProcess(gammaProcess);
295 G4ForceCondition* condition;
296
297 //electron
298
299 G4ProcessManager* eProcessManager = new G4ProcessManager(electron);
300 electron->SetProcessManager(eProcessManager);
301 eProcessManager->AddProcess(theeminusMultipleScattering);
302 eProcessManager->AddProcess(theeminusIonisation);
303 eProcessManager->AddProcess(theeminusBremsstrahlung);
304
305 //positron
306
307 G4ProcessManager* pProcessManager = new G4ProcessManager(positron);
308 positron->SetProcessManager(pProcessManager);
309 pProcessManager->AddProcess(theeplusMultipleScattering);
310 pProcessManager->AddProcess(theeplusIonisation);
311 pProcessManager->AddProcess(theeplusBremsstrahlung);
312 pProcessManager->AddProcess(theeplusAnnihilation);
313
314 //--------------
315 // set ordering
316 //--------------
317
318
319 eProcessManager->
320 SetProcessOrdering(theeminusMultipleScattering, idxAlongStep,1);
321 eProcessManager->
322 SetProcessOrdering(theeminusIonisation, idxAlongStep,2);
323
324 eProcessManager->
325 SetProcessOrdering(theeminusMultipleScattering, idxPostStep,1);
326 eProcessManager->
327 SetProcessOrdering(theeminusIonisation, idxPostStep,2);
328 eProcessManager->
329 SetProcessOrdering(theeminusBremsstrahlung, idxPostStep,3);
330
331
332
333 pProcessManager->SetProcessOrderingToFirst(theeplusAnnihilation, idxAtRest);
334 pProcessManager->
335 SetProcessOrdering(theeplusMultipleScattering, idxAlongStep,1);
336 pProcessManager->
337 SetProcessOrdering(theeplusIonisation, idxAlongStep,2);
338
339 pProcessManager->
340 SetProcessOrdering(theeplusMultipleScattering, idxPostStep,1);
341 pProcessManager->
342 SetProcessOrdering(theeplusIonisation, idxPostStep,2);
343 pProcessManager->
344 SetProcessOrdering(theeplusBremsstrahlung, idxPostStep,3);
345 pProcessManager->
346 SetProcessOrdering(theeplusAnnihilation, idxPostStep,4);
347
348 // G4LowEnergyIonisation IonisationProcess;
349 // eProcessManager->AddProcess(&IonisationProcess);
350 // eProcessManager->SetProcessOrdering(&IonisationProcess,idxAlongStep,1);
351 // eProcessManager->SetProcessOrdering(&IonisationProcess,idxPostStep, 1);
352
353 // G4LowEnergyBremsstrahlung BremstrahlungProcess;
354 // eProcessManager->AddProcess(&BremstrahlungProcess);
355 // eProcessManager->SetProcessOrdering(&BremstrahlungProcess,idxAlongStep,1);
356 // eProcessManager->SetProcessOrdering(&BremstrahlungProcess,idxPostStep, 1);
357
358 // G4eIonisation IonisationPlusProcess;
359 // pPositronProcessManager->AddProcess(&IonisationPlusProcess);
360 // pProcessManager->
361 // SetProcessOrdering(&IonisationPlusProcess,idxAlongStep,1);
362 // pProcessManager->SetProcessOrdering(&IonisationPlusProcess,idxPostStep,1);
363
364
365
366 // Create a DynamicParticle
367
368 G4double eEnergy = initEnergy*MeV;
369 G4ParticleMomentum eDirection(initX,initY,initZ);
370 G4DynamicParticle dynamicGamma(G4Gamma::Gamma(),eDirection,eEnergy);
371
372 dynamicGamma.DumpInfo(0);
373
374 // Track
375
376 G4ThreeVector aPosition(0.,0.,0.);
377 G4double aTime = 0. ;
378
379 G4Track* gTrack = new G4Track(&dynamicGamma,aTime,aPosition);
380
381 G4GRSVolume* touche = new G4GRSVolume(physicalFrame, NULL, aPosition);
382 gTrack->SetTouchable(touche);
383
384
385 // Step
386
387 G4Step* step = new G4Step();
388 step->SetTrack(gTrack);
389
390 G4StepPoint* aPoint = new G4StepPoint();
391 aPoint->SetPosition(aPosition);
392 aPoint->SetMaterial(material);
393 G4double safety = 10000.*cm;
394 aPoint->SetSafety(safety);
395 step->SetPreStepPoint(aPoint);
396
397 // Check applicability
398
399 if (! (gammaProcess->IsApplicable(*gamma)))
400 {
401 G4Exception("Not Applicable");
402 }
403 else
404 {
405 G4cout<< "applicability OK" << endl;
406 }
407
408 // Initialize the physics tables (in which material?)
409
410 gammaProcess->BuildPhysicsTable(*gamma);
411
412 theeminusMultipleScattering->BuildPhysicsTable(*electron);
413 theeminusIonisation->BuildPhysicsTable(*electron);
414 theeminusBremsstrahlung->BuildPhysicsTable(*electron);
415 theeplusMultipleScattering->BuildPhysicsTable(*positron);
416 theeplusIonisation->BuildPhysicsTable(*positron);
417 theeplusBremsstrahlung->BuildPhysicsTable(*positron);
418 theeplusAnnihilation->BuildPhysicsTable(*positron) ;
419
420 G4cout<< "table OK" << endl;
421
422 // Test GetMeanFreePath()
423
424 G4Material* apttoMaterial ;
425 G4String MaterialName ;
426
427 G4double minArg = 100*eV,maxArg = 100*GeV, argStp;
428 const G4int pntNum = 300;
429 G4double Tkin[pntNum+1];
430 G4double meanFreePath=0. ;
431
432 argStp = (std::log10(maxArg)-std::log10(minArg))/pntNum;
433
434 for(G4int d = 0; d < pntNum+1; d++)
435 {
436 Tkin[d] = std::pow(10,(std::log10(minArg) + d*argStp));
437 }
438
439 G4double sti = 1.*mm;
440 step->SetStepLength(sti);
441
442 // for ( G4int J = 0 ; J < nMaterials ; J++ )
443 // {
444 apttoMaterial = (*theMaterialTable)[materialId] ;
445 MaterialName = apttoMaterial->GetName() ;
446 logicalFrame->SetMaterial(apttoMaterial);
447
448 gTrack->SetStep(step);
449
450 G4LowEnergyGammaConversion* gammaLowEProcess =
451 (G4LowEnergyGammaConversion*) gammaProcess;
452 G4GammaConversion* gammaStdProcess =
453 (G4GammaConversion*) gammaProcess;
454
455
456 for (G4int i=0 ; i<pntNum; i++)
457 {
458 dynamicGamma.SetKineticEnergy(Tkin[i]);
459 if (processType == 1)
460 {
461 meanFreePath=gammaLowEProcess
462 ->GetMeanFreePath(*gTrack, sti, condition);
463 }
464 else
465 {
466 meanFreePath=gammaStdProcess
467 ->GetMeanFreePath(*gTrack, sti, condition);
468 }
469
470 ntuple3->column("kinen",Tkin[i]);
471 ntuple3->column("mfp",meanFreePath/cm);
472 ntuple3->dumpData();
473
474 // G4cout << meanFreePath/cm << G4endl;
475
476 }
477 G4cout << "Mean Free Path OK" << G4endl;
478
479 // --------- Test the DoIt
480
481 G4cout << "DoIt in " << material->GetName() << G4endl;
482
483
484 dynamicGamma.SetKineticEnergy(eEnergy);
485 for (G4int iter=0; iter<nIterations; iter++)
486 {
487
488 step->SetStepLength(1*micrometer);
489
490 G4cout << "Iteration = " << iter
491 << " - Step Length = "
492 << step->GetStepLength()/mm << " mm "
493 << G4endl;
494
495 gTrack->SetStep(step);
496
497 // G4cout << "Iteration = " << iter
498 // << " - Step Length = "
499 // << step->GetStepLength()/mm << " mm "
500 // << G4endl;
501
502 //G4cout << gTrack->GetStep()->GetStepLength()/mm
503 // << G4endl;
504
505 G4VParticleChange* dummy;
506 dummy = gammaProcess->PostStepDoIt(*gTrack, *step);
507
508 G4ParticleChange* particleChange = (G4ParticleChange*) dummy;
509
510 // Primary physical quantities
511
512 G4double energyChange = particleChange->GetEnergyChange();
513
514 G4double dedx = initEnergy - energyChange ;
515 G4double dedxNow = dedx / (step->GetStepLength());
516
517 G4ThreeVector eChange =
518 particleChange->CalcMomentum(energyChange,
519 (*particleChange->GetMomentumChange()),
520 particleChange->GetMassChange());
521
522 G4double pxChange = eChange.x();
523 G4double pyChange = eChange.y();
524 G4double pzChange = eChange.z();
525 G4double pChange =
526 std::sqrt(pxChange*pxChange + pyChange*pyChange + pzChange*pzChange);
527
528 G4double xChange = particleChange->GetPositionChange()->x();
529 G4double yChange = particleChange->GetPositionChange()->y();
530 G4double zChange = particleChange->GetPositionChange()->z();
531 G4double thetaChange = particleChange->GetPositionChange()->theta();
532
533 G4cout << "---- Primary after the step ---- " << G4endl;
534
535 // G4cout << "Position (x,y,z) = "
536 // << xChange << " "
537 // << yChange << " "
538 // << zChange << " "
539 // << G4endl;
540
541 G4cout << "---- Energy: " << energyChange/MeV << " MeV, "
542 << "(px,py,pz): ("
543 << pxChange/MeV << ","
544 << pyChange/MeV << ","
545 << pzChange/MeV << ") MeV"
546 << G4endl;
547
548 G4cout << "---- Energy loss (dE) = " << dedx/keV << " keV" << G4endl;
549 // G4cout << "Stopping power (dE/dx)=" << dedxNow << G4endl;
550
551 // Secondaries
552
553 ntuple1->column("eprimary", initEnergy);
554 ntuple1->column("energyf", energyChange);
555 ntuple1->column("de", dedx);
556 ntuple1->column("dedx", dedxNow);
557 ntuple1->column("pxch", xChange);
558 ntuple1->column("pych", pyChange);
559 ntuple1->column("pzch", pzChange);
560 ntuple1->column("pch", zChange);
561 ntuple1->column("thetach", thetaChange);
562 ntuple1->dumpData();
563
564 // Secondaries physical quantities
565
566 hNSec->accumulate(particleChange->GetNumberOfSecondaries());
567 hDebug->accumulate(particleChange->GetLocalEnergyDeposit());
568
569 G4cout << " secondaries " <<
570 particleChange->GetNumberOfSecondaries() << G4endl;
571
572 for (G4int i = 0; i < (particleChange->GetNumberOfSecondaries()); i++)
573 {
574 // The following two items should be filled per event, not
575 // per secondary; filled here just for convenience, to avoid
576 // complicated logic to dump ntuple when there are no secondaries
577
578 G4Track* finalParticle = particleChange->GetSecondary(i) ;
579
580 G4double e = finalParticle->GetTotalEnergy();
581 G4double eKin = finalParticle->GetKineticEnergy();
582 G4double px = (finalParticle->GetMomentum()).x();
583 G4double py = (finalParticle->GetMomentum()).y();
584 G4double pz = (finalParticle->GetMomentum()).z();
585 G4double theta = (finalParticle->GetMomentum()).theta();
586 G4double p = std::sqrt(px*px+py*py+pz*pz);
587
588 if (e > initEnergy)
589 {
590 G4cout << "WARNING: eFinal > eInit " << G4endl;
591 // << e
592 // << " > " initEnergy
593
594 }
595
596 G4String particleName =
597 finalParticle->GetDefinition()->GetParticleName();
598 G4cout << "==== Final "
599 << particleName << " "
600 << "energy: " << e/MeV << " MeV, "
601 << "eKin: " << eKin/MeV << " MeV, "
602 << "(px,py,pz): ("
603 << px/MeV << ","
604 << py/MeV << ","
605 << pz/MeV << ") MeV "
606 << G4endl;
607
608 hEKin->accumulate(eKin);
609 hP->accumulate(p);
610
611 G4int partType;
612 if (particleName == "e-") partType = 1;
613 else if (particleName == "e+") partType = 2;
614 else if (particleName == "gamma") partType = 3;
615
616 // Fill the secondaries ntuple
617
618 ntuple2->column("event",iter);
619 ntuple2->column("eprimary",initEnergy);
620 ntuple2->column("px", px);
621 ntuple2->column("py", py);
622 ntuple2->column("pz", pz);
623 ntuple2->column("p", p);
624 ntuple2->column("e", e);
625 ntuple2->column("theta", theta);
626 ntuple2->column("ekin", eKin);
627 ntuple2->column("type", partType);
628
629 ntuple2->dumpData();
630
631 delete particleChange->GetSecondary(i);
632 }
633
634 particleChange->Clear();
635
636 }
637
638
639 G4cout << "Iteration number: " << iter << G4endl;
640 hbookManager->write();
641 delete hbookManager;
642
643 delete step;
644
645
646 cout << "END OF THE MAIN PROGRAM" << G4endl;
647}
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
Note: See TracBrowser for help on using the repository browser.