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

Last change on this file since 1316 was 1316, checked in by garnier, 14 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.