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

Last change on this file since 1228 was 1199, checked in by garnier, 15 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.