source: trunk/source/processes/electromagnetic/lowenergy/test/G4ComplexTest.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: 24.6 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: G4ComplexTest.cc,v 1.21 2006/06/29 19:43:49 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: G4ComplexTest
36// This test provide AlongStepDoIt and PostStepDoIt
37// tests for electromagnetic processes. The input
38// data have to be describe in ASCII file
39//
40// Author: V.Ivanchenko on base of Maria Grazia Pia tests
41//
42// Creation date: 8 May 2001
43//
44// Modifications:
45//
46// -------------------------------------------------------------------
47
48#include "globals.hh"
49#include "G4ios.hh"
50#include <fstream>
51#include <iomanip>
52
53#include "G4ProductionCuts.hh"
54#include "G4ProductionCutsTable.hh"
55#include "G4Material.hh"
56#include "G4MaterialCutsCouple.hh"
57#include "G4VContinuousDiscreteProcess.hh"
58#include "G4ProcessManager.hh"
59
60#include "G4LowEnergyIonisation.hh"
61#include "G4LowEnergyBremsstrahlung.hh"
62#include "G4LowEnergyCompton.hh"
63#include "G4LowEnergyGammaConversion.hh"
64#include "G4LowEnergyPhotoElectric.hh"
65#include "G4LowEnergyRayleigh.hh"
66
67#include "G4hLowEnergyIonisation.hh"
68
69#include "G4eIonisation.hh"
70#include "G4eBremsstrahlung.hh"
71#include "G4ComptonScattering.hh"
72#include "G4GammaConversion.hh"
73#include "G4PhotoElectricEffect.hh"
74
75#include "G4eplusAnnihilation.hh"
76
77#include "G4MuIonisation.hh"
78#include "G4MuBremsstrahlung.hh"
79#include "G4MuPairProduction.hh"
80
81#include "G4hIonisation.hh"
82
83#include "G4MultipleScattering.hh"
84
85#include "G4EnergyLossTables.hh"
86#include "G4ParticleChange.hh"
87#include "G4ParticleChange.hh"
88#include "G4DynamicParticle.hh"
89#include "G4AntiProton.hh"
90#include "G4Proton.hh"
91#include "G4Electron.hh"
92#include "G4Positron.hh"
93#include "G4Gamma.hh"
94
95#include "G4Box.hh"
96#include "G4PVPlacement.hh"
97#include "G4VPhysicalVolume.hh"
98#include "G4LogicalVolume.hh"
99#include "G4RunManager.hh"
100
101#include "G4Step.hh"
102#include "G4GRSVolume.hh"
103
104#include "G4UnitsTable.hh"
105
106// New Histogramming (from AIDA and Anaphe):
107#include <memory> // for the auto_ptr(T>
108
109#include "AIDA/AIDA.h"
110
111
112#include "G4Timer.hh"
113
114int main(int argc,char** argv)
115{
116
117 // -------------------------------------------------------------------
118 // Setup
119
120 G4int nEvt = 100;
121 G4int nPart =-1;
122 G4String nameMat = "Si";
123 G4int nProcess = 0;
124 G4bool usepaw = false;
125 G4bool postDo = true;
126 G4bool lowE = true;
127 G4int verbose = 0;
128 G4double gEnergy = 0.1*MeV;
129 G4String hFile = "";
130 G4double theStep = 1.0*micrometer;
131 G4double range = 1.0*micrometer;
132 G4double cutG = 1.0*micrometer;
133 G4double cutE = 1.0*micrometer;
134 G4Material* material = 0;
135 G4String name[6] = {"Ionisation", "Bremsstrahlung", "Compton",
136 "GammaConversion", "PhotoElectric", "Raylaigh"};
137
138
139 G4cout.setf( ios::scientific, ios::floatfield );
140
141 // -------------------------------------------------------------------
142 // Control on input
143
144 if(argc < 2) {
145 G4cout << "Input file is not specified! Exit" << G4endl;
146 exit(1);
147 }
148
149 ifstream* fin = new ifstream();
150 string fname = argv[1];
151 fin->open(fname.c_str());
152 if( !fin->is_open()) {
153 G4cout << "Input file <" << fname << "> does not exist! Exit" << G4endl;
154 exit(1);
155 }
156
157 // -------------------------------------------------------------------
158 //--------- Materials definition ---------
159
160 G4Material* m;
161 m = new G4Material("Be", 4., 9.01*g/mole, 1.848*g/cm3);
162 m = new G4Material("Graphite",6., 12.00*g/mole, 2.265*g/cm3 );
163 m = new G4Material("Al", 13., 26.98*g/mole, 2.7 *g/cm3);
164 m = new G4Material("Si", 14., 28.055*g/mole, 2.33*g/cm3);
165 m = new G4Material("LAr", 18., 39.95*g/mole, 1.393*g/cm3);
166 m = new G4Material("Fe", 26., 55.85*g/mole, 7.87*g/cm3);
167 m = new G4Material("Cu", 29., 63.55*g/mole, 8.96*g/cm3);
168 m = new G4Material("W", 74., 183.85*g/mole, 19.30*g/cm3);
169 m = new G4Material("Pb", 82., 207.19*g/mole, 11.35*g/cm3);
170 m = new G4Material("U", 92., 238.03*g/mole, 18.95*g/cm3);
171
172 G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
173 G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
174 G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
175 G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
176 G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole);
177
178 m = new G4Material("O2", 8., 16.00*g/mole, 1.1*g/cm3);
179
180 m = new G4Material ("Water" , 1.*g/cm3, 2);
181 m->AddElement(H,2);
182 m->AddElement(O,1);
183
184 m = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
185 m->AddElement(H,6);
186 m->AddElement(C,2);
187
188 m = new G4Material ("CsI" , 4.53*g/cm3, 2);
189 m->AddElement(Cs,1);
190 m->AddElement(I,1);
191
192 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
193
194 G4int nMaterials = G4Material::GetNumberOfMaterials();
195
196 G4cout << "Available materials are: " << G4endl;
197 G4int mat;
198 for (mat = 0; mat < nMaterials; mat++) {
199 G4cout << mat << ") " << (*theMaterialTable)[mat]->GetName() << G4endl;
200 }
201 material = (*theMaterialTable)[0];
202
203 G4cout << "Available processes are: " << G4endl;
204 for (mat = 0; mat < 6; mat++) {
205 G4cout << mat << ") " << name[mat] << G4endl;
206 }
207
208 // Particle definitions
209
210 G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
211 G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
212 G4ParticleDefinition* positron = G4Positron::PositronDefinition();
213 G4ParticleDefinition* proton = G4Proton::ProtonDefinition();
214 G4ParticleDefinition* antiproton = G4AntiProton::AntiProtonDefinition();
215 G4ParticleDefinition* part = gamma;
216
217 // Geometry
218
219 G4double initX = 0.;
220 G4double initY = 0.;
221 G4double initZ = 1.;
222 G4double dimX = 100.0*cm;
223 G4double dimY = 100.0*cm;
224 G4double dimZ = 100.0*cm;
225
226 G4Box* sFrame = new G4Box ("Box",dimX, dimY, dimZ);
227 G4LogicalVolume* lFrame = new G4LogicalVolume(sFrame,material,"Box",0,0,0);
228 G4PVPlacement* pFrame = new G4PVPlacement(0,G4ThreeVector(),"Box",
229 lFrame,0,false,0);
230 G4RunManager* rm = new G4RunManager();
231 G4cout << "World is defined " << G4endl;
232 rm->GeometryHasBeenModified();
233 rm->DefineWorldVolume(pFrame);
234
235 // -------------------------------------------------------------------
236 // ---- Read input file
237 G4cout << "Available commands are: " << G4endl;
238 G4cout << "#events" << G4endl;
239 G4cout << "#particle" << G4endl;
240 G4cout << "#energy(MeV)" << G4endl;
241 G4cout << "#cutG(mm)" << G4endl;
242 G4cout << "#cutE(mm)" << G4endl;
243 G4cout << "#range(mm)" << G4endl;
244 G4cout << "#step(mm)" << G4endl;
245 G4cout << "#material" << G4endl;
246 G4cout << "#process" << G4endl;
247 G4cout << "#domain" << G4endl;
248 G4cout << "#test" << G4endl;
249 G4cout << "#paw" << G4endl;
250 G4cout << "#verbose" << G4endl;
251 G4cout << "#run" << G4endl;
252 G4cout << "#exit" << G4endl;
253 G4cout << pFrame << G4endl;
254
255 G4ProcessManager *elecManager, *positManager, *gammaManager,
256 *protManager, *aprotManager;
257
258 elecManager = new G4ProcessManager(electron);
259 electron->SetProcessManager(elecManager);
260
261 positManager = new G4ProcessManager(positron);
262 positron->SetProcessManager(positManager);
263
264 gammaManager = new G4ProcessManager(gamma);
265 gamma->SetProcessManager(gammaManager);
266
267 protManager = new G4ProcessManager(proton);
268 proton->SetProcessManager(protManager);
269
270 aprotManager = new G4ProcessManager(antiproton);
271 antiproton->SetProcessManager(aprotManager);
272
273 G4eIonisation* elecSTion = 0;
274 G4eBremsstrahlung* elecSTbr = 0;
275 G4LowEnergyIonisation* elecLEion = 0;
276 G4LowEnergyBremsstrahlung* elecLEbr = 0;
277 G4bool ionis = true;
278
279 string line, line1;
280 G4bool end = true;
281 for(G4int run=0; run<100; run++) {
282 do {
283 (*fin) >> line;
284 G4cout << "Next line " << line << G4endl;
285 if(line == "#events") {
286 (*fin) >> nEvt;
287 if(nEvt < 1) nEvt = 1;
288 } else if(line == "#particle") {
289 (*fin) >> nPart;
290 } else if(line == "#energy(MeV)") {
291 (*fin) >> gEnergy;
292 gEnergy *= MeV;
293 } else if(line == "#cutG(mm)") {
294 (*fin) >> cutG;
295 cutG *= mm;
296 } else if(line == "#cutE(mm)") {
297 (*fin) >> cutE;
298 cutE *= mm;
299 } else if(line == "#range(mm)") {
300 (*fin) >> range;
301 range *= mm;
302 } else if(line == "#step(mm)") {
303 (*fin) >> theStep;
304 theStep *= mm;
305 } else if(line == "#material") {
306 (*fin) >> nameMat;
307 } else if(line == "#process") {
308 (*fin) >> nProcess;
309 } else if(line == "#domain") {
310 (*fin) >> line1;
311 if(line1 == "lowenergy") {lowE = true;}
312 else {lowE = false;}
313 } else if(line == "#test") {
314 (*fin) >> line1;
315 if(line1 == "PostStep") {postDo = true;}
316 else {postDo = false;}
317 } else if(line == "#paw") {
318 usepaw = true;
319 (*fin) >> hFile;
320 } else if(line == "#run") {
321 break;
322 } else if(line == "#verbose") {
323 (*fin) >> verbose;
324 } else if(line == "#exit") {
325 end = false;
326 break;
327 }
328 } while(end);
329
330 if(!end) break;
331
332 G4cout << "###### Start new run # " << run << " #####" << G4endl;
333 if(nPart == 0) {
334 part = gamma;
335 } else if(nPart == 1) {
336 part = electron;
337 } else if(nPart == 2) {
338 part = positron;
339 } else if(nPart == 3) {
340 part = proton;
341 } else if(nPart == 4) {
342 part = antiproton;
343 } else {
344 G4cout << "Particle #" << nPart
345 << " is absent in the list of particles: Exit" << G4endl;
346 break;
347 }
348 if(nProcess < 0 || nProcess > 5) {
349 G4cout << "Process #" << nProcess
350 << " is absent in the list of processes: Exit" << G4endl;
351 break;
352 }
353
354 for (mat = 0; mat < nMaterials; mat++) {
355 material = (*theMaterialTable)[mat];
356 if(nameMat == material->GetName()) break;
357 }
358
359 G4cout << "The particle: " << part->GetParticleName() << G4endl;
360 G4cout << "The energy: " << gEnergy/MeV << " MeV" << G4endl;
361 G4cout << "The material: " << material->GetName() << G4endl;
362 G4cout << "The cut on e-:" << cutE/mm << " mm" << G4endl;
363 G4cout << "The cut on g: " << cutG/mm << " mm" << G4endl;
364 G4cout << "The step: " << theStep/mm << " mm" << G4endl;
365 if(postDo && lowE) {
366 G4cout << "Test of PostStepDoIt for " << name[nProcess]
367 << " for lowenergy" << G4endl;
368 } else if(postDo && !lowE) {
369 G4cout << "Test of PostStepDoIt for " << name[nProcess]
370 << " for standard" << G4endl;
371 } else if(!postDo && !lowE) {
372 G4cout << "Test of AlongStepDoIt for " << name[nProcess]
373 << " for standard" << G4endl;
374 } else if(!postDo && lowE) {
375 G4cout << "Test of AlongStepDoIt for " << name[nProcess]
376 << " for lowenergy" << G4endl;
377 }
378
379 // -------------------------------------------------------------------
380 // ---- HBOOK initialization
381
382 // Creating the analysis factory
383 std::auto_ptr< AIDA::IAnalysisFactory > af( AIDA_createAnalysisFactory() );
384
385 // Creating the tree factory
386 std::auto_ptr< AIDA::ITreeFactory > tf( af->createTreeFactory() );
387
388 // Creating a tree mapped to a new hbook file.
389 std::auto_ptr< AIDA::ITree > tree( tf->create( hFile,"hbook" ,false,false ) );
390 G4cout << "Tree store : " << tree->storeName() << G4endl;
391
392 // Creating a tuple factory, whose tuples will be handled by the tree
393 // std::auto_ptr< AIDA::ITupleFactory > tpf( af->createTupleFactory( *tree ) );
394
395 AIDA::IHistogram1D* hist[4];
396 //AIDA::ITuple* ntuple1 = 0;
397 //AIDA::ITuple* ntuple2 = 0;
398
399 if(usepaw) {
400
401 // ---- primary ntuple ------
402 // If using Anaphe HBOOK implementation, there is a limitation on the length of the
403 // variable names in a ntuple
404 //ntuple1 = tpf->create( "100", "Primary", "float ekin, dedx" );
405 //ntuple2 = tpf->create( "101", "Secondary", "float ekin, dedx" );
406
407
408 // Creating a histogram factory, whose histograms will be handled by the tree
409 std::auto_ptr< AIDA::IHistogramFactory > hf( af->createHistogramFactory( *tree ) );
410
411 // Creating an 1-dimensional histogram in the root directory of the tree
412
413 hist[0] = hf->createHistogram1D("11","Kinetic Energy (T/T0)", 50,0.,1.0);
414 hist[1] = hf->createHistogram1D("12","Momentum (MeV/c)", 50,0.,gEnergy*0.1/MeV);
415 hist[2] = hf->createHistogram1D("13","Number of secondaries", 20,-0.5,19.5);
416 hist[3] = hf->createHistogram1D("14","Energy deposition (MeV)", 50,0.,gEnergy*0.1/MeV);
417
418 G4cout<< "Histograms is initialised" << G4endl;
419 }
420
421 G4Timer* timer = new G4Timer();
422 timer->Start();
423 G4ProductionCutsTable* cutsTable = G4ProductionCutsTable::GetProductionCutsTable();
424
425 G4ProductionCuts* cuts = cutsTable->GetDefaultProductionCuts();
426 cuts->SetProductionCut(cutG, 0);
427 cuts->SetProductionCut(cutE, 1);
428 cuts->SetProductionCut(cutE, 2);
429 G4cout << "Cuts are defined " << G4endl;
430
431
432 cutsTable->UpdateCoupleTable();
433 (G4ProductionCutsTable::GetProductionCutsTable())->DumpCouples();
434 const G4MaterialCutsCouple* theCouple = cutsTable->GetMaterialCutsCouple(material,cuts);
435 // Processes
436
437 G4VDiscreteProcess* dProcess;
438 G4VContinuousDiscreteProcess* cdProcess;
439 dProcess = 0;
440 cdProcess = 0;
441
442 G4cout << "Start BuildPhysicsTable" << G4endl;;
443
444 if(lowE) {
445 if(ionis) {
446 elecLEion = new G4LowEnergyIonisation();
447 elecLEbr = new G4LowEnergyBremsstrahlung();
448 elecManager->AddProcess(elecLEion);
449 elecManager->AddProcess(elecLEbr);
450 elecLEion->BuildPhysicsTable(*electron);
451 elecLEbr->BuildPhysicsTable(*electron);
452 ionis = false;
453 }
454 if(nPart == 0) {
455 if(nProcess == 2) dProcess = new G4LowEnergyCompton();
456 if(nProcess == 3) dProcess = new G4LowEnergyGammaConversion();
457 if(nProcess == 4) dProcess = new G4LowEnergyPhotoElectric();
458 if(nProcess == 5) dProcess = new G4LowEnergyRayleigh();
459 if(dProcess) {
460 gammaManager->AddProcess(dProcess);
461 dProcess->BuildPhysicsTable(*gamma);
462 }
463 } else if(nPart == 1) {
464 if(nProcess == 0) cdProcess = elecLEion;
465 if(nProcess == 1) cdProcess = elecLEbr;
466 } else if(nPart == 3) {
467 if(nProcess == 0) {
468 cdProcess = new G4hLowEnergyIonisation();
469 protManager->AddProcess(cdProcess);
470 cdProcess->BuildPhysicsTable(*proton);
471 }
472 } else if(nPart == 4) {
473 if(nProcess == 0) {
474 cdProcess = new G4hLowEnergyIonisation();
475 aprotManager->AddProcess(cdProcess);
476 cdProcess->BuildPhysicsTable(*antiproton);
477 }
478 }
479
480 } else {
481 if(ionis) {
482 elecSTion = new G4eIonisation();
483 elecSTbr = new G4eBremsstrahlung();
484 elecManager->AddProcess(elecSTion);
485 elecManager->AddProcess(elecSTbr);
486 elecSTion->BuildPhysicsTable(*electron);
487 elecSTbr->BuildPhysicsTable(*electron);
488 ionis = false;
489 }
490 if(nPart == 0) {
491 if(nProcess == 2) dProcess = new G4ComptonScattering();
492 if(nProcess == 3) dProcess = new G4GammaConversion();
493 if(nProcess == 4) dProcess = new G4PhotoElectricEffect();
494 if(dProcess) {
495 gammaManager->AddProcess(dProcess);
496 dProcess->BuildPhysicsTable(*gamma);
497 }
498 } else if(nPart == 1) {
499 if(nProcess == 0) cdProcess = elecSTion;
500 if(nProcess == 1) cdProcess = elecSTbr;
501 } else if(nPart == 2) {
502 G4eIonisation* pSTion = new G4eIonisation();
503 G4eBremsstrahlung* pSTbr = new G4eBremsstrahlung();
504 if(nProcess == 0) cdProcess = pSTion;
505 if(nProcess == 1) cdProcess = pSTbr;
506 } else if(nPart == 3) {
507 if(nProcess == 0) {
508 cdProcess = new G4hIonisation();
509 protManager->AddProcess(cdProcess);
510 cdProcess->BuildPhysicsTable(*proton);
511 }
512 } else if(nPart == 4) {
513 if(nProcess == 0) {
514 cdProcess = new G4hIonisation();
515 aprotManager->AddProcess(cdProcess);
516 cdProcess->BuildPhysicsTable(*antiproton);
517 }
518 }
519 }
520
521 G4cout << "Physics tables are built" << G4endl;;
522
523 timer->Stop();
524 G4cout << " " << *timer << G4endl;
525 delete timer;
526
527
528 // Control on processes
529 if(postDo && !dProcess && !cdProcess) {
530 G4cout << "Discret Process is not found out! Exit" << G4endl;
531 break;
532 }
533 if(!postDo && !cdProcess) {
534 G4cout << "Continues Discret Process is not found out! Exit" << G4endl;
535 break;
536 }
537
538 // Create a DynamicParticle
539
540 G4ParticleMomentum gDir(initX,initY,initZ);
541 G4DynamicParticle dParticle(part,gDir,gEnergy);
542
543 // Track
544 G4ThreeVector aPosition(0.,0.,0.);
545 G4double aTime = 0. ;
546
547 G4Track* gTrack;
548 gTrack = new G4Track(&dParticle,aTime,aPosition);
549
550 // Step
551
552 G4Step* step;
553 step = new G4Step();
554 step->SetTrack(gTrack);
555
556 G4StepPoint *aPoint, *bPoint;
557 aPoint = new G4StepPoint();
558 aPoint->SetPosition(aPosition);
559 aPoint->SetMaterial(material);
560 aPoint->SetMaterialCutsCouple(theCouple);
561 G4double safety = 10000.*cm;
562 aPoint->SetSafety(safety);
563 step->SetPreStepPoint(aPoint);
564
565 bPoint = aPoint;
566 G4ThreeVector bPosition(0.,0.,theStep);
567 bPoint->SetPosition(bPosition);
568 step->SetPostStepPoint(bPoint);
569 step->SetStepLength(theStep);
570
571 // --------- Test the DoIt
572 G4int nElectrons = 0;
573 G4int nPositrons = 0;
574 G4int nPhotons = 0;
575 G4double rmax = 0.0;
576 G4double de = 0.0;
577 G4double de2 = 0.0;
578
579 G4cout << "dProcess= " << dProcess << " cdProcess= " << cdProcess;
580 if(dProcess) G4cout << " name= " << dProcess->GetProcessName() << G4endl;
581 if(cdProcess) G4cout << " name= " << cdProcess->GetProcessName() << G4endl;
582
583
584 timer = new G4Timer();
585 timer->Start();
586
587 for (G4int iter=0; iter<nEvt; iter++) {
588
589 gTrack->SetStep(step);
590
591 if(verbose) {
592 G4cout << "Iteration = " << iter
593 << " - Step Length = "
594 << step->GetStepLength()/mm << " mm "
595 << G4endl;
596 }
597
598 G4VParticleChange* dummy = 0;
599 if(postDo) {
600 if(dProcess) dummy = dProcess->PostStepDoIt(*gTrack, *step);
601 if(cdProcess) {
602 dummy = cdProcess->PostStepDoIt(*gTrack, *step);
603 }
604 } else {
605 dummy = cdProcess->AlongStepDoIt(*gTrack, *step);
606 }
607 G4ParticleChange* particleChange = (G4ParticleChange*) dummy;
608
609 // Primary physical quantities
610
611 G4double energyChange = particleChange->GetEnergyChange();
612 G4double deltaE = gEnergy - energyChange ;
613 G4double dedx = deltaE / (step->GetStepLength());
614
615 G4ThreeVector change = particleChange->CalcMomentum(energyChange,
616 *(particleChange->GetMomentumChange()),
617 part->GetPDGMass());
618
619 G4double pxChange = change.x();
620 G4double pyChange = change.y();
621 G4double pzChange = change.z();
622
623 if(verbose) {
624 G4cout << "---- Primary after the step ---- " << G4endl;
625 G4cout << "---- Energy: " << energyChange/MeV << " MeV, "
626 << "(px,py,pz): ("
627 << pxChange/MeV << ","
628 << pyChange/MeV << ","
629 << pzChange/MeV << ") MeV"
630 << G4endl;
631 G4cout << "---- Energy loss (dE) = " << deltaE/keV << " keV;"
632 << "Stopping power (dE/dx)=" << dedx*mm/keV << " keV/mm"
633 << "; rmax(mm)= " << rmax << G4endl;
634 }
635
636 // Primary
637
638 G4int nsec = particleChange->GetNumberOfSecondaries();
639 /*
640 if(ntuple1) {
641 ntuple1->column("epri", gEnergy/MeV);
642 ntuple1->column("efin", energyChange/MeV);
643 ntuple1->column("dedx", dedx*mm/MeV);
644 ntuple1->column("nsec", nsec);
645 ntuple1->column("nele", nElectrons);
646 ntuple1->column("npho", nPhotons);
647 ntuple1->dumpData();
648 }
649 */
650 de += deltaE;
651 de2 += deltaE*deltaE;
652
653 // Secondaries physical quantities
654
655 if(usepaw) {
656 hist[2]->fill((float)nsec, (float)1.0);
657 hist[3]->fill((float)(particleChange->GetLocalEnergyDeposit()/MeV), (float)1.0);
658 }
659
660 for (G4int i = 0; i<nsec; i++) {
661 // The following two items should be filled per event, not
662 // per secondary; filled here just for convenience, to avoid
663 // complicated logic to dump ntuple when there are no secondaries
664
665 G4Track* finalParticle = particleChange->GetSecondary(i) ;
666
667 G4double e = finalParticle->GetTotalEnergy();
668 G4double eKin = finalParticle->GetKineticEnergy();
669 G4double px = (finalParticle->GetMomentum()).x();
670 G4double py = (finalParticle->GetMomentum()).y();
671 G4double pz = (finalParticle->GetMomentum()).z();
672 //G4double theta= (finalParticle->GetMomentum()).theta();
673 G4double p = std::sqrt(px*px + py*py + pz*pz);
674
675 if (eKin > gEnergy) {
676 G4cout << "WARNING: eFinal > eInit in event #" << iter << G4endl;
677 }
678
679 G4String partName = finalParticle->GetDefinition()->GetParticleName();
680 if(verbose) {
681 G4cout << "==== Final "
682 << partName << " "
683 << "E= " << e/MeV << " MeV, "
684 << "eKin: " << eKin/MeV << " MeV, "
685 << "(px,py,pz): ("
686 << px/MeV << ","
687 << py/MeV << ","
688 << pz/MeV << ") MeV,"
689 << " p= " << p << " MeV"
690 << G4endl;
691 }
692
693 if(usepaw) {
694 hist[0]->fill((float)(eKin/gEnergy), (float)1.0);
695 hist[1]->fill((float)(p/MeV), (float)1.0);
696 }
697
698 G4int partType = 0;
699 if (partName == "e-") {
700 partType = 1;
701 nElectrons++;
702
703 } else if (partName == "e+") {
704 partType = 2;
705 nPositrons++;
706
707 } else if (partName == "gamma") {
708 partType = 0;
709 nPhotons++;
710 }
711
712 // Fill the secondaries ntuple
713 /*
714 if(ntuple2) {
715 ntuple2->column("eprimary",gEnergy);
716 ntuple2->column("px", px);
717 ntuple2->column("py", py);
718 ntuple2->column("pz", pz);
719 ntuple2->column("p", p);
720 ntuple2->column("e", e);
721 ntuple2->column("theta", theta);
722 ntuple2->column("ekin", eKin);
723 ntuple2->column("type", partType);
724 ntuple2->dumpData();
725 }
726 */
727 delete particleChange->GetSecondary(i);
728 }
729
730 particleChange->Clear();
731
732 }
733 G4cout << "###### Statistics:" << G4endl;
734 G4cout << "Average number of secondary electrons= "
735 << (G4double)nElectrons/(G4double)nEvt << G4endl;
736 G4cout << "Average number of secondary positrons= "
737 << (G4double)nPositrons/(G4double)nEvt << G4endl;
738 G4cout << "Average number of secondary photons= "
739 << (G4double)nPhotons/(G4double)nEvt << G4endl;
740 G4double x = de/(G4double)nEvt;
741 G4double y = de2/(G4double)nEvt - x*x;
742 if(0.0 < y) y = std::sqrt(y);
743 G4cout << "Average energy deposition(MeV)= "
744 << x/MeV << " +- " << y/MeV << G4endl;
745
746 timer->Stop();
747 G4cout << " " << *timer << G4endl;
748 delete timer;
749
750 if(usepaw) {
751 tree->commit();
752 std::cout << "Closing the tree..." << std::endl;
753 tree->close();
754 G4cout << "# hbook is writed" << G4endl;
755 }
756
757 G4cout << "###### End of run # " << run << " ######" << G4endl;
758
759 } while(end);
760 G4cout << "###### End of test #####" << G4endl;
761}
762
763
764
765
766
767
768
769
770
771
772
773
Note: See TracBrowser for help on using the repository browser.