source: trunk/source/processes/electromagnetic/lowenergy/test/G4MeanFreePathTest.cc@ 1201

Last change on this file since 1201 was 1199, checked in by garnier, 16 years ago

nvx fichiers dans CVS
File size: 18.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: G4MeanFreePathTest
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 on base of Complex test
39//
40// Creation date: 17 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
55#include "G4LowEnergyIonisation.hh"
56#include "G4LowEnergyBremsstrahlung.hh"
57#include "G4LowEnergyCompton.hh"
58#include "G4LowEnergyGammaConversion.hh"
59#include "G4LowEnergyPhotoElectric.hh"
60#include "G4LowEnergyRayleigh.hh"
61#include "G4hLowEnergyIonisation.hh"
62
63#include "G4eIonisation.hh"
64#include "G4eBremsstrahlung.hh"
65#include "G4ComptonScattering.hh"
66#include "G4GammaConversion.hh"
67#include "G4PhotoElectricEffect.hh"
68
69#include "G4eplusAnnihilation.hh"
70
71#include "G4MuIonisation.hh"
72#include "G4MuBremsstrahlung.hh"
73#include "G4MuPairProduction.hh"
74
75#include "G4hIonisation.hh"
76
77#include "G4MultipleScattering.hh"
78
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
92#include "G4Step.hh"
93#include "G4GRSVolume.hh"
94
95// New Histogramming (from AIDA and Anaphe):
96#include <memory> // for the auto_ptr(T>
97
98#include "AIDA/IAnalysisFactory.h"
99
100#include "AIDA/ITreeFactory.h"
101#include "AIDA/ITree.h"
102
103#include "AIDA/IHistogramFactory.h"
104#include "AIDA/IHistogram1D.h"
105#include "AIDA/IHistogram2D.h"
106//#include "AIDA/IHistogram3D.h"
107
108#include "AIDA/ITupleFactory.h"
109#include "AIDA/ITuple.h"
110
111#include "G4UnitsTable.hh"
112
113int main(int argc,char** argv)
114{
115
116 // -------------------------------------------------------------------
117 // Setup
118
119 G4int nEvt = 100;
120 G4int nPart =-1;
121 G4String nameMat = "Si";
122 G4int nProcess = 3;
123 G4bool usepaw = false;
124 G4bool lowE = true;
125 G4int verbose = 0;
126 G4double emin = 0.0*MeV;
127 G4double emax = 100.0*MeV;
128 G4int nbin = 1000;
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* ma[15];
161 /*
162 ma[0] = new G4Material("Be", 4., 9.01*g/mole, 1.848*g/cm3);
163 ma[1] = new G4Material("Graphite",6., 12.00*g/mole, 2.265*g/cm3 );
164 ma[2] = new G4Material("Al", 13., 26.98*g/mole, 2.7 *g/cm3);
165 ma[3] = new G4Material("Si", 14., 28.055*g/mole, 2.33*g/cm3);
166 ma[4] = new G4Material("LAr", 18., 39.95*g/mole, 1.393*g/cm3);
167 ma[5] = new G4Material("Fe", 26., 55.85*g/mole, 7.87*g/cm3);
168 ma[6] = new G4Material("Cu", 29., 63.55*g/mole, 8.96*g/cm3);
169 ma[7] = new G4Material("W", 74., 183.85*g/mole, 19.30*g/cm3);
170 ma[8] = new G4Material("Pb", 82., 207.19*g/mole, 11.35*g/cm3);
171 ma[9] = new G4Material("U", 92., 238.03*g/mole, 18.95*g/cm3);
172 */
173 G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
174 G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
175 G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
176 G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
177 G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole);
178
179 //ma[10] = new G4Material("O2", 8., 16.00*g/mole, 1.1*g/cm3);
180 ma[11] = new G4Material ("Water" , 1.*g/cm3, 2);
181 ma[11]->AddElement(H,2);
182 ma[11]->AddElement(O,1);
183
184 ma[12] = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
185 ma[12]->AddElement(H,6);
186 ma[12]->AddElement(C,2);
187
188 ma[13] = new G4Material ("CsI" , 4.53*g/cm3, 2);
189 ma[13]->AddElement(Cs,1);
190 ma[13]->AddElement(I,1);
191
192 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
193
194 G4int nMaterials = G4Material::GetNumberOfMaterials();
195 G4cout << "Available materials are: " << G4endl;
196 G4int mat;
197 for (mat = 0; mat < nMaterials; mat++) {
198 G4cout << mat << ") " << (*theMaterialTable)[mat]->GetName() << G4endl;
199 }
200
201 G4cout << "Available processes are: " << G4endl;
202 for (mat = 0; mat < 6; mat++) {
203 G4cout << mat << ") " << name[mat] << G4endl;
204 }
205
206 // Particle definitions
207
208 G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
209 G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
210 G4ParticleDefinition* positron = G4Positron::PositronDefinition();
211 G4ParticleDefinition* proton = G4Proton::ProtonDefinition();
212 G4ParticleDefinition* antiproton = G4AntiProton::AntiProtonDefinition();
213 G4ParticleDefinition* part = gamma;
214
215 G4LowEnergyIonisation* eionle = 0;
216 G4LowEnergyBremsstrahlung* ebrle = 0;
217 G4LowEnergyCompton* comple = 0;
218 G4LowEnergyGammaConversion* convle = 0;
219 G4LowEnergyPhotoElectric* pele = 0;
220 G4LowEnergyRayleigh* rayle = 0;
221 G4eIonisation* eionst = 0;
222 G4eBremsstrahlung* ebrst = 0;
223 G4ComptonScattering* compst = 0;
224 G4GammaConversion* convst = 0;
225 G4PhotoElectricEffect* pest = 0;
226
227 G4cout << "Process is initialized" << G4endl;;
228
229 // Geometry
230
231 G4double initX = 0.;
232 G4double initY = 0.;
233 G4double initZ = 1.;
234 G4double dimX = 100.0*cm;
235 G4double dimY = 100.0*cm;
236 G4double dimZ = 100.0*cm;
237
238 G4Box* sFrame = new G4Box ("Box",dimX, dimY, dimZ);
239 G4LogicalVolume* lFrame = new G4LogicalVolume(sFrame,material,"Box", 0, 0, 0);
240 G4PVPlacement* pFrame = new G4PVPlacement(0,G4ThreeVector(),"Box",lFrame,0,false,0);
241 if(pFrame) G4cout << "Geometry: " << pFrame->GetName() << G4endl;
242
243 // -------------------------------------------------------------------
244 // ---- Read input file
245 G4cout << "Available commands are: " << G4endl;
246 G4cout << "#events" << G4endl;
247 G4cout << "#particle" << G4endl;
248 G4cout << "#emin(MeV)" << G4endl;
249 G4cout << "#emax(MeV)" << G4endl;
250 G4cout << "#nbin" << G4endl;
251 G4cout << "#cutG(mm)" << G4endl;
252 G4cout << "#cutE(mm)" << G4endl;
253 G4cout << "#range(mm)" << G4endl;
254 G4cout << "#step(mm)" << G4endl;
255 G4cout << "#material" << G4endl;
256 G4cout << "#process" << G4endl;
257 G4cout << "#domain" << G4endl;
258 G4cout << "#paw" << G4endl;
259 G4cout << "#verbose" << G4endl;
260 G4cout << "#run" << G4endl;
261 G4cout << "#exit" << G4endl;
262
263 G4ProcessManager *elecManager, *positManager, *gammaManager, *protManager, *aprotManager;
264
265 elecManager = new G4ProcessManager(electron);
266 electron->SetProcessManager(elecManager);
267
268 positManager = new G4ProcessManager(positron);
269 positron->SetProcessManager(positManager);
270
271 gammaManager = new G4ProcessManager(gamma);
272 gamma->SetProcessManager(gammaManager);
273
274 protManager = new G4ProcessManager(proton);
275 proton->SetProcessManager(protManager);
276
277 aprotManager = new G4ProcessManager(antiproton);
278 antiproton->SetProcessManager(aprotManager);
279
280 // G4bool ionis = true;
281
282 string line, line1;
283 G4bool end = true;
284 for(G4int run=0; run<100; run++) {
285 do {
286 (*fin) >> line;
287 G4cout << "Next line " << line << G4endl;
288 if(line == "#events") {
289 (*fin) >> nEvt;
290 if(nEvt < 1) nEvt = 1;
291 } else if(line == "#particle") {
292 (*fin) >> nPart;
293 } else if(line == "#emin(MeV)") {
294 (*fin) >> emin;
295 emin *= MeV;
296 } else if(line == "#emax(MeV)") {
297 (*fin) >> emax;
298 emax *= MeV;
299 } else if(line == "#nbin") {
300 (*fin) >> nbin;
301 } else if(line == "#cutG(mm)") {
302 (*fin) >> cutG;
303 cutG *= mm;
304 } else if(line == "#cutE(mm)") {
305 (*fin) >> cutE;
306 cutE *= mm;
307 } else if(line == "#range(mm)") {
308 (*fin) >> range;
309 range *= mm;
310 } else if(line == "#step(mm)") {
311 (*fin) >> theStep;
312 theStep *= mm;
313 } else if(line == "#material") {
314 (*fin) >> nameMat;
315 } else if(line == "#process") {
316 (*fin) >> nProcess;
317 } else if(line == "#domain") {
318 (*fin) >> line1;
319 if(line1 == "lowenergy") {lowE = true;}
320 else {lowE = false;}
321 } else if(line == "#paw") {
322 usepaw = true;
323 (*fin) >> hFile;
324 } else if(line == "#run") {
325 break;
326 } else if(line == "#verbose") {
327 (*fin) >> verbose;
328 } else if(line == "#exit") {
329 end = false;
330 break;
331 }
332 } while(end);
333
334 if(!end) break;
335
336 G4cout << "###### Start new run # " << run << " #####" << G4endl;
337 if(nPart == 0) {
338 part = gamma;
339 } else if(nPart == 1) {
340 part = electron;
341 } else if(nPart == 2) {
342 part = positron;
343 } else if(nPart == 3) {
344 part = proton;
345 } else if(nPart == 4) {
346 part = antiproton;
347 } else {
348 G4cout << "Particle #" << nPart
349 << " is absent in the list of particles: Exit" << G4endl;
350 break;
351 }
352 if(nProcess < 0 || nProcess > 5) {
353 G4cout << "Process #" << nProcess
354 << " is absent in the list of processes: Exit" << G4endl;
355 break;
356 }
357
358 for (mat = 0; mat < nMaterials; mat++) {
359 material = (*theMaterialTable)[mat];
360 if(nameMat == material->GetName()) break;
361 }
362
363 G4cout << "The particle: " << part->GetParticleName() << G4endl;
364 G4cout << "The material: " << material->GetName() << G4endl;
365 G4cout << "The cut on e-:" << cutE/mm << " mm" << G4endl;
366 G4cout << "The cut on g: " << cutG/mm << " mm" << G4endl;
367 G4cout << "The step: " << theStep/mm << " mm" << G4endl;
368
369 // -------------------------------------------------------------------
370 // ---- HBOOK initialization
371
372 // Creating the analysis factory
373 std::auto_ptr< IAnalysisFactory > af( AIDA_createAnalysisFactory() );
374
375 // Creating the tree factory
376 std::auto_ptr< ITreeFactory > tf( af->createTreeFactory() );
377
378 // Creating a tree mapped to a new hbook file.
379 std::auto_ptr< ITree > tree( tf->create( hFile, false,false, "hbook" ) );
380 std::cout << "Tree store : " << tree->storeName() << std::endl;
381
382 // Creating a tuple factory, whose tuples will be handled by the tree
383 // std::auto_ptr< ITupleFactory > tpf(af->createTupleFactory( *tree ));
384
385 IHistogram1D* hist[6];
386
387
388 // ---- Book a histogram and ntuples
389 G4cout << "Hbook file name: <" << hFile << ">" << G4endl;
390 G4double bin = (emax - emin) / (G4double)nbin;
391
392 // Creating a histogram factory, whose histograms will be
393 // handled by the tree
394
395 if(usepaw) {
396
397 std::auto_ptr< IHistogramFactory > hf(af->createHistogramFactory(*tree));
398
399 // Creating an 1-dimensional histogram in the root directory of the tree
400
401 hist[0] = hf->create1D("1","MeanFreePath (mm)",
402 nbin,emin/MeV,emax/MeV);
403 /*
404 hist[1] = hf->create1D("2","Bremsstrahlung (E in MeV)",
405 nbin,emin/MeV,emax/MeV);
406 hist[2] = hf->create1D("3","Compton (E in MeV)",
407 nbin,emin/MeV,emax/MeV);
408 hist[3] = hf->create1D("4","GammaConversion (E in MeV)",
409 nbin,emin/MeV,emax/MeV);
410 hist[4] = hf->create1D("5","PhotoElectric (E in MeV)",
411 nbin,emin/MeV,emax/MeV);
412 hist[5] = hf->create1D("6","Raylaigh (E in MeV)",
413 nbin,emin/MeV,emax/MeV);
414 */
415 }
416
417 gamma->SetCuts(cutG);
418 electron->SetCuts(cutE);
419 positron->SetCuts(cutE);
420
421 G4cout << "Gamma cut in energy(keV) = "
422 << gamma->GetEnergyThreshold(material)/keV
423 << G4endl;
424 G4cout << "Electron cut in energy(keV) = "
425 << electron->GetEnergyThreshold(material)/keV
426 << G4endl;
427
428 // Processes
429
430 if(lowE && !eionle) {
431
432 eionle = new G4LowEnergyIonisation();
433 //eionle->SetLowEnergyLimit(0.2*MeV);
434 ebrle = new G4LowEnergyBremsstrahlung();
435 comple = new G4LowEnergyCompton();
436 convle = new G4LowEnergyGammaConversion();
437 pele = new G4LowEnergyPhotoElectric();
438 rayle = new G4LowEnergyRayleigh();
439
440 elecManager->AddProcess(eionle);
441 elecManager->AddProcess(ebrle);
442 gammaManager->AddProcess(comple);
443 gammaManager->AddProcess(convle);
444 gammaManager->AddProcess(pele);
445 gammaManager->AddProcess(rayle);
446
447 eionle->BuildPhysicsTable(*electron);
448 G4cout << "Electron ionisation are built" << G4endl;
449 ebrle->BuildPhysicsTable(*electron);
450 comple->BuildPhysicsTable(*gamma);
451 convle->BuildPhysicsTable(*gamma);
452 pele->BuildPhysicsTable(*gamma);
453 rayle->BuildPhysicsTable(*gamma);
454
455 } else if(!lowE && !eionst) {
456
457 eionst = new G4eIonisation();
458 ebrst = new G4eBremsstrahlung();
459 compst = new G4ComptonScattering();
460 convst = new G4GammaConversion();
461 pest = new G4PhotoElectricEffect();
462
463 elecManager->AddProcess(eionst);
464 elecManager->AddProcess(ebrst);
465 gammaManager->AddProcess(compst);
466 gammaManager->AddProcess(convst);
467 gammaManager->AddProcess(pest);
468
469 eionst->BuildPhysicsTable(*electron);
470 ebrst->BuildPhysicsTable(*electron);
471 compst->BuildPhysicsTable(*gamma);
472 convst->BuildPhysicsTable(*gamma);
473 pest->BuildPhysicsTable(*gamma);
474 }
475
476 G4cout << "Physics tables are built" << G4endl;
477
478 // Create a DynamicParticle
479
480 G4ParticleMomentum gDir(initX,initY,initZ);
481 G4double gEnergy = emax;
482 G4DynamicParticle dParticle(part,gDir,gEnergy);
483
484 // Track
485 G4ThreeVector aPosition(0.,0.,0.);
486 G4double aTime = 0. ;
487
488 G4Track* gTrack;
489 gTrack = new G4Track(&dParticle,aTime,aPosition);
490
491 // Step
492
493 G4Step* step;
494 step = new G4Step();
495 step->SetTrack(gTrack);
496
497 G4StepPoint *aPoint, *bPoint;
498 aPoint = new G4StepPoint();
499 aPoint->SetPosition(aPosition);
500 aPoint->SetMaterial(material);
501 G4double safety = 10000.*cm;
502 aPoint->SetSafety(safety);
503 step->SetPreStepPoint(aPoint);
504
505 bPoint = aPoint;
506 G4ThreeVector bPosition(0.,0.,theStep);
507 bPoint->SetPosition(bPosition);
508 step->SetPostStepPoint(bPoint);
509 step->SetStepLength(theStep);
510 G4ForceCondition a = NotForced;
511
512 for (G4int iter=0; iter<nbin; iter++) {
513
514 G4double e = emin + ((G4double)iter + 0.5)*bin;
515 gTrack->SetStep(step);
516 gTrack->SetKineticEnergy(e);
517 G4double x = 0.0;
518
519 if(lowE) {
520
521 if(nProcess == 0) x = eionle->GetMeanFreePath(*gTrack,theStep,&a);
522 if(nProcess == 1) x = ebrle->GetMeanFreePath(*gTrack,theStep,&a);
523
524 /*
525 if(nProcess == 2) x = comple->GetMeanFreePath(*gTrack,theStep,&a);
526 if(nProcess == 3) x = convle->GetMeanFreePath(*gTrack,theStep,&a);
527 if(nProcess == 4) x = pele->GetMeanFreePath(*gTrack,theStep,&a);
528 if(nProcess == 5) x = rayle->GetMeanFreePath(*gTrack,theStep,&a);
529 */
530 } else {
531
532 if(nProcess == 0) x = eionst->GetMeanFreePath(*gTrack,theStep,&a);
533 if(nProcess == 1) x = ebrst->GetMeanFreePath(*gTrack,theStep,&a);
534 if(nProcess == 2) x = compst->GetMeanFreePath(*gTrack,theStep,&a);
535 if(nProcess == 3) x = convst->GetMeanFreePath(*gTrack,theStep,&a);
536 if(nProcess == 4) x = pest->GetMeanFreePath(*gTrack,theStep,&a);
537 }
538
539 if(verbose) {
540 G4cout << "Iteration = " << iter
541 << " E = " << e/MeV << " MeV; MeanFreePath= "
542 << x/mm << " mm " << G4endl;
543 }
544
545 if(x > 1000.0*meter) x = 0.0*meter;
546
547 if(usepaw) hist[0]->fill(e/MeV,x/mm);
548 }
549 if(usepaw) {
550 tree->commit();
551 std::cout << "Closing the tree..." << std::endl;
552 tree->close();
553 G4cout << "# hbook is writed" << G4endl;
554 }
555 G4cout << "###### End of run # " << run << " ######" << G4endl;
556
557
558 } while(end);
559 G4cout << "###### End of test #####" << G4endl;
560}
561
562
563
564
565
566
567
568
569
570
571
572
Note: See TracBrowser for help on using the repository browser.