source: trunk/source/processes/electromagnetic/lowenergy/test/G4RayleighTest.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: 15.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// $Id: G4RayleighTest.cc,v 1.10 2006/06/29 19:44:32 gunter Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// -------------------------------------------------------------------
31//      GEANT 4 class file --- Copyright CERN 1998
32//      CERN Geneva Switzerland
33//
34//
35//      File name:     G4RayleighTest
36//
37//      Author:        Maria Grazia Pia
38//
39//      Creation date: 2 May 2001
40//
41//      Modifications:
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 "G4LowEnergyBremsstrahlung.hh"
54#include "G4eBremsstrahlung.hh"
55#include "G4LowEnergyRayleigh.hh"
56#include "G4LowEnergyIonisation.hh"
57#include "G4eIonisation.hh"
58#include "G4VeLowEnergyLoss.hh"
59#include "G4EnergyLossTables.hh"
60#include "G4VParticleChange.hh"
61#include "G4ParticleChange.hh"
62#include "G4DynamicParticle.hh"
63#include "G4Electron.hh"
64#include "G4Positron.hh"
65#include "G4Gamma.hh"
66
67#include "G4Box.hh"
68#include "G4PVPlacement.hh"
69
70#include "G4Step.hh"
71#include "G4GRSVolume.hh"
72
73#include "G4UnitsTable.hh"
74#include "CLHEP/Hist/TupleManager.h"
75#include "CLHEP/Hist/HBookFile.h"
76#include "CLHEP/Hist/Histogram.h"
77#include "CLHEP/Hist/Tuple.h"
78
79HepTupleManager* hbookManager;
80
81int main()
82{
83
84  // Setup
85
86  G4int nIterations = 100000;
87  G4int materialId = 3;
88  G4int test = 0;
89
90  G4cout.setf( ios::scientific, ios::floatfield );
91
92  // -------------------------------------------------------------------
93
94  // ---- HBOOK initialization
95
96
97  hbookManager = new HBookFile("rayleightest.hbook", 58);
98  assert (hbookManager != 0);
99 
100  // ---- Book a histogram and ntuples
101  G4cout<<"Hbook file name: "<<((HBookFile*) hbookManager)->filename()<<endl;
102 
103  // ---- primary ntuple ------
104  HepTuple* ntuple1 = hbookManager->ntuple("Primary Ntuple");
105  assert (ntuple1 != 0);
106 
107  // ---- secondary ntuple ------
108  HepTuple* ntuple2 = hbookManager->ntuple("Secondary Ntuple");
109  assert (ntuple2 != 0);
110   
111  // ---- secondaries histos ----
112  HepHistogram* hEKin;
113  hEKin = hbookManager->histogram("Kinetic Energy", 100,0.,200.);
114  assert (hEKin != 0); 
115 
116  HepHistogram* hP;
117  hP = hbookManager->histogram("Momentum", 100,0.,1000.);
118  assert (hP != 0); 
119 
120  HepHistogram* hNSec;
121  hNSec = hbookManager->histogram("Number of secondaries", 40,0.,40.);
122  assert (hNSec != 0); 
123 
124  HepHistogram* hDebug;
125  hDebug = hbookManager->histogram("Debug", 100,0.,200.);
126  assert (hDebug != 0); 
127 
128
129  //--------- Materials definition ---------
130
131  G4Material* Be = new G4Material("Beryllium",    4.,  9.01*g/mole, 1.848*g/cm3);
132  G4Material* Graphite = new G4Material("Graphite",6., 12.00*g/mole, 2.265*g/cm3 );
133  G4Material* Al  = new G4Material("Aluminium", 13., 26.98*g/mole, 2.7 *g/cm3);
134  G4Material* Si  = new G4Material("Silicon",   14., 28.055*g/mole, 2.33*g/cm3);
135  G4Material* LAr = new G4Material("LArgon",   18., 39.95*g/mole, 1.393*g/cm3);
136  G4Material* Fe  = new G4Material("Iron",      26., 55.85*g/mole, 7.87*g/cm3);
137  G4Material* Cu  = new G4Material("Copper",    29., 63.55*g/mole, 8.96*g/cm3);
138  G4Material*  W  = new G4Material("Tungsten", 74., 183.85*g/mole, 19.30*g/cm3);
139  G4Material* Pb  = new G4Material("Lead",      82., 207.19*g/mole, 11.35*g/cm3);
140  G4Material*  U  = new G4Material("Uranium", 92., 238.03*g/mole, 18.95*g/cm3);
141
142  G4Element*   H  = new G4Element ("Hydrogen", "H", 1. ,  1.01*g/mole);
143  G4Element*   O  = new G4Element ("Oxygen"  , "O", 8. , 16.00*g/mole);
144  G4Element*   C  = new G4Element ("Carbon"  , "C", 6. , 12.00*g/mole);
145  G4Element*  Cs  = new G4Element ("Cesium"  , "Cs", 55. , 132.905*g/mole);
146  G4Element*   I  = new G4Element ("Iodide"  , "I", 53. , 126.9044*g/mole);
147
148  G4Material*  maO = new G4Material("Oxygen", 8., 16.00*g/mole, 1.1*g/cm3);
149
150  G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2);
151  water->AddElement(H,2);
152  water->AddElement(O,1);
153
154  G4Material* ethane = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
155  ethane->AddElement(H,6);
156  ethane->AddElement(C,2);
157 
158  G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2);
159  csi->AddElement(Cs,1);
160  csi->AddElement(I,1);
161
162
163  // Interactive set-up
164
165  G4cout << "How many interactions? " << G4endl;
166  G4cin >> nIterations;
167
168  if (nIterations <= 0) G4Exception("Wrong input");
169
170  G4double initEnergy = 1*MeV; 
171  G4double initX = 0.; 
172  G4double initY = 0.; 
173  G4double initZ = 1.;
174
175  G4cout << "Enter the initial particle energy E (MeV)" << G4endl; 
176  G4cin >> initEnergy ;
177
178  initEnergy = initEnergy * MeV;
179
180  if (initEnergy  <= 0.) G4Exception("Wrong input");
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
194  G4cout << "Which material? " << G4endl;
195  G4cin >> materialId;
196
197  G4Material* material = (*theMaterialTable)[materialId] ;
198
199  G4cout << "The selected material is: "
200         << material->GetName()
201         << G4endl;
202
203  G4double dimX = 1*mm;
204  G4double dimY = 1*mm;
205  G4double dimZ = 1*mm;
206 
207  // Geometry
208
209  G4Box* theFrame = new G4Box ("Frame",dimX, dimY, dimZ);
210 
211  G4LogicalVolume* logicalFrame = new G4LogicalVolume(theFrame,
212                                                      (*theMaterialTable)[materialId],
213                                                      "LFrame", 0, 0, 0);
214  logicalFrame->SetMaterial(material); 
215 
216  G4PVPlacement* physicalFrame = new G4PVPlacement(0,G4ThreeVector(),
217                                                   "PFrame",logicalFrame,0,false,0);
218 
219  // Particle definitions
220
221  G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
222  G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
223  G4ParticleDefinition* positron = G4Positron::PositronDefinition();
224 
225  gamma->SetCuts(1e-3*mm);
226  electron->SetCuts(1e-3*mm);
227  positron->SetCuts(1e-3*mm);
228 
229  // Processes
230
231  G4VContinuousDiscreteProcess* bremProcess = new G4LowEnergyBremsstrahlung;
232  G4VContinuousDiscreteProcess* ioniProcess = new G4LowEnergyIonisation;
233
234  G4ProcessManager* eProcessManager = new G4ProcessManager(electron);
235  electron->SetProcessManager(eProcessManager);
236  eProcessManager->AddProcess(bremProcess);
237   
238  G4ProcessManager* positronProcessManager = new G4ProcessManager(positron);
239  positron->SetProcessManager(positronProcessManager);
240  positronProcessManager->AddProcess(bremProcess);
241 
242  // Initialize the physics tables
243  bremProcess->BuildPhysicsTable(*electron);
244  ioniProcess->BuildPhysicsTable(*electron);
245
246  // Photon process
247  G4VDiscreteProcess* photonProcess = new G4LowEnergyRayleigh;
248
249  G4ProcessManager* gProcessManager = new G4ProcessManager(gamma);
250  gamma->SetProcessManager(gProcessManager);
251  gProcessManager->AddProcess(photonProcess);
252  photonProcess->BuildPhysicsTable(*gamma);
253
254  // Create a DynamicParticle 
255 
256  G4double gEnergy = initEnergy*MeV;
257  G4ParticleMomentum gDirection(initX,initY,initZ);
258  G4DynamicParticle dynamicPhoton(G4Gamma::Gamma(),gDirection,gEnergy);
259
260  // Track
261
262  G4ThreeVector aPosition(0.,0.,0.);
263  G4ThreeVector newPosition(0.,0.,1.*mm);
264  G4double aTime = 0. ;
265
266  G4Track* gTrack = new G4Track(&dynamicPhoton,aTime,aPosition);
267
268  // do I really need this?
269
270  G4GRSVolume* touche = new G4GRSVolume(physicalFrame, 0, aPosition);   
271  gTrack->SetTouchable(touche);
272 
273 // Step
274
275  G4Step* step = new G4Step(); 
276  step->SetTrack(gTrack);
277
278  G4StepPoint* aPoint = new G4StepPoint();
279  aPoint->SetPosition(aPosition);
280  aPoint->SetMaterial(material);
281  G4double safety = 10000.*cm;
282  aPoint->SetSafety(safety);
283  step->SetPreStepPoint(aPoint);
284  G4StepPoint* newPoint = new G4StepPoint();
285  newPoint->SetPosition(newPosition);
286  newPoint->SetMaterial(material);
287  step->SetPostStepPoint(newPoint);
288
289  // Check applicability
290 
291  if (! (photonProcess->IsApplicable(*gamma)))
292    {
293      G4Exception("Not Applicable");
294    }
295
296  // --------- Test the DoIt
297
298  G4cout << "DoIt in material " << material->GetName() << G4endl;
299
300  G4int iteration = 0;   
301 
302
303  for (G4int iter=0; iter<nIterations; iter++)
304    {
305      step->SetStepLength(1*micrometer);
306
307      gTrack->SetStep(step); 
308 
309      G4cout  <<  "Iteration = "  <<  iter
310              << "  -  Step Length = " 
311              << step->GetStepLength()/mm << " mm "
312              << G4endl;
313
314      G4VParticleChange* dummy;
315      dummy = photonProcess->PostStepDoIt(*gTrack, *step);
316
317      G4ParticleChange* particleChange = (G4ParticleChange*) dummy;
318     
319      // Primary physical quantities
320
321      G4double energyChange = particleChange->GetEnergyChange();
322      G4double dedx = initEnergy - energyChange ;
323      G4double dedxNow = dedx / (step->GetStepLength());
324     
325      G4ThreeVector eChange = particleChange->CalcMomentum(energyChange,
326                                                           (*particleChange->GetMomentumChange()),
327                                                           particleChange->GetMassChange());
328
329      G4double pxChange  = eChange.x();
330      G4double pyChange  = eChange.y();
331      G4double pzChange  = eChange.z();
332      G4double pChange   = std::sqrt(pxChange*pxChange + pyChange*pyChange + pzChange*pzChange);
333     
334      G4double xChange = particleChange->GetPositionChange()->x();
335      G4double yChange = particleChange->GetPositionChange()->y();
336      G4double zChange = particleChange->GetPositionChange()->z();
337      G4double thetaChange = particleChange->GetPositionChange()->theta();
338
339      G4cout << "---- Primary after the step ---- " << G4endl;
340 
341      //      G4cout << "Position (x,y,z) = "
342      //             << xChange << "  "
343      //             << yChange << "   "
344      //             << zChange << "   "
345      //             << G4endl;
346
347      G4cout << "---- Energy: " << energyChange/MeV << " MeV,  " 
348             << "(px,py,pz): ("
349             << pxChange/MeV << ","
350             << pyChange/MeV << "," 
351             << pzChange/MeV << ") MeV"
352             << G4endl;
353
354      G4cout << "---- Energy loss (dE) = " << dedx/keV << " keV" << G4endl;
355      //      G4cout << "Stopping power (dE/dx)=" << dedxNow << G4endl;
356     
357      // Secondaries
358
359      ntuple1->column("eprimary", initEnergy);
360      ntuple1->column("energyf", energyChange);
361      ntuple1->column("de", dedx);
362      ntuple1->column("dedx", dedxNow);
363      ntuple1->column("pxch", xChange);
364      ntuple1->column("pych", pyChange);
365      ntuple1->column("pzch", pzChange);
366      ntuple1->column("pch", zChange); 
367      ntuple1->column("thetach", thetaChange); 
368     
369      // Secondaries physical quantities
370     
371      hNSec->accumulate(particleChange->GetNumberOfSecondaries());
372      hDebug->accumulate(particleChange->GetLocalEnergyDeposit());
373     
374      G4int nElectrons = 0;
375      G4int nPositrons = 0;
376      G4int nPhotons = 0;
377
378      for (G4int i = 0; i < (particleChange->GetNumberOfSecondaries()); i++) 
379        {
380          // The following two items should be filled per event, not
381          // per secondary; filled here just for convenience, to avoid
382          // complicated logic to dump ntuple when there are no secondaries
383         
384          G4Track* finalParticle = particleChange->GetSecondary(i) ;
385         
386          G4double e    = finalParticle->GetTotalEnergy();
387          G4double eKin = finalParticle->GetKineticEnergy();
388          G4double px   = (finalParticle->GetMomentum()).x();
389          G4double py   = (finalParticle->GetMomentum()).y();
390          G4double pz   = (finalParticle->GetMomentum()).z();
391          G4double theta   = (finalParticle->GetMomentum()).theta();
392          G4double p   = std::sqrt(px*px+py*py+pz*pz);
393
394          if (e > initEnergy)
395            {
396              G4cout << "WARNING: eFinal > eInit " << G4endl;
397                //           << e
398                //                   << " > " initEnergy
399                 
400            }
401
402          G4String particleName = finalParticle->GetDefinition()->GetParticleName();
403          G4cout  << "==== Final " 
404                  <<  particleName  <<  " " 
405                  << "energy: " <<  e/MeV  <<  " MeV,  " 
406                  << "eKin: " <<  eKin/MeV  <<  " MeV, " 
407                  << "(px,py,pz): ("
408                  <<  px/MeV  <<  "," 
409                  <<  py/MeV  <<  ","
410                  <<  pz/MeV  << ") MeV "
411                  <<  G4endl;   
412         
413          hEKin->accumulate(eKin);
414          hP->accumulate(p);
415         
416          G4int partType = 0;
417          if (particleName == "e-") 
418            {
419              partType = 1;
420              nElectrons++;
421            }
422          else if (particleName == "e+") 
423            {
424              partType = 2;
425              nPositrons++;
426            }
427          else if (particleName == "gamma") 
428            {
429              partType = 3;
430              nPhotons++;
431            }
432          // Fill the secondaries ntuple
433          ntuple2->column("eprimary",initEnergy);
434          ntuple2->column("px", px);
435          ntuple2->column("py", py);
436          ntuple2->column("pz", pz);
437          ntuple2->column("p", p);
438          ntuple2->column("e", e);
439          ntuple2->column("theta", theta);
440          ntuple2->column("ekin", eKin);
441          ntuple2->column("type", partType);
442         
443          ntuple2->dumpData(); 
444         
445          delete particleChange->GetSecondary(i);
446        }
447
448      ntuple1->column("nelectrons",nElectrons);
449      ntuple1->column("npositrons",nPositrons);
450      ntuple1->column("nphotons",nPhotons);
451      ntuple1->dumpData(); 
452                 
453      particleChange->Clear();
454     
455    } 
456 
457 
458  cout  << "Iteration number: "  <<  G4endl;
459  hbookManager->write();
460  delete hbookManager;
461 
462  // delete materials and elements
463  //  delete Be;
464  //  delete Graphite;
465  //  delete Al;
466  //  delete Si;
467  //  delete LAr;
468  //  delete Fe;
469  //  delete Cu;
470  //  delete W;
471  //  delete Pb;
472  //  delete U;
473  //  delete H;
474  //  delete maO;
475  //  delete C;
476  //  delete Cs;
477  //  delete I;
478  //  delete O;
479  //  delete water;
480  //  delete ethane;
481  //  delete csi;
482  //  delete step;
483  //  delete touche;
484  //  delete Be;
485  //  delete Graphite;
486  //  delete Al;
487  //  delete Si;
488  //  delete LAr;
489  //  delete Fe;
490  //  delete Cu;
491  //  delete W;
492  //  delete Pb;
493  //  delete U;
494  //  delete H;
495  //  delete maO;
496  //  delete C;
497  //  delete Cs;
498  //  delete I;
499  //  delete O;
500  //  delete water;
501  //  delete ethane;
502  //  delete csi;
503  delete step;
504  //  delete touche;
505
506  cout << "END OF THE MAIN PROGRAM" << G4endl;
507}
508
509
510
511
512
513
514
515
516
517
518
519
Note: See TracBrowser for help on using the repository browser.