source: trunk/source/processes/electromagnetic/utils/src/G4LossTableBuilder.cc@ 891

Last change on this file since 891 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 7.3 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// $Id: G4LossTableBuilder.cc,v 1.24 2007/02/16 11:59:35 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4LossTableBuilder
35//
36// Author: Vladimir Ivanchenko
37//
38// Creation date: 03.01.2002
39//
40// Modifications:
41//
42// 23-01-03 V.Ivanchenko Cut per region
43// 21-07-04 V.Ivanchenko Fix problem of range for dedx=0
44// 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko)
45// 07-12-04 Fix of BuildDEDX table (V.Ivanchenko)
46// 27-03-06 Add bool options isIonisation (V.Ivanchenko)
47// 16-01-07 Fill new (not old) DEDX table (V.Ivanchenko)
48// 12-02-07 Use G4LPhysicsFreeVector for the inverse range table (V.Ivanchenko)
49//
50// Class Description:
51//
52// -------------------------------------------------------------------
53//
54//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56
57#include "G4LossTableBuilder.hh"
58#include "G4PhysicsTable.hh"
59#include "G4PhysicsLogVector.hh"
60#include "G4PhysicsTableHelper.hh"
61#include "G4LPhysicsFreeVector.hh"
62
63//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64
65G4LossTableBuilder::G4LossTableBuilder()
66{}
67
68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69
70G4LossTableBuilder::~G4LossTableBuilder()
71{}
72
73//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74
75void G4LossTableBuilder::BuildDEDXTable(G4PhysicsTable* dedxTable,
76 const std::vector<G4PhysicsTable*>& list)
77{
78 size_t n_processes = list.size();
79 if(1 >= n_processes) return;
80
81 size_t n_vectors = (list[0])->length();
82 if(0 >= n_vectors) return;
83
84 G4bool b;
85
86 G4PhysicsVector* pv = (*(list[0]))[0];
87 size_t nbins = pv->GetVectorLength();
88 G4double elow = pv->GetLowEdgeEnergy(0);
89 G4double ehigh = pv->GetLowEdgeEnergy(nbins);
90 for (size_t i=0; i<n_vectors; i++) {
91
92 pv = new G4PhysicsLogVector(elow, ehigh, nbins);
93 for (size_t j=0; j<nbins; j++) {
94 G4double dedx = 0.0;
95 G4double energy = pv->GetLowEdgeEnergy(j);
96
97 for (size_t k=0; k<n_processes; k++) {
98 dedx += ((*(list[k]))[i])->GetValue(energy, b);
99 }
100 pv->PutValue(j, dedx);
101 G4PhysicsTableHelper::SetPhysicsVector(dedxTable, i, pv);
102 }
103 }
104}
105
106//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107
108void G4LossTableBuilder::BuildRangeTable(const G4PhysicsTable* dedxTable,
109 G4PhysicsTable* rangeTable,
110 G4bool isIonisation)
111// Build range table from the energy loss table
112{
113 size_t n_vectors = dedxTable->length();
114 if(!n_vectors) return;
115
116 G4bool b;
117 size_t n = 100;
118 G4double del = 1.0/(G4double)n;
119
120 for (size_t i=0; i<n_vectors; i++) {
121
122 if (rangeTable->GetFlag(i) || !isIonisation) {
123 G4PhysicsVector* pv = (*dedxTable)[i];
124 size_t nbins = pv->GetVectorLength();
125 size_t bin0 = 0;
126 G4double elow = pv->GetLowEdgeEnergy(0);
127 G4double ehigh = pv->GetLowEdgeEnergy(nbins);
128 G4double dedx1 = pv->GetValue(elow, b);
129
130 if(dedx1 == 0.0) {
131 for (size_t k=1; k<nbins; k++) {
132 bin0++;
133 elow = pv->GetLowEdgeEnergy(k);
134 dedx1 = pv->GetValue(elow, b);
135 if(dedx1 > 0.0) break;
136 }
137 nbins -= bin0;
138 }
139
140 G4PhysicsLogVector* v = new G4PhysicsLogVector(elow, ehigh, nbins);
141
142 G4double range = 2.*elow/dedx1;
143 //G4double range = elow/dedx1;
144 v->PutValue(0,range);
145 G4double energy1 = elow;
146
147 for (size_t j=1; j<nbins; j++) {
148
149 G4double energy2 = pv->GetLowEdgeEnergy(j+bin0);
150 G4double dedx2 = pv->GetValue(energy2, b);
151 G4double de = (energy2 - energy1) * del;
152 G4double energy = energy1 - de*0.5;
153
154 G4bool yes = true;
155 //G4bool yes = false;
156 if(dedx1 < DBL_MIN || dedx2 < DBL_MIN) yes = false;
157
158 G4double fac, f;
159
160 if(yes) fac = std::log(dedx2/dedx1)/std::log(energy2/energy1);
161 else fac = (dedx2 - dedx1)/(energy2 - energy1);
162
163 for (size_t k=0; k<n; k++) {
164 energy += de;
165 if(yes) f = dedx1*std::exp(fac*std::log(energy/energy1));
166 else f = dedx1 + fac*(energy - energy1);
167 if(f > DBL_MIN) range += de/f;
168 }
169 // G4cout << "Range i= " <<i << " j= " << j << G4endl;
170 v->PutValue(j,range);
171 energy1 = energy2;
172 dedx1 = dedx2;
173 }
174 G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
175 }
176 }
177}
178
179//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
180
181void G4LossTableBuilder::BuildInverseRangeTable(const G4PhysicsTable* rangeTable,
182 G4PhysicsTable* invRangeTable,
183 G4bool isIonisation)
184// Build inverse range table from the energy loss table
185{
186 size_t n_vectors = rangeTable->length();
187 if(!n_vectors) return;
188 G4bool b;
189
190 for (size_t i=0; i<n_vectors; i++) {
191
192 if (invRangeTable->GetFlag(i) || !isIonisation) {
193 G4PhysicsVector* pv = (*rangeTable)[i];
194 size_t nbins = pv->GetVectorLength();
195 G4double elow = pv->GetLowEdgeEnergy(0);
196 G4double ehigh = pv->GetLowEdgeEnergy(nbins-1);
197 G4double rlow = pv->GetValue(elow, b);
198 G4double rhigh = pv->GetValue(ehigh, b);
199
200 rhigh *= std::exp(std::log(rhigh/rlow)/((G4double)(nbins-1)));
201
202
203 G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(nbins,rlow,rhigh);
204 for (size_t j=0; j<nbins; j++) {
205 G4double e = pv->GetLowEdgeEnergy(j);
206 G4double r = pv->GetValue(e, b);
207 v->PutValues(j,r,e);
208 }
209
210 G4PhysicsTableHelper::SetPhysicsVector(invRangeTable, i, v);
211 }
212 }
213}
214
215//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
216
Note: See TracBrowser for help on using the repository browser.