source: trunk/source/processes/electromagnetic/lowenergy/test/G4mfpTest.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: 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.