source: trunk/source/processes/electromagnetic/lowenergy/test/G4mfpTest.cc @ 1350

Last change on this file since 1350 was 1350, checked in by garnier, 13 years ago

update to last version 4.9.4
File size: 16.4 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// $Id: G4mfpTest.cc,v 1.4 2006/06/29 19:44:48 gunter Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// -------------------------------------------------------------------
31//
32// File name:     G4mfpTest
33//
34// Author:        Maria Grazia Pia
35//
36// Description:   Tests MeanFreePath in electromagnetic processes
37//                Output: ntuple with MeanFreePath and cross sections for
38//                lowenergy and standard processes
39//
40// Modifications:
41// 2 August 2001  MGP            Created
42//
43// -------------------------------------------------------------------
44
45#include "globals.hh"
46#include "G4ios.hh"
47#include <fstream>
48#include <iomanip>
49
50#include "G4Material.hh"
51#include "G4VContinuousDiscreteProcess.hh"
52#include "G4ProcessManager.hh"
53#include "G4VeLowEnergyLoss.hh"
54#include "G4EnergyLossTables.hh"
55#include "G4VParticleChange.hh"
56#include "G4ParticleChange.hh"
57#include "G4DynamicParticle.hh"
58#include "G4Electron.hh"
59#include "G4Positron.hh"
60#include "G4Gamma.hh"
61#include "G4Box.hh"
62#include "G4PVPlacement.hh"
63#include "G4Step.hh"
64#include "G4GRSVolume.hh"
65#include "Randomize.hh"
66#include "G4UnitsTable.hh"
67
68#include "G4LowEnergyCompton.hh"
69#include "G4LowEnergyPolarizedCompton.hh"
70#include "G4ComptonScattering.hh"
71
72#include "G4LowEnergyPhotoElectric.hh"
73#include "G4PhotoElectricEffect.hh"
74
75#include "G4LowEnergyGammaConversion.hh"
76#include "G4GammaConversion.hh"
77
78#include "G4LowEnergyRayleigh.hh"
79#include "G4LowEnergyBremsstrahlung.hh"
80#include "G4eBremsstrahlung.hh"
81
82#include "G4LowEnergyIonisation.hh"
83#include "G4eIonisation.hh"
84
85#include "CLHEP/Hist/TupleManager.h"
86#include "CLHEP/Hist/HBookFile.h"
87#include "CLHEP/Hist/Histogram.h"
88#include "CLHEP/Hist/Tuple.h"
89
90#include "G4LowEnergyGammaConversionMG.hh"
91#include "G4LowEnergyRayleighMG.hh"
92#include "G4LowEnergyComptonMG.hh"
93#include "G4LowEnergyPhotoElectricMG.hh"
94
95
96HepTupleManager* hbookManager;
97
98int main()
99{
100  // Setup
101
102  G4int nIterations = 100000;
103  G4int materialId = 3;
104  G4int test = 0;
105
106  G4cout.setf( ios::scientific, ios::floatfield );
107
108  // -------------------------------------------------------------------
109
110  // ---- HBOOK initialization
111
112
113  hbookManager = new HBookFile("mfp.hbook", 58);
114  assert (hbookManager != 0);
115 
116  // ---- Book a histogram and ntuples
117  G4cout<<"Hbook file name: "<<((HBookFile*) hbookManager)->filename()<<endl;
118 
119  // ---- primary ntuple ------
120  HepTuple* ntuple1 = hbookManager->ntuple("MeanFreePath Ntuple");
121  assert (ntuple1 != 0);
122 
123  //--------- Material definitions ---------
124
125  G4Material* Be = new G4Material("Beryllium",    4.,  9.01*g/mole, 1.848*g/cm3);
126  G4Material* Graphite = new G4Material("Graphite",6., 12.00*g/mole, 2.265*g/cm3 );
127  G4Material* Al  = new G4Material("Aluminium", 13., 26.98*g/mole, 2.7 *g/cm3);
128  G4Material* Si  = new G4Material("Silicon",   14., 28.055*g/mole, 2.33*g/cm3);
129  G4Material* LAr = new G4Material("LArgon",   18., 39.95*g/mole, 1.393*g/cm3);
130  G4Material* Fe  = new G4Material("Iron",      26., 55.85*g/mole, 7.87*g/cm3);
131  G4Material* Cu  = new G4Material("Copper",    29., 63.55*g/mole, 8.96*g/cm3);
132  G4Material*  W  = new G4Material("Tungsten", 74., 183.85*g/mole, 19.30*g/cm3);
133  G4Material* Pb  = new G4Material("Lead",      82., 207.19*g/mole, 11.35*g/cm3);
134  G4Material*  U  = new G4Material("Uranium", 92., 238.03*g/mole, 18.95*g/cm3);
135
136  G4Element*   H  = new G4Element ("Hydrogen", "H", 1. ,  1.01*g/mole);
137  G4Element*   O  = new G4Element ("Oxygen"  , "O", 8. , 16.00*g/mole);
138  G4Element*   C  = new G4Element ("Carbon"  , "C", 6. , 12.00*g/mole);
139  G4Element*  Cs  = new G4Element ("Cesium"  , "Cs", 55. , 132.905*g/mole);
140  G4Element*   I  = new G4Element ("Iodide"  , "I", 53. , 126.9044*g/mole);
141  G4Element*   N  = new G4Element("Nitrogen" , "N" , 7., 14.01*g/mole);
142
143  G4Material*  maO = new G4Material("Oxygen", 8., 16.00*g/mole, 1.1*g/cm3);
144
145  G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2);
146  water->AddElement(H,2);
147  water->AddElement(O,1);
148
149  G4Material* ethane = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
150  ethane->AddElement(H,6);
151  ethane->AddElement(C,2);
152 
153  G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2);
154  csi->AddElement(Cs,1);
155  csi->AddElement(I,1);
156
157  G4double fractionMass;
158  G4Material* air = new G4Material("Air", 1.290*mg/cm3, 2);
159  air->AddElement(N, fractionMass=0.7);
160  air->AddElement(O, fractionMass=0.3);
161
162  // Interactive set-up
163
164  G4cout << "How many interactions? " << G4endl;
165  G4cin >> nIterations;
166
167  if (nIterations <= 0) G4Exception("Wrong input");
168
169  G4cout << "1) Compton 2) Photoelectric 3) GammaConversion 4) Rayleigh" << G4endl;
170  G4cout << "5) Bremsstrahlung 6) Ionisation" << G4endl;
171  G4int type;
172  G4cin >> type;
173
174  if (nIterations <= 0) G4Exception("Wrong input");
175  if (type < 1 || type > 6) G4Exception("Wrong input");
176
177  G4double initEnergy = 1*MeV; 
178  G4double initX = 0.; 
179  G4double initY = 0.; 
180  G4double initZ = 1.;
181
182  static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
183
184 G4int nMaterials = G4Material::GetNumberOfMaterials();
185
186  G4cout << "Available materials are: " << G4endl;
187  for (G4int mat = 0; mat < nMaterials; mat++)
188    {
189      G4cout << mat << ") "
190             << (*theMaterialTable)[mat]->GetName()
191             << G4endl;
192    }
193  G4cout << "Which material? " << G4endl;
194  G4cin >> materialId;
195
196  G4Material* material = (*theMaterialTable)[materialId] ;
197
198  G4cout << "The selected material is: "
199         << material->GetName()
200         << G4endl;
201
202  G4double dimX = 1*mm;
203  G4double dimY = 1*mm;
204  G4double dimZ = 1*mm;
205 
206  // Geometry
207
208  G4Box* theFrame = new G4Box ("Frame",dimX, dimY, dimZ);
209 
210  G4LogicalVolume* logicalFrame = new G4LogicalVolume(theFrame,
211                                                      (*theMaterialTable)[materialId],
212                                                      "LFrame", 0, 0, 0);
213  logicalFrame->SetMaterial(material); 
214 
215  G4PVPlacement* physicalFrame = new G4PVPlacement(0,G4ThreeVector(),
216                                                   "PFrame",logicalFrame,0,false,0);
217 
218  // Particle definitions
219
220  G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
221  G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
222  G4ParticleDefinition* positron = G4Positron::PositronDefinition();
223 
224  gamma->SetCuts(1e-3*mm);
225  electron->SetCuts(1e-3*mm);
226  positron->SetCuts(1e-3*mm);
227 
228  // Electrons
229  G4ProcessManager* eProcessManager = new G4ProcessManager(electron);
230  electron->SetProcessManager(eProcessManager);
231
232 // Bremsstrahlung
233  G4LowEnergyBremsstrahlung* bremsstrahlung = new G4LowEnergyBremsstrahlung;
234  eProcessManager->AddProcess(bremsstrahlung);
235  bremsstrahlung->BuildPhysicsTable(*electron);
236  G4eBremsstrahlung* bremsstrahlungStd = new G4eBremsstrahlung;
237  eProcessManager->AddProcess(bremsstrahlungStd);
238  bremsstrahlungStd->BuildPhysicsTable(*electron);
239
240  // Ionisation
241  G4LowEnergyIonisation* ionisation = new G4LowEnergyIonisation;
242  eProcessManager->AddProcess(ionisation);
243  ionisation->BuildPhysicsTable(*electron);
244  G4eIonisation* ionisationStd = new G4eIonisation;
245  eProcessManager->AddProcess(ionisationStd);
246  ionisationStd->BuildPhysicsTable(*electron);
247
248  // Positrons   
249  G4ProcessManager* positronProcessManager = new G4ProcessManager(positron);
250  positron->SetProcessManager(positronProcessManager);
251  positronProcessManager->AddProcess(bremsstrahlung);
252 
253  // Photons
254  G4ProcessManager* gProcessManager = new G4ProcessManager(gamma);
255  gamma->SetProcessManager(gProcessManager);
256
257  // Compton
258  G4LowEnergyCompton* compton = new G4LowEnergyCompton;
259  G4LowEnergyComptonMG* comptonMG = new G4LowEnergyComptonMG;
260  G4ComptonScattering* comptonStd = new G4ComptonScattering;
261  gProcessManager->AddProcess(compton);
262  compton->BuildPhysicsTable(*gamma);
263  gProcessManager->AddProcess(comptonMG);
264  comptonMG->BuildPhysicsTable(*gamma);
265  gProcessManager->AddProcess(comptonStd);
266  comptonStd->BuildPhysicsTable(*gamma);
267
268  // Photoelectric
269  G4LowEnergyPhotoElectric* photoelectric = new G4LowEnergyPhotoElectric;
270  gProcessManager->AddProcess(photoelectric);
271  photoelectric->BuildPhysicsTable(*gamma);
272  G4PhotoElectricEffect* photoelectricStd = new G4PhotoElectricEffect;
273  gProcessManager->AddProcess(photoelectricStd);
274  photoelectricStd->BuildPhysicsTable(*gamma);
275  G4LowEnergyPhotoElectricMG* photoelectricMG = new G4LowEnergyPhotoElectricMG;
276  gProcessManager->AddProcess(photoelectricMG);
277  photoelectricMG->BuildPhysicsTable(*gamma);
278
279  // GammaConversion
280  G4LowEnergyGammaConversion* gammaConversion = new G4LowEnergyGammaConversion;
281  gProcessManager->AddProcess(gammaConversion);
282  gammaConversion->BuildPhysicsTable(*gamma);
283  G4GammaConversion* gammaConversionStd = new G4GammaConversion;
284  gProcessManager->AddProcess(gammaConversionStd);
285  gammaConversionStd->BuildPhysicsTable(*gamma);
286  G4LowEnergyGammaConversionMG* gammaConversionMG = new G4LowEnergyGammaConversionMG;
287  gProcessManager->AddProcess(gammaConversionMG);
288  gammaConversionMG->BuildPhysicsTable(*gamma);
289
290  // Rayleigh
291  G4LowEnergyRayleigh* rayleigh = new G4LowEnergyRayleigh;
292  gProcessManager->AddProcess(rayleigh);
293  rayleigh->BuildPhysicsTable(*gamma);
294  G4LowEnergyRayleighMG* rayleighMG = new G4LowEnergyRayleighMG;
295  gProcessManager->AddProcess(rayleighMG);
296  rayleighMG->BuildPhysicsTable(*gamma);
297
298  // --------- Test GetMeanFreePath
299
300  G4cout << "Mean Free Path in material " << material->GetName() << G4endl;
301
302  G4double eMin = 250. * eV;
303  G4double eMax = 200. * GeV;
304  G4double diff = eMax - eMin;
305
306  for (G4int iter=0; iter<nIterations; iter++)
307    {
308      G4double range = G4UniformRand();
309      G4double gEnergy;
310      if (range < 0.1)
311        {
312          eMin = 250. * eV;
313          eMax = 10. * keV;
314          diff = eMax - eMin;
315          G4double random =  G4UniformRand();
316          gEnergy = eMin + diff * random;
317        }
318      else if (range < 0.2)
319        {
320          eMin = 10. * keV;
321          eMax = 100. * keV;
322          diff = eMax - eMin;
323          G4double random =  G4UniformRand();
324          gEnergy = eMin + diff * random;
325        }
326      else if (range < 0.3)
327        {
328          eMin = 100. * keV;
329          eMax = 1. * MeV;
330          diff = eMax - eMin;
331          G4double random =  G4UniformRand();
332          gEnergy = eMin + diff * random;
333        }
334      else if (range < 0.4)
335        {
336          eMin = 1. * MeV;
337          eMax = 10. * MeV;
338          diff = eMax - eMin;
339          G4double random =  G4UniformRand();
340          gEnergy = eMin + diff * random;
341        }
342      else if (range < 0.5)
343        {
344          eMin = 10. * MeV;
345          eMax = 100. * MeV;
346          diff = eMax - eMin;
347          G4double random =  G4UniformRand();
348          gEnergy = eMin + diff * random;
349        }
350      else if (range < 0.6)
351        {
352          eMin = 100. * MeV;
353          eMax = 1. * GeV;
354          diff = eMax - eMin;
355          G4double random =  G4UniformRand();
356          gEnergy = eMin + diff * random;
357        }
358      else if (range < 0.8)
359        {
360          eMin = 1. * GeV;
361          eMax = 100. * GeV;
362          diff = eMax - eMin;
363          G4double random =  G4UniformRand();
364          gEnergy = eMin + diff * random;
365        }
366         else
367        {
368          gEnergy = eMin + G4UniformRand();
369        }
370
371      if (type == 3 && gEnergy < 1.0220*MeV) gEnergy = 1.02200 * MeV;
372
373      // Create a DynamicParticle 
374      G4ParticleMomentum gDirection(initX,initY,initZ);
375      G4DynamicParticle dynamicPhoton(G4Gamma::Gamma(),gDirection,gEnergy);
376     
377      // Track
378     
379      G4ThreeVector aPosition(0.,0.,0.);
380      G4ThreeVector newPosition(0.,0.,1.*mm);
381      G4double aTime = 0. ;
382     
383      G4Track* gTrack = new G4Track(&dynamicPhoton,aTime,aPosition);
384     
385      // do I really need this?
386     
387      G4GRSVolume* touche = new G4GRSVolume(physicalFrame, 0, aPosition);   
388      gTrack->SetTouchable(touche);
389     
390      // Step
391     
392      G4Step* step = new G4Step(); 
393      step->SetTrack(gTrack);     
394      G4StepPoint* aPoint = new G4StepPoint();
395      aPoint->SetPosition(aPosition);
396      aPoint->SetMaterial(material);
397      G4double safety = 100.*cm;
398      aPoint->SetSafety(safety);
399      step->SetPreStepPoint(aPoint);
400      G4StepPoint* newPoint = new G4StepPoint();
401      newPoint->SetPosition(newPosition);
402      newPoint->SetMaterial(material);
403      newPoint->SetSafety(safety);
404      step->SetPostStepPoint(newPoint);
405      step->SetStepLength(1*micrometer);
406      gTrack->SetStep(step); 
407
408      G4ForceCondition* force = new G4ForceCondition;
409
410      G4double mfp;
411      G4double mfpStd;
412      G4double mfpMG = 0.;
413
414      if (type == 1)
415        {
416          mfp = compton->GetMeanFreePath(*gTrack,0.,force);
417          mfpStd = comptonStd->GetMeanFreePath(*gTrack,0.,force);
418          mfpMG = comptonMG->DumpMeanFreePath(*gTrack,0.,force);
419        }
420      else if (type == 2)
421        {
422          mfp = photoelectric->GetMeanFreePath(*gTrack,0.,force);
423          mfpStd = photoelectricStd->GetMeanFreePath(*gTrack,0.,force);
424          mfpMG = photoelectricMG->DumpMeanFreePath(*gTrack,0.,force);
425        }
426        else if (type == 3)
427        {
428          mfp = gammaConversion->GetMeanFreePath(*gTrack,0.,force);
429          mfpStd = gammaConversionStd->GetMeanFreePath(*gTrack,0.,force);
430          mfpMG = gammaConversionMG->DumpMeanFreePath(*gTrack,0.,force);
431        }
432      else if (type == 4)
433        {
434          mfp = rayleigh->GetMeanFreePath(*gTrack,0.,force);
435          //  mfpStd = rayleighStd->GetMeanFreePath(*gTrack,0.,force);
436          mfpMG = rayleighMG->DumpMeanFreePath(*gTrack,0.,force);
437        }
438      else if (type == 5)
439        {
440          mfp = bremsstrahlung->GetMeanFreePath(*gTrack,0.,force);
441          mfpStd = bremsstrahlungStd->GetMeanFreePath(*gTrack,0.,force);
442          //  mfpMG = bremsstrahlungMG->DumpMeanFreePath(*gTrack,0.,force);
443        }
444      else 
445        {
446          mfp = ionisation->GetMeanFreePath(*gTrack,0.,force);
447          mfpStd = ionisationStd->GetMeanFreePath(*gTrack,0.,force);
448          //  mfpMG = ionisationMG->DumpMeanFreePath(*gTrack,0.,force);
449        }
450
451      G4double sigma = 0.;
452      if (mfp > 0.) sigma = 1./mfp;
453      G4double sigmaMG = 0.;
454      if (mfpMG > 0.) sigmaMG = 1./mfpMG;
455      G4double sigmaStd = 0.;
456      if (mfpStd > 0.) sigmaStd = 1./mfpStd;
457
458      G4double kineticEnergy = gTrack->GetKineticEnergy();
459
460      G4cout << iter << ") e = " << gEnergy
461             << " - mfp = " << mfp
462             << ";  mfpMG = " << mfpMG << G4endl;
463
464      ntuple1->column("e", gEnergy);
465      ntuple1->column("ekin", kineticEnergy);
466      ntuple1->column("mfp", mfp);
467      ntuple1->column("mfpmg", mfpMG); 
468      ntuple1->column("mfpstd", mfpStd); 
469      ntuple1->column("sigma", sigma);
470      ntuple1->column("sigmamg", sigmaMG); 
471      ntuple1->column("sigmastd", sigmaStd);
472      ntuple1->dumpData();     
473    } 
474 
475
476  cout  << "End of iteration "  <<  G4endl;
477  hbookManager->write();
478  delete hbookManager;
479 
480  // delete materials and elements
481  //  delete Be;
482  //  delete Graphite;
483  //  delete Al;
484  //  delete Si;
485  //  delete LAr;
486  //  delete Fe;
487  //  delete Cu;
488  //  delete W;
489  //  delete Pb;
490  //  delete U;
491  //  delete H;
492  //  delete maO;
493  //  delete C;
494  //  delete Cs;
495  //  delete I;
496  //  delete O;
497  //  delete water;
498  //  delete ethane;
499  //  delete csi;
500  //  delete step;
501  //  delete touche;
502  //  delete Be;
503  //  delete Graphite;
504  //  delete Al;
505  //  delete Si;
506  //  delete LAr;
507  //  delete Fe;
508  //  delete Cu;
509  //  delete W;
510  //  delete Pb;
511  //  delete U;
512  //  delete H;
513  //  delete maO;
514  //  delete C;
515  //  delete Cs;
516  //  delete I;
517  //  delete O;
518  //  delete water;
519  //  delete ethane;
520  //  delete csi;
521  // delete step;
522  //  delete touche;
523
524  cout << "END OF THE MAIN PROGRAM" << G4endl;
525}
526
527
528
529
530
531
532
533
534
535
536
537
Note: See TracBrowser for help on using the repository browser.