source: trunk/source/processes/electromagnetic/standard/test/G4ICRU73QOModelTest.cc@ 1330

Last change on this file since 1330 was 1316, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 10.0 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// ------------------------------------------------------------
28//
29// To print eloss by using of ICRU73QOModel
30//
31#include "G4Material.hh"
32#include "G4NistManager.hh"
33
34#include "G4BetheBlochModel.hh"
35#include "G4BraggModel.hh"
36#include "G4ICRU73QOModel.hh"
37#include "G4VEmModel.hh"
38#include "G4LossTableManager.hh"
39#include "G4ParticleTable.hh"
40#include "G4hIonisation.hh"
41#include "G4PhysicsLogVector.hh"
42
43#include "globals.hh"
44#include "G4UnitsTable.hh"
45
46#include "G4Proton.hh"
47#include "G4AntiProton.hh"
48#include "G4MuonPlus.hh"
49#include "G4MuonMinus.hh"
50
51#include <vector>
52
53int main() {
54
55 G4UnitDefinition::BuildUnitsTable();
56
57 G4NistManager* nist = G4NistManager::Instance();
58 G4Material* material = nist->FindOrBuildMaterial("G4_H");
59
60 G4cout << *(G4Material::GetMaterialTable()) << G4endl;
61 const std::vector<G4String> mnames = nist->GetNistMaterialNames();
62
63 G4ParticleTable* partTable = G4ParticleTable::GetParticleTable();
64
65 // initialise proton processes (models)
66 //
67 G4ParticleDefinition* part = G4AntiProton::AntiProton();
68 G4ParticleDefinition* prot = G4Proton::Proton();
69 //G4ParticleDefinition* part = G4MuonMinus::MuonMinus();
70
71 G4cout << "pName = " << part->GetParticleName() << G4endl;
72 G4String pName = part->GetParticleName();
73
74 partTable->SetReadiness();
75
76 G4VEmModel* bragg = new G4BraggModel();
77 G4VEmModel* bethe1 = new G4BetheBlochModel();
78 G4VEmModel* bethe2 = new G4BetheBlochModel();
79 G4VEmModel* qo = new G4ICRU73QOModel();
80
81 G4DataVector v;
82 v.resize(100);
83
84 bragg->Initialise(prot,v);
85 bethe1->Initialise(prot,v);
86 bethe2->Initialise(part,v);
87 qo->Initialise(part,v);
88
89 G4double Emin = 1.1*keV;
90 G4double Emax = 100.01*MeV;
91 G4double Ecut = 1*GeV;
92
93 size_t nBinTab = 35;
94
95 G4PhysicsLogVector* pVector;
96 pVector = new G4PhysicsLogVector(Emin, Emax, nBinTab);
97
98 G4double Energy = 0.;
99 G4double ELos0 = 0.;
100 G4double ELos1 = 0.;
101 G4double ELos2 = 0.;
102 G4double ELos3 = 0.;
103
104 G4String mat = material->GetName();
105 G4String asciiFileName = pName + "_QOEloss_diff_" + mat + ".out";
106
107 std::ofstream asciiFile(asciiFileName, std::ios::out);
108 if(asciiFile.is_open()) {
109 asciiFile << " Energy(Mev)\t J(MeV/mm) (bragg bethe quantosc)" << G4endl;
110 } else {
111 G4cout << "ERROR file <" << asciiFileName << "> is not opened" << G4endl;
112 exit(1);
113 }
114 G4double fact1[98];
115 G4double fact2[98];
116 G4double fact3[98];
117
118 Energy = 2.*MeV;
119 part = G4AntiProton::AntiProton();
120 pName = part->GetParticleName();
121 G4cout << "### E(MeV)= " << Energy/MeV << " for " << pName
122 << std::setprecision(4)
123 << G4endl;
124 G4cout << " const G4double factorBethe[98] = { 1.0, " << G4endl;
125 for(G4int j=0; j<98; ++j) {
126 const G4Material* mat = nist->FindOrBuildMaterial(mnames[j]);
127 G4int Z = G4int(mat->GetZ());
128 ELos0 = bragg->ComputeDEDXPerVolume(mat,prot,Energy,Ecut);
129 ELos1 = bethe1->ComputeDEDXPerVolume(mat,prot,Energy,Ecut);
130 ELos2 = bethe2->ComputeDEDXPerVolume(mat,part,Energy,Ecut);
131 ELos3 = qo->ComputeDEDXPerVolume(mat,part,Energy,Ecut);
132 fact1[j]= ELos2/ELos3;
133 if(ELos2 > ELos0) { fact1[j] *= ELos0/ELos1; }
134 fact2[j]= ELos3*fact1[j]/ELos0;
135 fact3[j]= ELos1/ELos2;
136 G4cout << fact1[j] << ", ";
137 if((Z/10)*10 == Z) { G4cout << " // " << Z-9 << " - " << Z << G4endl; }
138 }
139 G4cout << " } " << G4endl;
140 G4cout << " ##### ICRU73(pbar)/Bragg(p) ratio at 2 MeV: " << G4endl;
141 for(G4int j=1; j<98; ++j) {
142 G4int Z = j;
143 G4cout << fact2[j] << ", ";
144 if((Z/10)*10 == Z) { G4cout << " // " << Z-9 << " - " << Z << G4endl; }
145 }
146 G4cout << " " << G4endl;
147 G4cout << " ##### p/pbar BetheBloch ratio at 2 MeV: " << G4endl;
148 for(G4int j=1; j<98; ++j) {
149 G4int Z = j;
150 G4cout << fact3[j] << ", ";
151 if((Z/10)*10 == Z) { G4cout << " // " << Z-9 << " - " << Z << G4endl; }
152 }
153 G4cout << " " << G4endl;
154 part = G4Proton::Proton();
155 pName = part->GetParticleName();
156 G4cout << "### E(MeV)= " << Energy/MeV << " for " << pName
157 << " Bethe/Bragg" <<G4endl;
158 G4cout << " const G4double factorBragg[98] = { 1.0, " << G4endl;
159 for(G4int j=0; j<98; ++j) {
160 const G4Material* mat = nist->FindOrBuildMaterial(mnames[j]);
161 G4int Z = G4int(mat->GetZ());
162 ELos1 = bragg->ComputeDEDXPerVolume(mat,prot,Energy,Ecut);
163 ELos2 = bethe1->ComputeDEDXPerVolume(mat,prot,Energy,Ecut);
164 fact2[j]= ELos2/ELos1;
165 G4cout << fact2[j] << ", ";
166 if((Z/10)*10 == Z) { G4cout << " // " << Z-9 << " - " << Z << G4endl; }
167 }
168 G4cout << " } " << G4endl;
169
170 Energy = 2.*MeV;
171 G4cout << "### E(MeV)= " << Energy/MeV << " for " << pName
172 << std::setprecision(5) << G4endl;
173 for(G4int j=0; j<98; ++j) {
174 const G4Material* mat = nist->FindOrBuildMaterial(mnames[j]);
175 G4double Z = mat->GetZ();
176 ELos1 = bragg->ComputeDEDXPerVolume(mat,prot,Energy,Ecut);
177 ELos2 = bethe1->ComputeDEDXPerVolume(mat,prot,Energy,Ecut);
178 ELos3 = qo->ComputeDEDXPerVolume(mat,part,Energy,Ecut);
179 G4cout << "Z= " << Z << " p bragg/bethe-1= " << ELos1/ELos2 -1
180 << " go(pbar)/bragg(p)-1= " << ELos3/ELos1 -1 << " "
181 << mnames[j] << " Bethe(MeV*cm^2/g)= "
182 << ELos2*g/(cm2*MeV*mat->GetDensity()) << G4endl;
183 }
184 Energy = 10.*MeV;
185 G4cout << "### E(MeV)= " << Energy/MeV << " for " << pName
186 << std::setprecision(5) << G4endl;
187 for(G4int j=0; j<98; ++j) {
188 const G4Material* mat = nist->FindOrBuildMaterial(mnames[j]);
189 G4double Z = mat->GetZ();
190 ELos1 = bragg->ComputeDEDXPerVolume(mat,prot,Energy,Ecut);
191 ELos2 = bethe1->ComputeDEDXPerVolume(mat,prot,Energy,Ecut);
192 ELos3 = qo->ComputeDEDXPerVolume(mat,part,Energy,Ecut);
193 G4cout << "Z= " << Z << " p bragg/bethe-1= " << ELos1/ELos2 -1
194 << " go(pbar)/bragg(p)-1= " << ELos3/ELos1 -1 << " "
195 << mnames[j] << " Bethe(MeV*cm^2/g)= "
196 << ELos2*g/(cm2*MeV*mat->GetDensity()) << G4endl;
197 }
198
199 part = G4AntiProton::AntiProton();
200 pName = part->GetParticleName();
201 Energy = 2.*MeV;
202 G4cout << "### E(MeV)= " << Energy/MeV << " for " << pName
203 << std::setprecision(5) << G4endl;
204 for(G4int j=0; j<98; ++j) {
205 const G4Material* mat = nist->FindOrBuildMaterial(mnames[j]);
206 G4double Z = mat->GetZ();
207 ELos1 = bragg->ComputeDEDXPerVolume(mat,prot,Energy,Ecut);
208 ELos2 = bethe2->ComputeDEDXPerVolume(mat,part,Energy,Ecut);
209 ELos3 = qo->ComputeDEDXPerVolume(mat,part,Energy,Ecut);
210 G4cout << "Z= " << Z << " bragg(p)/bethe(pbar)-1= " << ELos1/ELos2 -1
211 << " pbar go/bethe-1= " << ELos3/ELos2 -1 << " "
212 << mnames[j] << " Bethe(MeV*cm^2/g)= "
213 << ELos2*g/(cm2*MeV*mat->GetDensity()) << G4endl;
214 }
215 Energy = 10.*MeV;
216 G4cout << "### E(MeV)= " << Energy/MeV << " for " << pName
217 << std::setprecision(5) << G4endl;
218 for(G4int j=0; j<98; ++j) {
219 const G4Material* mat = nist->FindOrBuildMaterial(mnames[j]);
220 G4double Z = mat->GetZ();
221 ELos1 = bragg->ComputeDEDXPerVolume(mat,prot,Energy,Ecut);
222 ELos2 = bethe2->ComputeDEDXPerVolume(mat,part,Energy,Ecut);
223 ELos3 = qo->ComputeDEDXPerVolume(mat,part,Energy,Ecut);
224 G4cout << "Z= " << Z << " bragg(p)/bethe(pbar)-1= " << ELos1/ELos2 -1
225 << " pbar go/bethe-1= " << ELos3/ELos2 -1 << " "
226 << mnames[j] << " Bethe(MeV*cm^2/g)= "
227 << ELos2*g/(cm2*MeV*mat->GetDensity()) << G4endl;
228 }
229
230 // Write to file
231 for (size_t i = 0; i < nBinTab+1; i++) {
232 Energy = pVector->GetLowEdgeEnergy(i);
233 ELos1 = bragg->ComputeDEDXPerVolume(material,part,Energy,Ecut);
234 ELos2 = bethe2->ComputeDEDXPerVolume(material,part,Energy,Ecut);
235 ELos3 = qo->ComputeDEDXPerVolume(material,part,Energy,Ecut);
236 asciiFile << std::setiosflags(std::ios::fixed)
237 << std::setprecision(5)
238 << std::setiosflags(std::ios::right)
239 << std::setw(10);
240 asciiFile << Energy;
241 asciiFile << "\t";
242 asciiFile << std::setiosflags(std::ios::fixed)
243 << std::setprecision(5)
244 << std::setiosflags(std::ios::right)
245 << std::setw(10);
246 asciiFile << ELos1;
247 asciiFile << "\t";
248 asciiFile << std::setiosflags(std::ios::fixed)
249 << std::setprecision(5)
250 << std::setiosflags(std::ios::right)
251 << std::setw(10);
252 asciiFile << ELos2;
253 asciiFile << "\t";
254 asciiFile << std::setiosflags(std::ios::fixed)
255 << std::setprecision(5)
256 << std::setiosflags(std::ios::right)
257 << std::setw(10);
258 asciiFile << ELos3
259 << G4endl;
260 }
261
262 delete pVector;
263
264 return EXIT_SUCCESS;
265}
Note: See TracBrowser for help on using the repository browser.