source: trunk/source/processes/electromagnetic/lowenergy/test/G4hLowEnergyTest.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: 13.5 KB
RevLine 
[1350]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.