source: trunk/source/processes/electromagnetic/lowenergy/test/G4hLowEnergyTest.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: 13.5 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: G4hLowEnergyTest.cc,v 1.10 2006/06/29 19:44:42 gunter Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// KaonMinusAtRestTest.cc
31// -------------------------------------------------------------------
32//      GEANT 4 class file --- Copyright CERN 1998
33//      CERN Geneva Switzerland
34//
35//
36//      File name:     G4LowEnergyTest
37//
38//      Author:        Christian Voelcker (from M. Maire)
39//
40//      Creation date: ?
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 "G4ProcessManager.hh"
53#include "G4hIonisation.hh"
54#include "G4hLowEnergyIonisation.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 "G4Proton.hh"
62#include "G4AntiProton.hh"
63
64#include "G4Box.hh"
65#include "G4PVPlacement.hh"
66
67#include "G4Step.hh"
68#include "G4GRSVolume.hh"
69
70#include "G4UnitsTable.hh"
71#include "CLHEP/Hist/TupleManager.h"
72#include "CLHEP/Hist/HBookFile.h"
73#include "CLHEP/Hist/Histogram.h"
74#include "CLHEP/Hist/Tuple.h"
75
76HepTupleManager* hbookManager;
77
78main()
79{
80
81  // Setup
82
83  G4int niter=3;
84  G4int imat=4;
85  G4int verboseLevel=1;
86  //G4int processID=9;
87
88  //G4cout << "How many interactions? [10], Which material? [3], which Verbose Level? [1]" << G4endl;
89  //G4cin >> niter >> imat >> verboseLevel;
90
91  //G4cout<<"which process?"<<G4endl<<std::setw(60)<<"[1] = G4LowEnergyPhotoElectric, [2] = G4LowEnergyCompton"<<G4endl;
92  //G4cout<<std::setw(60)<<"[3] = G4LowEnergyRayleigh, [4] = G4LowEnergyGammaconversion"<<G4endl;
93  //G4cout<<std::setw(60)<<"[5] = G4LowEnergyBremstrahlung"<<"[6] = G4LowEnergyIonisation"<<G4endl;
94
95  //G4cin >> processID;
96
97  G4double InitEnergy = 1e-3, InitX = 0., InitY = 0., InitZ = 1.;
98  //G4cout<<"Enter the initial particle energy E and its direction"<<G4endl;
99  //G4cin >> InitEnergy >> InitX >> InitY >> InitZ;
100
101  G4cout.setf( std::ios::scientific, std::ios::floatfield );
102  // -------------------------------------------------------------------
103
104  // ---- HBOOK initialization
105
106
107    hbookManager = new HBookFile("hionisation.hbook", 58);
108    assert (hbookManager != 0);
109 
110  // ---- Book a histogram and ntuples
111  G4cout<<"Hbook file name: "<<((HBookFile*) hbookManager)->filename()<<G4endl;
112 
113  // ---- primary ntuple ------
114  HepTuple* ntuple1 = hbookManager->ntuple("Primary Ntuple");
115  assert (ntuple1 != 0);
116 
117  // ---- secondary ntuple ------
118  HepTuple* ntuple2 = hbookManager->ntuple("Secondary Ntuple");
119  assert (ntuple2 != 0);
120   
121  // ---- secondaries histos ----
122  HepHistogram* hEKin;
123  hEKin = hbookManager->histogram("Kinetic Energy", 100,0.,200.);
124  assert (hEKin != 0); 
125 
126  HepHistogram* hP;
127  hP = hbookManager->histogram("Momentum", 100,0.,1000.);
128  assert (hP != 0); 
129 
130  HepHistogram* hNSec;
131  hNSec = hbookManager->histogram("Number of secondaries", 40,0.,40.);
132  assert (hNSec != 0); 
133 
134  HepHistogram* hDebug;
135  hDebug = hbookManager->histogram("Debug", 100,0.,200.);
136  assert (hDebug != 0); 
137 
138
139  //--------- Materials definition ---------
140
141  G4Material* Be = new G4Material("Beryllium",    4.,  9.01*g/mole, 1.848*g/cm3);
142  G4Material* Graphite = new G4Material("Graphite",6., 12.00*g/mole, 2.265*g/cm3 );
143  G4Material* Al  = new G4Material("Aluminium", 13., 26.98*g/mole, 2.7 *g/cm3);
144  G4Material* Si  = new G4Material("Silicon",   14., 28.055*g/mole, 2.33*g/cm3);
145  G4Material* LAr = new G4Material("LArgon",   18., 39.95*g/mole, 1.393*g/cm3);
146  G4Material* Fe  = new G4Material("Iron",      26., 55.85*g/mole, 7.87*g/cm3);
147  G4Material* Cu  = new G4Material("Copper",    29., 63.55*g/mole, 8.96*g/cm3);
148  G4Material*  W  = new G4Material("Tungsten", 74., 183.85*g/mole, 19.30*g/cm3);
149  G4Material* Pb  = new G4Material("Lead",      82., 207.19*g/mole, 11.35*g/cm3);
150  G4Material*  U  = new G4Material("Uranium", 92., 238.03*g/mole, 18.95*g/cm3);
151
152  G4Element*   H  = new G4Element ("Hydrogen", "H", 1. ,  1.01*g/mole);
153  G4Element*   O  = new G4Element ("Oxygen"  , "O", 8. , 16.00*g/mole);
154  G4Element*   C  = new G4Element ("Carbon"  , "C", 6. , 12.00*g/mole);
155  G4Element*  Cs  = new G4Element ("Cesium"  , "Cs", 55. , 132.905*g/mole);
156  G4Element*   I  = new G4Element ("Iodide"  , "I", 53. , 126.9044*g/mole);
157
158  G4Material*  maO = new G4Material("Oxygen", 8., 16.00*g/mole, 1.1*g/cm3);
159
160  G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2);
161  water->AddElement(H,2);
162  water->AddElement(O,1);
163
164  G4Material* ethane = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
165  ethane->AddElement(H,6);
166  ethane->AddElement(C,2);
167 
168  G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2);
169  csi->AddElement(Cs,1);
170  csi->AddElement(I,1);
171
172  static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
173
174  G4cout<<"The material is: "<<(*theMaterialTable)(imat)->GetName()<<G4endl;
175
176  // Geometry definitions
177  G4Box* theFrame = new G4Box ("Frame",92*mm, 92*mm, 92*mm);
178 
179  G4LogicalVolume* LogicalFrame = new G4LogicalVolume(theFrame,
180                                                      (*theMaterialTable)(imat),
181                                                      "LFrame", 0, 0, 0);
182 
183  G4PVPlacement* PhysicalFrame = new G4PVPlacement(0,G4ThreeVector(),
184                                                   "PFrame",LogicalFrame,0,false,0);
185 
186  // the center-of-mass of the cube should be located at the origin!
187
188  //--------- Particle definition ---------
189
190  G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
191  G4ParticleDefinition* positron = G4Positron::PositronDefinition();
192  G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
193  G4ParticleDefinition* proton = G4Proton::ProtonDefinition();
194  G4ParticleDefinition* antiproton = G4AntiProton::AntiProtonDefinition();
195 
196  //my add 8:45
197  electron->SetCuts(1e-6*mm);
198 
199  //--------- Processes definition ---------
200
201  G4ProcessManager* theProtonProcessManager = new G4ProcessManager(proton);
202  proton->SetProcessManager(theProtonProcessManager);
203
204  G4ProcessManager* theAntiProtonProcessManager = new G4ProcessManager(antiproton);
205  antiproton->SetProcessManager(theAntiProtonProcessManager);
206
207  G4hLowEnergyIonisation* hIonisationProcess = new G4hLowEnergyIonisation;
208  theProtonProcessManager->AddProcess(hIonisationProcess);
209  theProtonProcessManager->SetProcessOrdering(hIonisationProcess,idxAlongStep,1);
210  theProtonProcessManager->SetProcessOrdering(hIonisationProcess,idxPostStep,1);
211
212  // ready for being exploited by antiproton
213  //G4hIonisation hIonisationProcess;
214  //theProtonProcessManager->AddProcess(&hIonisationProcess);
215  //theProtonProcessManager->SetProcessOrdering(&hIonisationProcess,idxAlongStep,1);
216  //theProtonProcessManager->SetProcessOrdering(&hIonisationProcess,idxPostStep,1);
217 
218  G4ForceCondition* condition;
219
220  // ------- set cut and Build CrossSection Tables -------
221  //
222  // -------- create 1 Dynamic Particle  ----
223
224  G4double pEnergy = InitEnergy*MeV;
225
226  G4ParticleMomentum pDirection(InitX,InitY,InitZ);
227 
228  G4DynamicParticle p(G4Proton::Proton(),pDirection,pEnergy);
229  G4DynamicParticle pbar(G4AntiProton::AntiProton(),pDirection,pEnergy);
230 
231
232  //--------- track definition (for this test ONLY!)------------
233
234  G4ThreeVector aPosition(0.,0.,0.);
235  G4double aTime = 0. ;
236
237  G4Track* ptrack;
238  G4Track* ptrackBar;
239
240  ptrack = new G4Track(&p,aTime,aPosition) ;
241 
242  ptrackBar = new G4Track(&pbar,aTime,aPosition) ;
243
244
245  // proton or antiproton?
246
247  G4Track& aTrack = (*ptrack) ;
248 
249  // do I really need this?
250
251  G4GRSVolume* touche = new G4GRSVolume(PhysicalFrame, NULL, aPosition);   
252  ptrack->SetTouchable(touche);
253 
254  // -------- create 1 Step (for this test only)---- 
255
256  G4Step* Step = new G4Step(); 
257  G4Step& aStep = (*Step);
258  Step->SetTrack(ptrack);
259 
260  // --------- check applicability
261  G4ParticleDefinition* ProtonDefinition = p.GetDefinition();
262  G4ParticleDefinition* AntiProtonDefinition = pbar.GetDefinition();
263
264  if(! (hIonisationProcess->IsApplicable(*ProtonDefinition)
265        || hIonisationProcess->IsApplicable(*AntiProtonDefinition)) )
266    {
267      G4Exception("FAIL: *** Not Applicable ***\n");
268    }
269
270  // Initialize the physics tables for ALL processes
271
272  hIonisationProcess->BuildPhysicsTable(*ProtonDefinition);
273  hIonisationProcess->SetPhysicsTableBining(0.5*keV, 3*MeV, 100);
274
275  //hIonisationProcess.BuildPhysicsTable(*antiproton);
276
277
278  G4Material* apttoMaterial ;
279  G4String MaterialName ;
280
281
282  // --------- Test the DoIt for the hIonization
283
284  apttoMaterial = (*theMaterialTable)(imat) ;
285 
286  LogicalFrame->SetMaterial(apttoMaterial); 
287
288  // PostStepDoIt calls
289  G4int iteration = 0;   
290 
291  G4VParticleChange* adummy;
292  G4Track* aFinalParticle;
293  G4String aParticleName; 
294
295  do {
296    adummy = hIonisationProcess->PostStepDoIt(aTrack, aStep);
297 
298    G4ParticleChange* aParticleChange = (G4ParticleChange*) adummy;
299 
300    // ------------ book primary physical quantities -------------
301    G4double pEnChange = 0, pxChange = 0, pyChange = 0, pzChange = 0, pChange = 0;
302
303    pEnChange = aParticleChange->GetEnergyChange();
304    pxChange  = aParticleChange->GetMomentumChange()->x();
305    pyChange  = aParticleChange->GetMomentumChange()->y();
306    pzChange  = aParticleChange->GetMomentumChange()->z();
307    pChange   = std::sqrt(pxChange*pxChange+pyChange*pyChange+pzChange*pzChange);
308
309    // ---- secondaries histos ----   
310    G4cout<<"E and p of the primary particle: "<<pEnChange<<"  "<<pxChange<<"  "
311          <<pyChange<<"  "<<pzChange<<G4endl;
312
313    ntuple1->column("ench", pEnChange);
314    ntuple1->column("pxch", pxChange);
315    ntuple1->column("pych", pyChange);
316    ntuple1->column("pzch", pzChange);
317    ntuple1->column("pch", pChange);
318    ntuple1->dumpData(); 
319   
320    // ------------ book secondaries physical quantities ---------
321   
322    G4double e = 0, eKin = 0, Px = 0, Py = 0, Pz = 0, P = 0, ShID = 0;
323   
324    hNSec->accumulate(aParticleChange->GetNumberOfSecondaries());
325    hDebug->accumulate(aParticleChange->GetLocalEnergyDeposit());
326   
327    for (G4int i = 0; i < (aParticleChange->GetNumberOfSecondaries()); i++) {
328
329      // The following two items should be filled per event, not
330      // per secondary; filled here just for convenience, to avoid
331      // complicated logic to dump ntuple when there are no secondaries
332
333      aFinalParticle = aParticleChange->GetSecondary(i) ;
334     
335      e    = aFinalParticle->GetTotalEnergy();
336      eKin = aFinalParticle->GetKineticEnergy();
337      Px   = (aFinalParticle->GetMomentum()).x();
338      Py   = (aFinalParticle->GetMomentum()).y();
339      Pz   = (aFinalParticle->GetMomentum()).z();
340      P    = std::sqrt(Px*Px+Py*Py+Pz*Pz);
341
342      aParticleName = aFinalParticle->GetDefinition()->GetParticleName();
343      G4cout << aParticleName << ": "
344             << " " << e << "  " 
345             << eKin << "  " 
346             << Px << "  " 
347             << Py << "  "
348             << Pz << " ***" << G4endl;   
349   
350      hEKin->accumulate(eKin);
351      hP->accumulate(std::sqrt(Px*Px+Py*Py+Pz*Pz));
352
353      G4int ptype;
354      if(aParticleName == "proton") ptype = 0;
355      else if(aParticleName == "e-") ptype = -1;
356      else if(aParticleName == "e+") ptype = 1;
357      else if(aParticleName == "antiproton") ptype = 2;
358
359      // Fill the secondaries ntuple
360      ntuple2->column("px", Px);
361      ntuple2->column("py", Py);
362      ntuple2->column("pz", Pz);
363      ntuple2->column("p", P);
364      ntuple2->column("e", e);
365      ntuple2->column("ekin", eKin);
366      ntuple2->column("ptype", ptype);
367
368      ntuple2->dumpData(); 
369
370      delete aParticleChange->GetSecondary(i);
371    }
372   
373    aParticleChange->Clear();
374    iteration++;
375    G4cout << "******* Iteration = " << iteration << G4endl;
376   
377  }  while (iteration < niter) ;
378
379  cout <<"Iteration number: " << G4endl;
380  hbookManager->write();
381  delete hbookManager;
382
383  // delete materials and elements
384  delete Be;
385  delete Graphite;
386  delete Al;
387  delete Si;
388  delete LAr;
389  delete Fe;
390  delete Cu;
391  delete W;
392  delete Pb;
393  delete U;
394  delete H;
395  delete maO;
396  delete C;
397  delete Cs;
398  delete I;
399  delete O;
400  delete water;
401  delete ethane;
402  delete csi;
403  delete Step;
404  delete touche;
405
406  cout<<"END OF THE MAIN PROGRAM"<<G4endl;
407}
408
409
410
411
412
413
414
415
416
417
418
419
420
Note: See TracBrowser for help on using the repository browser.