source: trunk/source/processes/electromagnetic/lowenergy/test/G4LowEnergyTest.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: 21.5 KB
RevLine 
[1350]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
83HepTupleManager* hbookManager;
84
85main()
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.01*g/mole, 1.848*g/cm3);
196 G4Material* Graphite = new G4Material("Graphite",6., 12.00*g/mole, 2.265*g/cm3 );
197 G4Material* Al = new G4Material("Aluminium", 13., 26.98*g/mole, 2.7 *g/cm3);
198 G4Material* LAr = new G4Material("LArgon", 18., 39.95*g/mole, 1.393*g/cm3);
199 G4Material* Fe = new G4Material("Iron", 26., 55.85*g/mole, 7.87*g/cm3);
200 G4Material* Cu = new G4Material("Copper", 29., 63.55*g/mole, 8.96*g/cm3);
201 G4Material* W = new G4Material("Tungsten", 74., 183.85*g/mole, 19.30*g/cm3);
202 G4Material* Pb = new G4Material("Lead", 82., 207.19*g/mole, 11.35*g/cm3);
203 G4Material* U = new G4Material("Uranium", 92., 238.03*g/mole, 18.95*g/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.00*g/mole, 1.1*g/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",92*mm, 92*mm, 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-4*MeV,1e5*MeV);
297 G4Electron::SetEnergyRange(2.5e-4*MeV,1e5*MeV);
298 G4Positron::SetEnergyRange(2.5e-4*MeV,1e5*MeV);
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 = 100*eV,maxArg = 100*GeV, 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(pxChange*pxChange+pyChange*pyChange+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(Px*Px+Py*Py+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(Px*Px+Py*Py+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(Px*Px+Py*Py+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.