source: trunk/source/processes/electromagnetic/lowenergy/test/G4DNAProcessTest.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: 8.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: G4DNAProcessTest.cc,v 1.6 2007/10/15 09:00:55 pia Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30///
31// -------------------------------------------------------------------
32// Author: Maria Grazia Pia
33//
34// Creation date: 6 August 2001
35//
36// Modifications:
37//
38// -------------------------------------------------------------------
39
40#include "globals.hh"
41#include "G4ios.hh"
42#include <fstream>
43#include <iomanip>
44
45#include "G4ParticleDefinition.hh"
46#include "G4ParticleTypes.hh"
47//#include "G4ParticleTable.hh"
48#include "G4ParticleMomentum.hh"
49#include "G4DynamicParticle.hh"
50#include "G4ThreeVector.hh"
51#include "G4Track.hh"
52#include "G4SystemOfUnits.hh"
53
54#include "G4Material.hh"
55#include "G4ProcessManager.hh"
56#include "G4VParticleChange.hh"
57#include "G4ParticleChange.hh"
58#include "G4PVPlacement.hh"
59#include "G4Step.hh"
60#include "G4GRSVolume.hh"
61#include "G4Box.hh"
62#include "G4PVPlacement.hh"
63
64#include "G4CrossSectionElasticScreenedRutherford.hh"
65#include "G4FinalStateElasticScreenedRutherford.hh"
66#include "G4FinalStateElasticBrennerZaider.hh"
67
68#include "G4CrossSectionExcitationEmfietzoglou.hh"
69#include "G4FinalStateExcitationEmfietzoglou.hh"
70
71#include "G4CrossSectionExcitationBorn.hh"
72#include "G4FinalStateExcitationBorn.hh"
73
74#include "G4FinalStateProduct.hh"
75#include "G4DummyFinalState.hh"
76#include "G4DNAProcess.hh"
77
78//typedef G4DNAProcess<G4eCrossSectionScreenedRutherford,G4DummyFinalState> G4MyProcess;
79
80//typedef G4DNAProcess<G4CrossSectionElasticScreenedRutherford,G4FinalStateElasticScreenedRutherford> G4MyProcess;
81
82//typedef G4DNAProcess<G4CrossSectionElasticScreenedRutherford,G4FinalStateElasticBrennerZaider> G4MyProcess;
83
84//typedef G4DNAProcess<G4CrossSectionExcitationEmfietzoglou,G4FinalStateExcitationEmfietzoglou> G4MyProcess;
85
86typedef G4DNAProcess<G4CrossSectionExcitationBorn,G4FinalStateExcitationBorn> G4MyProcess;
87
88int main()
89{
90 // G4cout.setf( ios::scientific, ios::floatfield );
91
92 G4MyProcess* process = new G4MyProcess;
93
94
95 // Particle definitions
96 G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
97 G4ParticleDefinition* proton = G4Proton::ProtonDefinition();
98
99 // Create a DynamicParticle
100
101 G4double initX = 0.;
102 G4double initY = 0.;
103 G4double initZ = 1.;
104
105 G4ParticleMomentum direction(initX,initY,initZ);
106
107
108 G4cout << "Enter energy in keV" << G4endl;
109 G4double energy;
110 G4cin >> energy;
111 energy = energy * keV;
112
113 // G4DynamicParticle dynamicParticle(electron,direction,energy);
114
115 G4DynamicParticle dynamicParticle(proton,direction,energy);
116
117 // dynamicParticle.DumpInfo(0);
118
119
120// Materials
121 G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
122 G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
123 G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2);
124 water->AddElement(H,2);
125 water->AddElement(O,1);
126
127 // Dump the material table
128 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
129 G4int nMaterials = G4Material::GetNumberOfMaterials();
130 G4cout << "Available materials are: " << G4endl;
131 for (G4int mat = 0; mat < nMaterials; mat++)
132 {
133 G4cout << mat << ") "
134 << (*materialTable)[mat]->GetName()
135 << G4endl;
136 }
137
138 G4double dimX = 1 * mm;
139 G4double dimY = 1 * mm;
140 G4double dimZ = 1 * mm;
141
142 // Geometry
143
144 G4Box* theFrame = new G4Box ("Frame",dimX, dimY, dimZ);
145 G4LogicalVolume* logicalFrame = new G4LogicalVolume(theFrame,
146 water,
147 "LFrame", 0, 0, 0);
148 logicalFrame->SetMaterial(water);
149 G4PVPlacement* physicalFrame = new G4PVPlacement(0,G4ThreeVector(),
150 "PFrame",logicalFrame,0,false,0);
151
152
153 // Track
154 G4ThreeVector position(0.,0.,0.);
155 G4double time = 0. ;
156
157 G4Track* track = new G4Track(&dynamicParticle,time,position);
158 // Do I really need this?
159 G4GRSVolume* touche = new G4GRSVolume(physicalFrame, 0, position);
160 // track->SetTouchable(touche);
161
162 G4Step* step = new G4Step();
163 step->SetTrack(track);
164
165 G4StepPoint* point = new G4StepPoint();
166 point->SetPosition(position);
167 point->SetMaterial(water);
168 G4double safety = 10000.*cm;
169 point->SetSafety(safety);
170 step->SetPreStepPoint(point);
171
172 G4StepPoint* newPoint = new G4StepPoint();
173 G4ThreeVector newPosition(0.,0.,0.05*mm);
174 newPoint->SetPosition(newPosition);
175 newPoint->SetMaterial(water);
176 step->SetPostStepPoint(newPoint);
177
178 step->SetStepLength(1*micrometer);
179
180 track->SetStep(step);
181
182
183 // Calculate mean free path and cross section
184
185 G4ForceCondition* force = new G4ForceCondition;
186
187 G4double mfp = process->DumpMeanFreePath(*track,0.1,force);
188 G4double cross = 0.;
189
190 if (mfp > 0.0) cross = 1. / mfp;
191
192 G4double atomicDensity = 3.34192e+19 / (1./ mm3);
193
194 cross = cross/atomicDensity;
195
196 G4cout << "MeanFreePath = " << mfp
197 << " - Inverse mean free path = " << 1./(mfp/um)
198 << " - Cross section = " << cross /(cm*cm) << " cm-2 "
199 << G4endl;
200
201 G4cout << "Before invoking PostStepDoIt " << G4endl;
202
203
204 G4int nTry = 5;
205
206 for (G4int iTry = 0; iTry < nTry; iTry++)
207 {
208 // Retrieve final state produced
209 G4VParticleChange* particleChange = process->PostStepDoIt(*track,*step);
210
211 G4int nSecondaries = particleChange->GetNumberOfSecondaries();
212
213 G4cout << iTry << ") Number of secondary particles produced = " << nSecondaries << G4endl;
214
215 for (G4int i = 0; i < nSecondaries; i++)
216 {
217 G4Track* finalParticle = particleChange->GetSecondary(i) ;
218
219 G4double e = finalParticle->GetTotalEnergy();
220 G4double eKin = finalParticle->GetKineticEnergy();
221 G4double px = (finalParticle->GetMomentum()).x();
222 G4double py = (finalParticle->GetMomentum()).y();
223 G4double pz = (finalParticle->GetMomentum()).z();
224 // G4double theta = (finalParticle->GetMomentum()).theta();
225 // G4double phi = (finalParticle->GetMomentum()).phi();
226 // G4double p = std::sqrt(px*px+py*py+pz*pz);
227 G4String particleName = finalParticle->GetDefinition()->GetParticleName();
228 G4cout << "==== Final "
229 << particleName << " "
230 << "energy: " << e/keV << " keV, "
231 << "eKin: " << eKin/keV << " keV, "
232 << "(px,py,pz): ("
233 << px/keV << ","
234 << py/keV << ","
235 << pz/keV << ") keV "
236 << G4endl;
237 }
238
239 G4ParticleChange* pChange = dynamic_cast<G4ParticleChange*>(particleChange);
240
241 G4double eFinal = pChange->GetEnergy();
242 const G4ThreeVector* dirFinal = pChange->GetMomentumDirection();
243 G4double pX = dirFinal->x();
244 G4double pY = dirFinal->y();
245 G4double pZ = dirFinal->z();
246
247 G4cout << "==== Final energy: " << eFinal/keV << " keV, "
248 << "Modified direction (px,py,pz): ("
249 << pX << ","
250 << pY << ","
251 << pZ << ") keV "
252 << G4endl;
253 // delete particleChange;
254 }
255
256 // delete dirFinal;
257 // delete track;
258 // delete step;
259 // delete process;
260
261
262 // G4cout << "END OF THE MAIN PROGRAM" << G4endl;
263}
264
265
266
267
268
269
270
271
Note: See TracBrowser for help on using the repository browser.