source: trunk/source/processes/electromagnetic/lowenergy/test/G4StoppingPowerTest.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: 27.2 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//
28// -------------------------------------------------------------------
29// GEANT 4 class file --- Copyright CERN 1998
30// CERN Geneva Switzerland
31//
32//
33// File name: G4StoppingPowerTest
34// This test provide cross sections
35// tests for electromagnetic processes. The input
36// data have to be describe in ASCII file
37//
38// Author: V.Ivanchenko
39//
40// Creation date: 23 May 2001
41//
42// Modifications:
43//
44// -------------------------------------------------------------------
45
46#include "globals.hh"
47#include "G4ios.hh"
48#include <fstream>
49#include <iomanip>
50
51#include "G4Material.hh"
52#include "G4VContinuousDiscreteProcess.hh"
53#include "G4ProcessManager.hh"
54#include "G4VParticleChange.hh"
55
56#include "G4LowEnergyIonisation.hh"
57#include "G4LowEnergyBremsstrahlung.hh"
58#include "G4LowEnergyCompton.hh"
59#include "G4LowEnergyGammaConversion.hh"
60#include "G4LowEnergyPhotoElectric.hh"
61#include "G4LowEnergyRayleigh.hh"
62#include "G4hLowEnergyIonisation.hh"
63
64#include "G4VeEnergyLoss.hh"
65#include "G4eIonisation.hh"
66#include "G4eBremsstrahlung.hh"
67#include "G4ComptonScattering.hh"
68#include "G4GammaConversion.hh"
69#include "G4PhotoElectricEffect.hh"
70
71#include "G4eplusAnnihilation.hh"
72
73#include "G4MuIonisation.hh"
74#include "G4MuBremsstrahlung.hh"
75#include "G4MuPairProduction.hh"
76
77#include "G4hIonisation.hh"
78#include "G4MultipleScattering.hh"
79#include "G4EnergyLossTables.hh"
80#include "G4ParticleChange.hh"
81#include "G4ParticleChange.hh"
82#include "G4DynamicParticle.hh"
83#include "G4AntiProton.hh"
84#include "G4Proton.hh"
85#include "G4Electron.hh"
86#include "G4Positron.hh"
87#include "G4Gamma.hh"
88#include "G4ForceCondition.hh"
89#include "G4Box.hh"
90#include "G4PVPlacement.hh"
91#include "G4Step.hh"
92#include "G4GRSVolume.hh"
93
94// New Histogramming (from AIDA and Anaphe):
95#include <memory> // for the auto_ptr(T>
96
97
98#include "AIDA/AIDA.h"
99
100#include "hTest/include/G4IonC12.hh"
101#include "hTest/include/G4IonAr40.hh"
102
103#include "G4Timer.hh"
104
105int main(int argc,char** argv)
106{
107 // HepTupleManager* hbookManager;
108
109 // -------------------------------------------------------------------
110 // Setup
111
112 G4int nPart =-1;
113 G4String nameMat = "Si";
114 G4int nProcess = 3;
115 G4bool usepaw = false;
116 G4bool fluct = false;
117 G4bool lowE = true;
118 G4int verbose = 0;
119 G4double emin = 0.01*MeV;
120 G4double emax = 100.0*MeV;
121 G4int nstatf = 10;
122 G4double xstatf = 1.0/(G4double)nstatf;
123 G4int nbin = 1000;
124 G4String hFile = "hbook.paw";
125 G4double theStep = 0.01*micrometer;
126 G4double range = 1.0*micrometer;
127 G4double cutG = 10.0*mm;
128 G4double cutE = 10.0*mm;
129 G4Material* material = 0;
130 G4String name[3] = {"Ionisation", "Bremsstrahlung",
131 "Ionisation+Bremsstrahlung"};
132 G4bool setBarkasOff = false;
133 G4bool setNuclearOff= true;
134
135 G4cout.setf( ios::scientific, ios::floatfield );
136
137
138 // -------------------------------------------------------------------
139 // Control on input
140
141 if(argc < 2) {
142 G4cout << "Input file is not specified! Exit" << G4endl;
143 exit(1);
144 }
145
146 ifstream* fin = new ifstream();
147 string fname = argv[1];
148 fin->open(fname.c_str());
149 if( !fin->is_open()) {
150 G4cout << "Input file <" << fname << "> does not exist! Exit" << G4endl;
151 exit(1);
152 }
153
154 // -------------------------------------------------------------------
155 //--------- Materials definition ---------
156
157 G4Material* ma[16];
158 ma[0] = new G4Material("Be", 4., 9.01*g/mole, 1.848*g/cm3);
159 ma[1] = new G4Material("Graphite",6., 12.00*g/mole, 2.265*g/cm3 );
160 ma[1]->SetChemicalFormula("Graphite");
161 ma[2] = new G4Material("Al", 13., 26.98*g/mole, 2.7 *g/cm3);
162 ma[3] = new G4Material("Si", 14., 28.055*g/mole, 2.33*g/cm3);
163
164 ma[4] = new G4Material("LAr", 18., 39.95*g/mole, 1.393*g/cm3);
165 ma[5] = new G4Material("Fe", 26., 55.85*g/mole, 7.87*g/cm3);
166 ma[6] = new G4Material("Cu", 29., 63.55*g/mole, 8.96*g/cm3);
167 ma[7] = new G4Material("W", 74., 183.85*g/mole, 19.30*g/cm3);
168 ma[8] = new G4Material("Pb",82., 207.19*g/mole, 11.35*g/cm3);
169 ma[9] = new G4Material("U", 92., 238.03*g/mole, 18.95*g/cm3);
170
171 G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
172 G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
173 G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
174 G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
175 G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole);
176
177 ma[10] = new G4Material("O2", 8., 16.00*g/mole, 1.1*g/cm3);
178 ma[10]->SetChemicalFormula("O_2");
179
180 ma[11] = new G4Material ("Water" , 1.*g/cm3, 2);
181 ma[11]->AddElement(H,2);
182 ma[11]->AddElement(O,1);
183 ma[11]->SetChemicalFormula("H_2O");
184
185 ma[12] = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
186 ma[12]->AddElement(H,6);
187 ma[12]->AddElement(C,2);
188 ma[12]->SetChemicalFormula("C_2H_6");
189
190 ma[13] = new G4Material ("CsI" , 4.53*g/cm3, 2);
191 ma[13]->AddElement(Cs,1);
192 ma[13]->AddElement(I,1);
193 ma[13]->SetChemicalFormula("CsI");
194
195 ma[14] = new G4Material("H2", 1., 1.00794*g/mole, 1.*g/cm3);
196 ma[14]->SetChemicalFormula("H_2");
197 ma[15] = new G4Material("Au", 79., 196.96655*g/mole, 19.32*g/cm3);
198
199
200 static const G4MaterialTable* theMaterialTable =
201 G4Material::GetMaterialTable();
202
203 G4int nMaterials = G4Material::GetNumberOfMaterials();
204 G4cout << "Available materials are: " << G4endl;
205 G4int mat;
206 for (mat = 0; mat < nMaterials; mat++) {
207 G4cout << mat << ") " << ma[mat]->GetName() << G4endl;
208 }
209
210 G4cout << "Available processes are: " << G4endl;
211 for (mat = 0; mat < 2; mat++) {
212 G4cout << mat << ") " << name[mat] << G4endl;
213 }
214
215 // Particle definitions
216
217 G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
218 G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
219 G4ParticleDefinition* positron = G4Positron::PositronDefinition();
220 G4ParticleDefinition* proton = G4Proton::ProtonDefinition();
221 G4ParticleDefinition* antiproton = G4AntiProton::AntiProtonDefinition();
222 G4ParticleDefinition* c12 = G4IonC12::IonC12Definition();
223 G4ParticleDefinition* ar40 = G4IonAr40::IonAr40Definition();
224 G4ParticleDefinition* part = electron;
225
226 G4hLowEnergyIonisation* hionle = 0;
227 G4hLowEnergyIonisation* ionle = 0;
228 G4LowEnergyIonisation* eionle = 0;
229 G4LowEnergyBremsstrahlung* ebrle = 0;
230 G4hIonisation* ionst = 0;
231 G4hIonisation* hionst = 0;
232 G4eIonisation* eionst = 0;
233 G4eBremsstrahlung* ebrst = 0;
234
235 G4cout << "Process is initialized" << G4endl;;
236
237 // Geometry
238
239 G4double initX = 0.;
240 G4double initY = 0.;
241 G4double initZ = 1.;
242 G4double dimX = 100.0*cm;
243 G4double dimY = 100.0*cm;
244 G4double dimZ = 100.0*cm;
245
246 G4Box* sFrame = new G4Box ("Box",dimX, dimY, dimZ);
247 G4LogicalVolume* lFrame = new G4LogicalVolume(sFrame,material,"Box",0,0,0);
248 G4PVPlacement* pFrame = new G4PVPlacement(0,G4ThreeVector(),"Box",
249 lFrame,0,false,0);
250
251 // -------------------------------------------------------------------
252 // ---- Read input file
253 G4cout << "Available commands are: " << G4endl;
254 G4cout << "#events" << G4endl;
255 G4cout << "#particle" << G4endl;
256 G4cout << "#emin(MeV)" << G4endl;
257 G4cout << "#emax(MeV)" << G4endl;
258 G4cout << "#nbin" << G4endl;
259 G4cout << "#cutG(mm)" << G4endl;
260 G4cout << "#cutE(mm)" << G4endl;
261 G4cout << "#range(mm)" << G4endl;
262 G4cout << "#step(mm)" << G4endl;
263 G4cout << "#material" << G4endl;
264 G4cout << "#process" << G4endl;
265 G4cout << "#domain" << G4endl;
266 G4cout << "#paw" << G4endl;
267 G4cout << "#verbose" << G4endl;
268 G4cout << "#run" << G4endl;
269 G4cout << "#exit" << G4endl;
270 G4cout << "#barkas" << G4endl;
271 G4cout << "#nuclear" << G4endl;
272 G4cout << pFrame << G4endl;
273
274 G4ProcessManager *gmanager, *elecManager, *positManager,
275 *protManager, *aprotManager, *ionManager;
276
277 string line, line1;
278 G4bool end = true;
279 for(G4int run=0; run<100; run++) {
280 do {
281 (*fin) >> line;
282 G4cout << "Next line " << line << G4endl;
283 if(line == "#particle") {
284 (*fin) >> nPart;
285 } else if(line == "#emin(MeV)") {
286 (*fin) >> emin;
287 emin *= MeV;
288 } else if(line == "#emax(MeV)") {
289 (*fin) >> emax;
290 emax *= MeV;
291 } else if(line == "#nbin") {
292 (*fin) >> nbin;
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 line1 = "";
311 (*fin) >> line1;
312 if(line1 == "lowenergy") {lowE = true;}
313 else {lowE = false;}
314 } else if(line == "#paw") {
315 usepaw = true;
316 (*fin) >> hFile;
317 } else if(line == "#run") {
318 break;
319 } else if(line == "#verbose") {
320 (*fin) >> verbose;
321 } else if(line == "#fluct") {
322 fluct = true;
323 } else if(line == "#exit") {
324 end = false;
325 break;
326 } else if(line == "#barkas") {
327 line1 = "";
328 (*fin) >> line1;
329 G4cout << line1 << G4endl;
330 if(line1 == "off") setBarkasOff = true;
331 if(line1 == "on") setBarkasOff = false;
332 } else if(line == "#nuclear") {
333 line1 = "";
334 (*fin) >> line1;
335 G4cout << line1 << G4endl;
336 if(line1 == "off") setNuclearOff = true;
337 if(line1 == "on") setNuclearOff = false;
338 }
339 } while(end);
340
341 if(!end) break;
342
343 G4cout << "###### Start new run # " << run << " #####" << G4endl;
344 if(nPart == 0) {
345 part = gamma;
346 } else if(nPart == 1) {
347 part = electron;
348 } else if(nPart == 2) {
349 part = positron;
350 } else if(nPart == 3) {
351 part = proton;
352 } else if(nPart == 4) {
353 part = antiproton;
354 } else if(nPart == 5) {
355 part = c12;
356 } else if(nPart == 6) {
357 part = ar40;
358 } else {
359 G4cout << "Particle #" << nPart
360 << " is absent in the list of particles: Exit" << G4endl;
361 end = false;
362 break;
363 }
364 if(nProcess < 0 || nProcess > 2) {
365 G4cout << "Process #" << nProcess
366 << " is absent in the list of processes: Exit" << G4endl;
367 end = false;
368 break;
369 }
370
371
372 for (mat = 0; mat < nMaterials; mat++) {
373 material = ma[mat];
374 if(nameMat == material->GetName()) break;
375 }
376
377 G4cout << "The particle: " << part->GetParticleName() << G4endl;
378 G4cout << "The material: " << material->GetName() << G4endl;
379 G4cout << "The cut on e-:" << cutE/mm << " mm" << G4endl;
380 G4cout << "The cut on g: " << cutG/mm << " mm" << G4endl;
381 G4cout << "The step: " << theStep/mm << " mm" << G4endl;
382
383 // -------------------------------------------------------------------
384 // ---- HBOOK initialization
385
386 G4double emin10 = std::log10(emin/MeV);
387 G4double emax10 = std::log10(emax/MeV);
388 G4double bin = (emax10 - emin10) / (G4double)(nbin-1);
389
390 // Creating the analysis factory
391 std::auto_ptr< AIDA::IAnalysisFactory > af( AIDA_createAnalysisFactory() );
392
393 // Creating the tree factory
394 std::auto_ptr< AIDA::ITreeFactory > tf( af->createTreeFactory() );
395
396 // Creating a tree mapped to a new hbook file.
397 std::auto_ptr< AIDA::ITree > tree( tf->create( hFile,"hbook",false,false ) );
398 std::cout << "Tree store : " << tree->storeName() << std::endl;
399
400 // Creating a tuple factory, whose tuples will be handled by the tree
401 std::auto_ptr< AIDA::ITupleFactory > tpf( af->createTupleFactory( *tree ) );
402
403 AIDA::IHistogram1D* hist[4];
404 AIDA::ITuple* ntuple1 = 0;
405
406 if(usepaw) {
407
408 // ---- primary ntuple ------
409 // If using Anaphe HBOOK implementation, there is a limitation on the length of the
410 // variable names in a ntuple
411 ntuple1 = tpf->create( "100", "tuple", "float ekin, dedx" );
412
413
414 // Creating a histogram factory, whose histograms will be handled by the tree
415 std::auto_ptr< AIDA::IHistogramFactory > hf( af->createHistogramFactory( *tree ) );
416
417 // Creating an 1-dimensional histogram in the root directory of the tree
418
419 hist[0] = hf->createHistogram1D("11","Stopping power (MeV*cm**2/g)",
420 nbin,emin10,emax10);
421 hist[1] = hf->createHistogram1D("12","Stopping power (MeV/mm)",
422 nbin,emin10,emax10);
423 hist[2] = hf->createHistogram1D("13","Step limit (mm)",
424 nbin,emin10,emax10);
425 hist[3] = hf->createHistogram1D("14","Number of secondaries",
426 nbin,emin10,emax10);
427
428
429 G4cout<< "Histograms is initialised" << G4endl;
430 }
431
432 gamma->SetCuts(cutG);
433 electron->SetCuts(cutE);
434 // positron->SetCuts(cutE);
435
436 // Processes - all new
437 G4bool success = false;
438
439 gmanager = new G4ProcessManager(gamma);
440 gmanager->AddDiscreteProcess(new G4LowEnergyPhotoElectric());
441 gmanager->AddDiscreteProcess(new G4LowEnergyCompton());
442 gmanager->AddDiscreteProcess(new G4LowEnergyGammaConversion());
443
444 if(lowE) {
445
446 elecManager = new G4ProcessManager(electron);
447 electron->SetProcessManager(elecManager);
448 eionle = new G4LowEnergyIonisation();
449 if(!fluct) eionle->SetEnlossFluc(false);
450 ebrle = new G4LowEnergyBremsstrahlung();
451 if(!fluct) ebrle->SetEnlossFluc(false);
452 elecManager->AddProcess(eionle);
453 elecManager->AddProcess(ebrle);
454 eionle->BuildPhysicsTable(*electron);
455 ebrle->BuildPhysicsTable(*electron);
456 if(nPart == 1) {
457 if(nProcess == 2) {
458 success = true;
459 }
460
461 } else if (nPart == 3) {
462 protManager = new G4ProcessManager(proton);
463 proton->SetProcessManager(protManager);
464 hionle = new G4hLowEnergyIonisation();
465 if(!fluct) hionle->SetEnlossFluc(false);
466 protManager->AddProcess(hionle);
467 if(setNuclearOff) hionle->SetNuclearStoppingOff();
468 if(!setNuclearOff) hionle->SetNuclearStoppingOn();
469 if(setBarkasOff) hionle->SetBarkasOff();
470 if(!setBarkasOff) hionle->SetBarkasOn();
471 hionle->SetVerboseLevel(verbose);
472 hionle->BuildPhysicsTable(*proton);
473 success = true;
474
475 } else if (nPart == 4) {
476 aprotManager = new G4ProcessManager(antiproton);
477 antiproton->SetProcessManager(aprotManager);
478 hionle = new G4hLowEnergyIonisation();
479 if(!fluct) hionle->SetEnlossFluc(false);
480 aprotManager->AddProcess(hionle);
481 if(setNuclearOff) hionle->SetNuclearStoppingOff();
482 if(!setNuclearOff) hionle->SetNuclearStoppingOn();
483 if(setBarkasOff) hionle->SetBarkasOff();
484 if(!setBarkasOff) hionle->SetBarkasOn();
485 hionle->SetVerboseLevel(verbose);
486 hionle->BuildPhysicsTable(*antiproton);
487 success = true;
488
489 } else if (nPart == 5 || nPart == 6) {
490 protManager = new G4ProcessManager(proton);
491 proton->SetProcessManager(protManager);
492 hionle = new G4hLowEnergyIonisation();
493 if(!fluct) hionle->SetEnlossFluc(false);
494 protManager->AddProcess(hionle);
495 if(setNuclearOff) hionle->SetNuclearStoppingOff();
496 if(!setNuclearOff) hionle->SetNuclearStoppingOn();
497 if(setBarkasOff) hionle->SetBarkasOff();
498 if(!setBarkasOff) hionle->SetBarkasOn();
499 hionle->SetVerboseLevel(verbose);
500 hionle->BuildPhysicsTable(*proton);
501 ionManager = new G4ProcessManager(part);
502 part->SetProcessManager(ionManager);
503 ionle = new G4hLowEnergyIonisation();
504 if(!fluct) ionle->SetEnlossFluc(false);
505 if(setNuclearOff) ionle->SetNuclearStoppingOff();
506 if(!setNuclearOff) ionle->SetNuclearStoppingOn();
507 if(setBarkasOff) ionle->SetBarkasOff();
508 if(!setBarkasOff) ionle->SetBarkasOn();
509 ionManager->AddProcess(ionle);
510 ionle->SetVerboseLevel(verbose);
511 ionle->BuildPhysicsTable(*part);
512 success = true;
513 }
514
515 } else {
516
517 elecManager = new G4ProcessManager(electron);
518 electron->SetProcessManager(elecManager);
519 eionst = new G4eIonisation();
520 if(!fluct) eionst->SetEnlossFluc(false);
521 elecManager->AddProcess(eionst);
522 ebrst = new G4eBremsstrahlung();
523 if(!fluct) ebrst->SetEnlossFluc(false);
524 elecManager->AddProcess(ebrst);
525 eionst->BuildPhysicsTable(*electron);
526 ebrst->BuildPhysicsTable(*electron);
527 if(nPart == 1) {
528 if(nProcess == 2) {
529 success = true;
530 }
531 } else if(nPart == 2) {
532 positManager = new G4ProcessManager(positron);
533 positron->SetProcessManager(positManager);
534 if(nProcess == 0) {
535 eionst = new G4eIonisation();
536 if(!fluct) eionst->SetEnlossFluc(false);
537 positManager->AddProcess(eionst);
538 eionst->BuildPhysicsTable(*positron);
539 success = true;
540 } else if(nProcess == 1) {
541 ebrst = new G4eBremsstrahlung();
542 if(!fluct) ebrst->SetEnlossFluc(false);
543 positManager->AddProcess(ebrst);
544 ebrst->BuildPhysicsTable(*positron);
545 success = true;
546 } else if(nProcess == 2) {
547 eionst = new G4eIonisation();
548 ebrst = new G4eBremsstrahlung();
549 if(!fluct) eionst->SetEnlossFluc(false);
550 if(!fluct) ebrst->SetEnlossFluc(false);
551 positManager->AddProcess(eionst);
552 positManager->AddProcess(ebrst);
553 eionst->BuildPhysicsTable(*positron);
554 ebrst->BuildPhysicsTable(*positron);
555 success = true;
556 }
557
558 } else if (nPart == 3) {
559 protManager = new G4ProcessManager(proton);
560 proton->SetProcessManager(protManager);
561 hionst = new G4hIonisation();
562 if(!fluct) hionst->SetEnlossFluc(false);
563 protManager->AddProcess(hionst);
564 // hionst->SetVerboseLevel(verbose);
565 hionst->BuildPhysicsTable(*proton);
566 success = true;
567
568 } else if (nPart == 4) {
569 aprotManager = new G4ProcessManager(antiproton);
570 antiproton->SetProcessManager(aprotManager);
571 hionst = new G4hIonisation();
572 if(!fluct) hionst->SetEnlossFluc(false);
573 aprotManager->AddProcess(hionst);
574 hionst->SetVerboseLevel(verbose);
575 hionst->BuildPhysicsTable(*antiproton);
576 success = true;
577
578 } else if (nPart == 5 || nPart == 6) {
579 protManager = new G4ProcessManager(proton);
580 proton->SetProcessManager(protManager);
581 hionst = new G4hIonisation();
582 if(!fluct) hionst->SetEnlossFluc(false);
583 protManager->AddProcess(hionst);
584 // hionst->SetVerboseLevel(verbose);
585 hionst->BuildPhysicsTable(*proton);
586 ionManager = new G4ProcessManager(part);
587 part->SetProcessManager(ionManager);
588 ionst = new G4hIonisation();
589 if(!fluct) ionst->SetEnlossFluc(false);
590 ionManager->AddProcess(ionst);
591 ionst->SetVerboseLevel(verbose);
592 ionst->BuildPhysicsTable(*part);
593 success = true;
594 }
595 }
596
597 if(success) G4cout << "Physics tables are built" << G4endl;
598 else G4cout << "Physics tables are not built!!!" << G4endl;
599
600 G4cout << "gCut(MeV)= " << gamma->GetEnergyThreshold(material)/MeV << G4endl;
601 G4cout << "eCut(MeV)= " << electron->GetEnergyThreshold(material)/MeV << G4endl;
602
603
604 // Create a DynamicParticle
605
606 G4ParticleMomentum gDir(initX,initY,initZ);
607 G4double gEnergy = emax;
608 G4DynamicParticle dParticle(part,gDir,gEnergy);
609
610 // Track
611 G4ThreeVector aPosition(0.,0.,0.);
612 G4double aTime = 0. ;
613
614 G4Track* gTrack;
615 gTrack = new G4Track(&dParticle,aTime,aPosition);
616
617 // Step
618
619 G4Step* step;
620 step = new G4Step();
621 step->SetTrack(gTrack);
622
623 G4StepPoint *aPoint, *bPoint;
624 aPoint = new G4StepPoint();
625 aPoint->SetPosition(aPosition);
626 aPoint->SetMaterial(material);
627 G4double safety = 10000.*cm;
628 aPoint->SetSafety(safety);
629 step->SetPreStepPoint(aPoint);
630
631 bPoint = aPoint;
632 G4ThreeVector bPosition(0.,0.,theStep);
633 bPoint->SetPosition(bPosition);
634 step->SetPostStepPoint(bPoint);
635 step->SetStepLength(theStep);
636
637 if(!fluct) {
638 nstatf = 1;
639 xstatf = 1.0;
640 }
641
642 G4Timer* timer = new G4Timer();
643 timer->Start();
644
645 G4double le = emin10 - bin;
646
647 for (G4int iter=0; iter<nbin; iter++) {
648
649 le += bin;
650 G4double e = std::pow(10.0,le) * MeV;
651 gTrack->SetStep(step);
652 gTrack->SetKineticEnergy(e);
653
654 for (G4int jj=0; jj<nstatf; jj++) {
655
656 G4double x = 0.0;
657 G4VParticleChange* aChange = 0;
658
659 if(lowE) {
660
661 if(nPart == 1) {
662 if(nProcess == 0) {
663 x = eionle->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
664 aChange = eionle->AlongStepDoIt(*gTrack,*step);
665
666 } else if(nProcess == 1) {
667 x = ebrle->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
668 aChange = ebrle->AlongStepDoIt(*gTrack,*step);
669 } else if(nProcess == 2) {
670 x = eionle->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
671 aChange = eionle->AlongStepDoIt(*gTrack,*step);
672 }
673
674 } else if (nPart == 3 ) {
675 x = hionle->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
676 aChange = hionle->AlongStepDoIt(*gTrack,*step);
677
678 } else if (nPart == 4) {
679 x = hionle->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
680 aChange = hionle->AlongStepDoIt(*gTrack,*step);
681
682 } else if (nPart == 5 || nPart == 6) {
683 x = ionle->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
684 aChange = ionle->AlongStepDoIt(*gTrack,*step);
685
686 }
687
688 } else {
689
690
691 if(nPart == 1) {
692 if(nProcess == 0) {
693 x = eionst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
694 aChange = eionst->AlongStepDoIt(*gTrack,*step);
695
696 } else if(nProcess == 1) {
697 x = ebrst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
698 aChange = ebrst->AlongStepDoIt(*gTrack,*step);
699 } else if(nProcess == 2) {
700 x = eionst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
701 aChange = eionst->AlongStepDoIt(*gTrack,*step);
702 }
703
704 } else if(nPart == 2) {
705 if(nProcess == 0) {
706 x = eionst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
707 aChange = eionst->AlongStepDoIt(*gTrack,*step);
708
709 } else if(nProcess == 1) {
710 x = ebrst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
711 aChange = ebrst->AlongStepDoIt(*gTrack,*step);
712 } else if(nProcess == 2) {
713 x = eionst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
714 aChange = eionst->AlongStepDoIt(*gTrack,*step);
715 }
716
717 } else if (nPart == 3) {
718 x = hionst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
719 aChange = hionst->AlongStepDoIt(*gTrack,*step);
720
721 } else if (nPart == 4) {
722 x = hionst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
723 aChange = hionst->AlongStepDoIt(*gTrack,*step);
724
725 } else if (nPart == 5 || nPart == 6) {
726 x = ionst->GetContinuousStepLimit(*gTrack,theStep,theStep,safety);
727 aChange = ionst->AlongStepDoIt(*gTrack,*step);
728 }
729 }
730
731 G4double delx = theStep;
732 G4double de = aChange->GetLocalEnergyDeposit();
733 G4int n = aChange->GetNumberOfSecondaries();
734
735
736 //G4cout << " de(MeV) = " << de/MeV << " n= " << n << G4endl;
737
738 if(n > 0) {
739 for(G4int i=0; i<n; i++) {
740 de += (aChange->GetSecondary(i))->GetKineticEnergy();
741 if(verbose) {
742 G4cout << "add "
743 << ((aChange->GetSecondary(i))->GetKineticEnergy())/eV
744 << " eV" << G4endl;
745 }
746 }
747 }
748 G4double st = de/(delx*(material->GetDensity()));
749 st *= gram/(cm*cm*MeV);
750 G4double s = de*mm/(delx*MeV);
751
752 G4cout << "E(MeV)= " << e/MeV << " dE/dx(MeV/mm)= " << s << G4endl;
753
754 if(verbose) {
755 G4cout << "Iteration = " << iter
756 << " E = " << e/MeV << " MeV; StepLimit= "
757 << x/mm << " mm; de= "
758 << de/eV << " eV; dE/dx= "
759 << st << " MeV*cm^2/g" << G4endl;
760 }
761
762 if(x > 1000.0*meter) x = 1000.0*meter;
763
764 if (usepaw) {
765 float st10 = -5.0;
766 if(st > 1.e-5) st10 = (float)log10(st);
767 if(verbose>1) {
768 G4cout << " de(MeV) = " << de/MeV
769 << G4endl;
770 G4cout << " n1= " << ntuple1->findColumn("ekin")
771 << " n2= " << ntuple1->findColumn("dedx")
772 << G4endl;
773 }
774 ntuple1->fill( ntuple1->findColumn("ekin"), (float)le);
775 ntuple1->fill( ntuple1->findColumn("dedx"), st10);
776 ntuple1->addRow();
777 // G4cout << "ntuple is filled " << G4endl;
778
779 hist[0]->fill(le,st*xstatf);
780 hist[1]->fill(le,s*xstatf);
781 hist[2]->fill(le,x*xstatf/mm);
782 hist[3]->fill(le,xstatf*(G4double)n);
783 }
784 }
785 }
786
787 timer->Stop();
788 G4cout << " " << *timer << G4endl;
789 delete timer;
790
791 // Committing the transaction with the tree
792 if(usepaw) {
793 std::cout << "Committing..." << std::endl;
794 tree->commit();
795 std::cout << "Closing the tree..." << std::endl;
796 tree->close();
797 }
798 G4cout << "###### End of run # " << run << " ######" << G4endl;
799
800
801 } while(end);
802 G4cout << "###### End of test #####" << G4endl;
803}
804
805#include "hTest/src/G4IonC12.cc"
806#include "hTest/src/G4IonAr40.cc"
Note: See TracBrowser for help on using the repository browser.