| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // -------------------------------------------
27 //
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4LossTableBuilder
32 //
33 // Author: Vladimir Ivanchenko
34 //
35 // Creation date: 03.01.2002
36 //
37 // Modifications:
38 //
39 // 23-01-03 V.Ivanchenko Cut per region
40 // 21-07-04 V.Ivanchenko Fix problem of range
41 // 08-11-04 Migration to new interface of Stor
42 // 07-12-04 Fix of BuildDEDX table (V.Ivanchen
43 // 27-03-06 Add bool options isIonisation (V.I
44 // 16-01-07 Fill new (not old) DEDX table (V.I
45 // 12-02-07 Use G4LPhysicsFreeVector for the i
46 // 24-06-09 Removed hidden bin in G4PhysicsVec
47 //
48 // Class Description:
49 //
50 // -------------------------------------------
51 //
52 //....oooOO0OOooo........oooOO0OOooo........oo
53 //....oooOO0OOooo........oooOO0OOooo........oo
54
55 #include "G4LossTableBuilder.hh"
56 #include "G4SystemOfUnits.hh"
57 #include "G4PhysicsTable.hh"
58 #include "G4PhysicsLogVector.hh"
59 #include "G4PhysicsTableHelper.hh"
60 #include "G4PhysicsFreeVector.hh"
61 #include "G4ProductionCutsTable.hh"
62 #include "G4MaterialCutsCouple.hh"
63 #include "G4Material.hh"
64 #include "G4VEmModel.hh"
65 #include "G4ParticleDefinition.hh"
66 #include "G4LossTableManager.hh"
67 #include "G4EmParameters.hh"
68
69 G4bool G4LossTableBuilder::baseMatFlag = false
70 std::vector<G4double>* G4LossTableBuilder::the
71 std::vector<G4int>* G4LossTableBuilder::the
72 std::vector<G4bool>* G4LossTableBuilder::the
73
74 //....oooOO0OOooo........oooOO0OOooo........oo
75
76 G4LossTableBuilder::G4LossTableBuilder(G4bool
77 : isInitializer(master)
78 {
79 theParameters = G4EmParameters::Instance();
80 if (nullptr == theFlag) {
81 theDensityFactor = new std::vector<G4doubl
82 theDensityIdx = new std::vector<G4int>;
83 theFlag = new std::vector<G4bool>;
84 }
85 }
86
87 //....oooOO0OOooo........oooOO0OOooo........oo
88
89 G4LossTableBuilder::~G4LossTableBuilder()
90 {
91 if (isInitializer) {
92 delete theDensityFactor;
93 delete theDensityIdx;
94 delete theFlag;
95 theDensityFactor = nullptr;
96 theDensityIdx = nullptr;
97 theFlag = nullptr;
98 }
99 }
100
101 //....oooOO0OOooo........oooOO0OOooo........oo
102
103 const std::vector<G4int>* G4LossTableBuilder::
104 {
105 return theDensityIdx;
106 }
107
108 //....oooOO0OOooo........oooOO0OOooo........oo
109
110 const std::vector<G4double>* G4LossTableBuilde
111 {
112 return theDensityFactor;
113 }
114
115 //....oooOO0OOooo........oooOO0OOooo........oo
116
117 G4bool G4LossTableBuilder::GetFlag(std::size_t
118 {
119 return (idx < theFlag->size()) ? (*theFlag)[
120 }
121
122 //....oooOO0OOooo........oooOO0OOooo........oo
123
124 G4bool G4LossTableBuilder::GetBaseMaterialFlag
125 {
126 return baseMatFlag;
127 }
128
129 //....oooOO0OOooo........oooOO0OOooo........oo
130
131 void
132 G4LossTableBuilder::BuildDEDXTable(G4PhysicsTa
133 const std::
134 {
135 InitialiseBaseMaterials(dedxTable);
136 std::size_t n_processes = list.size();
137 if(1 >= n_processes) { return; }
138
139 std::size_t nCouples = dedxTable->size();
140 //G4cout << "Nproc= " << n_processes << " nC
141 // << dedxTable->size() << G4endl;
142 if(0 >= nCouples) { return; }
143
144 for (std::size_t i=0; i<nCouples; ++i) {
145 auto pv0 = static_cast<G4PhysicsLogVector*
146 //if (0 == i) G4cout << i << ". pv0=" << p
147 if(pv0 == nullptr) { continue; }
148 std::size_t npoints = pv0->GetVectorLength
149 auto pv = new G4PhysicsLogVector(*pv0);
150 for (std::size_t j=0; j<npoints; ++j) {
151 G4double dedx = 0.0;
152 for (std::size_t k=0; k<n_processes; ++k
153 const G4PhysicsVector* pv1 = (*(list[k]))[i]
154 //if (0 == i) G4cout << " " << k << ". p
155 dedx += (*pv1)[j];
156 }
157 pv->PutValue(j, dedx);
158 }
159 if(splineFlag) { pv->FillSecondDerivatives
160 G4PhysicsTableHelper::SetPhysicsVector(ded
161 }
162 //G4cout << "### G4LossTableBuilder::BuildDE
163 //G4cout << *dedxTable << G4endl;
164 }
165
166 //....oooOO0OOooo........oooOO0OOooo........oo
167
168 void G4LossTableBuilder::BuildRangeTable(const
169 G4Phy
170 // Build range table from the energy loss tabl
171 {
172 //G4cout << "### G4LossTableBuilder::BuildRa
173 //G4cout << *const_cast<G4PhysicsTable*>(ded
174 const std::size_t nCouples = dedxTable->size
175 if(0 >= nCouples) { return; }
176
177 const std::size_t n = 100;
178 const G4double del = 1.0/(G4double)n;
179
180 for (std::size_t i=0; i<nCouples; ++i) {
181 auto pv = static_cast<G4PhysicsLogVector*>
182 if((pv == nullptr) || (isBaseMatActive &&
183 std::size_t npoints = pv->GetVectorLength(
184 std::size_t bin0 = 0;
185 G4double elow = pv->Energy(0);
186 G4double ehigh = pv->Energy(npoints-1);
187 G4double dedx1 = (*pv)[0];
188
189 // protection for specific cases dedx=0
190 if(dedx1 == 0.0) {
191 for (std::size_t k=1; k<npoints; ++k) {
192 ++bin0;
193 elow = pv->Energy(k);
194 dedx1 = (*pv)[k];
195 if(dedx1 > 0.0) { break; }
196 }
197 npoints -= bin0;
198 }
199
200 // initialisation of a new vector
201 if(npoints < 3) { npoints = 3; }
202
203 delete (*rangeTable)[i];
204 G4PhysicsLogVector* v;
205 if(0 == bin0) { v = new G4PhysicsLogVector
206 else { v = new G4PhysicsLogVector(elow, eh
207
208 // assumed dedx proportional to beta
209 G4double energy1 = v->Energy(0);
210 G4double range = 2.*energy1/dedx1;
211 /*
212 G4cout << "New Range vector Npoints=" << v
213 << " coupleIdx=" << i << " spline=" << v-
214 << " Elow=" << v->GetMinEnergy() <<" Ehig
215 << " DEDX(Elow)=" << dedx1 << " R(Elow)="
216 */
217 v->PutValue(0,range);
218
219 for (std::size_t j=1; j<npoints; ++j) {
220
221 G4double energy2 = v->Energy(j);
222 G4double de = (energy2 - energy1) *
223 G4double energy = energy2 + de*0.5;
224 G4double sum = 0.0;
225 std::size_t idx = j - 1;
226 for (std::size_t k=0; k<n; ++k) {
227 energy -= de;
228 dedx1 = pv->Value(energy, idx);
229 if(dedx1 > 0.0) { sum += de/dedx1; }
230 }
231 range += sum;
232 /*
233 if(energy < 10.)
234 G4cout << "j= " << j << " e1= " << energy1 <
235 << " n= " << n << " range=" << range<
236 */
237 v->PutValue(j,range);
238 energy1 = energy2;
239 }
240 if(splineFlag) { v->FillSecondDerivatives(
241 G4PhysicsTableHelper::SetPhysicsVector(ran
242 }
243 //G4cout << "### Range table" << G4endl;
244 //G4cout << *rangeTable << G4endl;
245 }
246
247 //....oooOO0OOooo........oooOO0OOooo........oo
248
249 void
250 G4LossTableBuilder::BuildInverseRangeTable(con
251 G4P
252 // Build inverse range table from the energy l
253 {
254 std::size_t nCouples = rangeTable->size();
255 if(0 >= nCouples) { return; }
256
257 for (std::size_t i=0; i<nCouples; ++i) {
258 G4PhysicsVector* pv = (*rangeTable)[i];
259 if((pv == nullptr) || (isBaseMatActive &&
260 std::size_t npoints = pv->GetVectorLength(
261
262 delete (*invRangeTable)[i];
263 auto v = new G4PhysicsFreeVector(npoints,
264
265 for (std::size_t j=0; j<npoints; ++j) {
266 G4double e = pv->Energy(j);
267 G4double r = (*pv)[j];
268 v->PutValues(j,r,e);
269 }
270 if (splineFlag) { v->FillSecondDerivatives
271 v->EnableLogBinSearch(theParameters->Numbe
272
273 G4PhysicsTableHelper::SetPhysicsVector(inv
274 }
275 //G4cout << "### Inverse range table" << G4e
276 //G4cout << *invRangeTable << G4endl;
277 }
278
279 //....oooOO0OOooo........oooOO0OOooo........oo
280
281 void G4LossTableBuilder::InitialiseBaseMateria
282 {
283 if(!isInitializer) { return; }
284 const G4ProductionCutsTable* theCoupleTable=
285 G4ProductionCutsTable::GetProductionCutsTa
286 std::size_t nCouples = theCoupleTable->GetTa
287 std::size_t nFlags = theFlag->size();
288 /*
289 G4cout << "### InitialiseBaseMaterials: nCou
290 << " nFlags=" << nFlags << " isInit:" << is
291 << " baseMat:" << baseMatFlag << G4endl;
292 */
293 // define base material flag
294 if(isBaseMatActive && !baseMatFlag) {
295 for(G4int i=0; i<(G4int)nCouples; ++i) {
296 if(nullptr != theCoupleTable->GetMateria
297 baseMatFlag = true;
298 isInitialized = false;
299 break;
300 }
301 }
302 }
303
304 if(nFlags != nCouples) { isInitialized = fal
305 if(isInitialized) { return; }
306
307 // reserve memory
308 theFlag->resize(nCouples, true);
309 theDensityFactor->resize(nCouples,1.0);
310 theDensityIdx->resize(nCouples, 0);
311
312 // define default flag and index of used mat
313 for (G4int i=0; i<(G4int)nCouples; ++i) {
314 (*theFlag)[i] = (nullptr == table) ? true
315 (*theDensityIdx)[i] = i;
316 }
317 isInitialized = true;
318 if (!baseMatFlag) { return; }
319
320 // use base materials
321 for (G4int i=0; i<(G4int)nCouples; ++i) {
322 // base material is needed only for a coup
323 // initialised and for which tables will b
324 auto couple = theCoupleTable->GetMaterialC
325 auto pcuts = couple->GetProductionCuts();
326 auto mat = couple->GetMaterial();
327 auto bmat = mat->GetBaseMaterial();
328
329 // base material exists - find it and chec
330 if(nullptr != bmat) {
331 for(G4int j=0; j<(G4int)nCouples; ++j) {
332 if(j == i) { continue; }
333 auto bcouple = theCoupleTable->GetMaterialCu
334
335 if(bcouple->GetMaterial() == bmat &&
336 bcouple->GetProductionCuts() == pcuts) {
337
338 // based couple exist in the same region
339 (*theDensityFactor)[i] = mat->GetDensity()
340 (*theDensityIdx)[i] = j;
341 (*theFlag)[i] = false;
342
343 // ensure that there will no double initia
344 (*theDensityFactor)[j] = 1.0;
345 (*theDensityIdx)[j] = j;
346 (*theFlag)[j] = true;
347 break;
348 }
349 }
350 }
351 }
352 }
353
354 //....oooOO0OOooo........oooOO0OOooo........oo
355
356 G4PhysicsTable*
357 G4LossTableBuilder::BuildTableForModel(G4Physi
358 G4VEmMo
359 const G
360 G4doubl
361 G4bool
362 {
363 // check input
364 G4PhysicsTable* table = G4PhysicsTableHelper
365 if (nullptr == table) { return table; }
366 if (aTable != nullptr && aTable != table) {
367 aTable->clearAndDestroy();
368 delete aTable;
369 }
370
371 InitialiseBaseMaterials(table);
372 G4int nbins = theParameters->NumberOfBinsPer
373
374 // Access to materials
375 const G4ProductionCutsTable* theCoupleTable=
376 G4ProductionCutsTable::GetProductionCu
377 std::size_t numOfCouples = theCoupleTable->G
378 /*
379 G4cout << " G4LossTableBuilder::BuildTable
380 << " isMaster=" << isInitializer << " model
381 << " " << part->GetParticleName() << G4end
382 */
383 G4PhysicsLogVector* aVector = nullptr;
384
385 for(G4int i=0; i<(G4int)numOfCouples; ++i) {
386 //G4cout << i << ". " << (*theFlag)[i] <<
387 if (table->GetFlag(i)) {
388
389 // create physics vector and fill it
390 auto couple = theCoupleTable->GetMateria
391 delete (*table)[i];
392
393 // if start from zero then change the sc
394 const G4Material* mat = couple->GetMater
395
396 G4double tmin = std::max(emin, model->Mi
397 if(0.0 >= tmin) { tmin = CLHEP::eV; }
398 G4int n = nbins;
399
400 if(tmin >= emax) {
401 aVector = nullptr;
402 } else {
403 n *= G4lrint(std::log10(emax/tmin));
404 n = std::max(n, 3);
405 aVector = new G4PhysicsLogVector(tmin,
406 }
407
408 if(nullptr != aVector) {
409 //G4cout << part->GetParticleName() <<
410 // << " emin= " << tmin << " emax=" << e
411 for(G4int j=0; j<=n; ++j) {
412 G4double e = aVector->Energy(j);
413 G4double y = model->Value(couple, part, e)
414 //G4cout << " " << j << ") E=" << e <
415 aVector->PutValue(j, y);
416 }
417 if(spline) { aVector->FillSecondDeriva
418 }
419 G4PhysicsTableHelper::SetPhysicsVector(t
420 }
421 }
422 /*
423 G4cout << "G4LossTableBuilder::BuildTableFor
424 << part->GetParticleName() << " and "
425 << " " << table << G4endl;
426 */
427 //G4cout << *table << G4endl;
428 return table;
429 }
430
431 //....oooOO0OOooo........oooOO0OOooo........oo
432