ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LossTableBuilder.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LossTableBuilder.cc
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 // 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 for dedx=0
41 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko)
42 // 07-12-04 Fix of BuildDEDX table (V.Ivanchenko)
43 // 27-03-06 Add bool options isIonisation (V.Ivanchenko)
44 // 16-01-07 Fill new (not old) DEDX table (V.Ivanchenko)
45 // 12-02-07 Use G4LPhysicsFreeVector for the inverse range table (V.Ivanchenko)
46 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
47 //
48 // Class Description:
49 //
50 // -------------------------------------------------------------------
51 //
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54 
55 #include "G4LossTableBuilder.hh"
56 #include "G4SystemOfUnits.hh"
57 #include "G4PhysicsTable.hh"
58 #include "G4PhysicsLogVector.hh"
59 #include "G4PhysicsTableHelper.hh"
60 #include "G4LPhysicsFreeVector.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 std::vector<G4double>* G4LossTableBuilder::theDensityFactor = nullptr;
70 std::vector<G4int>* G4LossTableBuilder::theDensityIdx = nullptr;
71 std::vector<G4bool>* G4LossTableBuilder::theFlag = nullptr;
72 #ifdef G4MULTITHREADED
73  G4Mutex G4LossTableBuilder::ltbMutex = G4MUTEX_INITIALIZER;
74 #endif
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 
79 {
81  splineFlag = true;
82  isInitialized = false;
83  if(isMaster || !theFlag) {
84 #ifdef G4MULTITHREADED
85  G4MUTEXLOCK(&ltbMutex);
86  if(isMaster || !theFlag) {
87 #endif
88  isMaster = true;
89  theDensityFactor = new std::vector<G4double>;
90  theDensityIdx = new std::vector<G4int>;
91  theFlag = new std::vector<G4bool>;
92  } else {
93  isMaster = false;
94 #ifdef G4MULTITHREADED
95  }
96  G4MUTEXUNLOCK(&ltbMutex);
97 #endif
98  }
99 }
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102 
104 {
105  if(isMaster) {
106  delete theDensityFactor;
107  delete theDensityIdx;
108  delete theFlag;
109  theDensityFactor = nullptr;
110  theDensityIdx = nullptr;
111  theFlag = nullptr;
112  }
113 }
114 
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116 
117 const std::vector<G4int>* G4LossTableBuilder::GetCoupleIndexes() const
118 {
119  return theDensityIdx;
120 }
121 
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 
124 const std::vector<G4double>* G4LossTableBuilder::GetDensityFactors() const
125 {
126  return theDensityFactor;
127 }
128 
129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
130 
132 {
133  if(theFlag->empty()) { InitialiseBaseMaterials(); }
134  return (idx < theFlag->size()) ? (*theFlag)[idx] : false;
135 }
136 
137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138 
139 void
141  const std::vector<G4PhysicsTable*>& list)
142 {
143  InitialiseBaseMaterials(dedxTable);
144  size_t n_processes = list.size();
145  //G4cout << "Nproc= " << n_processes << " Ncoup= "
146  //<< dedxTable->size() << G4endl;
147  if(1 >= n_processes) { return; }
148 
149  size_t nCouples = dedxTable->size();
150  if(0 >= nCouples) { return; }
151 
152  for (size_t i=0; i<nCouples; ++i) {
153  G4PhysicsLogVector* pv0 = static_cast<G4PhysicsLogVector*>((*(list[0]))[i]);
154  if(pv0) {
155  size_t npoints = pv0->GetVectorLength();
156  G4PhysicsLogVector* pv = new G4PhysicsLogVector(*pv0);
157  pv->SetSpline(splineFlag);
158  for (size_t j=0; j<npoints; ++j) {
159  G4double dedx = 0.0;
160  for (size_t k=0; k<n_processes; ++k) {
161  G4PhysicsVector* pv1 = (*(list[k]))[i];
162  dedx += (*pv1)[j];
163  }
164  pv->PutValue(j, dedx);
165  }
166  if(splineFlag) { pv->FillSecondDerivatives(); }
167  G4PhysicsTableHelper::SetPhysicsVector(dedxTable, i, pv);
168  }
169  }
170 }
171 
172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
173 
175  G4PhysicsTable* rangeTable,
176  G4bool useBM)
177 // Build range table from the energy loss table
178 {
179  size_t nCouples = dedxTable->size();
180  if(0 >= nCouples) { return; }
181 
182  size_t n = 100;
183  G4double del = 1.0/(G4double)n;
184 
185  for (size_t i=0; i<nCouples; ++i) {
186  G4PhysicsLogVector* pv = static_cast<G4PhysicsLogVector*>((*dedxTable)[i]);
187  if((pv == nullptr) || (useBM && !(*theFlag)[i])) { continue; }
188  size_t npoints = pv->GetVectorLength();
189  size_t bin0 = 0;
190  G4double elow = pv->Energy(0);
191  G4double ehigh = pv->Energy(npoints-1);
192  G4double dedx1 = (*pv)[0];
193 
194  //G4cout << "i= " << i << "npoints= " << npoints << " dedx1= "
195  //<< dedx1 << G4endl;
196 
197  // protection for specific cases dedx=0
198  if(dedx1 == 0.0) {
199  for (size_t k=1; k<npoints; ++k) {
200  ++bin0;
201  elow = pv->Energy(k);
202  dedx1 = (*pv)[k];
203  if(dedx1 > 0.0) { break; }
204  }
205  npoints -= bin0;
206  }
207  //G4cout<<"New Range vector" << G4endl;
208  //G4cout<<"nbins= "<<npoints-1<<" elow= "<<elow<<" ehigh= "<<ehigh
209  // <<" bin0= " << bin0 <<G4endl;
210 
211  // initialisation of a new vector
212  if(npoints < 2) { npoints = 2; }
213 
214  delete (*rangeTable)[i];
216  if(0 == bin0) { v = new G4PhysicsLogVector(*pv); }
217  else { v = new G4PhysicsLogVector(elow, ehigh, npoints-1); }
218 
219  // dedx is exact zero cannot build range table
220  if(2 == npoints) {
221  v->PutValue(0,1000.);
222  v->PutValue(1,2000.);
223  G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
224  return;
225  }
226  v->SetSpline(splineFlag);
227 
228  // assumed dedx proportional to beta
229  G4double energy1 = v->Energy(0);
230  G4double range = 2.*energy1/dedx1;
231  //G4cout << "range0= " << range << G4endl;
232  v->PutValue(0,range);
233 
234  for (size_t j=1; j<npoints; ++j) {
235 
236  G4double energy2 = v->Energy(j);
237  G4double de = (energy2 - energy1) * del;
238  G4double energy = energy2 + de*0.5;
239  G4double sum = 0.0;
240  //G4cout << "j= " << j << " e1= " << energy1 << " e2= " << energy2
241  // << " n= " << n << G4endl;
242  for (size_t k=0; k<n; ++k) {
243  energy -= de;
244  dedx1 = pv->Value(energy);
245  if(dedx1 > 0.0) { sum += de/dedx1; }
246  }
247  range += sum;
248  v->PutValue(j,range);
249  energy1 = energy2;
250  }
251  if(splineFlag) { v->FillSecondDerivatives(); }
252  G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
253  }
254 }
255 
256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
257 
258 void
260  G4PhysicsTable* invRangeTable,
261  G4bool useBM)
262 // Build inverse range table from the energy loss table
263 {
264  size_t nCouples = rangeTable->size();
265  if(0 >= nCouples) { return; }
266 
267  for (size_t i=0; i<nCouples; ++i) {
268  G4PhysicsVector* pv = (*rangeTable)[i];
269  if((pv == nullptr) || (useBM && !(*theFlag)[i])) { continue; }
270  size_t npoints = pv->GetVectorLength();
271  G4double rlow = (*pv)[0];
272  G4double rhigh = (*pv)[npoints-1];
273 
274  delete (*invRangeTable)[i];
275  G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(npoints,rlow,rhigh);
276  v->SetSpline(splineFlag);
277 
278  for (size_t j=0; j<npoints; ++j) {
279  G4double e = pv->Energy(j);
280  G4double r = (*pv)[j];
281  v->PutValues(j,r,e);
282  }
283  if(splineFlag) { v->FillSecondDerivatives(); }
284 
285  G4PhysicsTableHelper::SetPhysicsVector(invRangeTable, i, v);
286  }
287 }
288 
289 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
290 
292 {
293  if(!isMaster) { return; }
294  const G4ProductionCutsTable* theCoupleTable=
296  size_t nCouples = theCoupleTable->GetTableSize();
297  size_t nFlags = theFlag->size();
298  if(isInitialized && nFlags == nCouples) { return; }
299 
300  isInitialized = true;
301  if(0 == nFlags) {
302  theDensityFactor->reserve(nCouples);
303  theDensityIdx->reserve(nCouples);
304  theFlag->reserve(nCouples);
305  }
306  for(size_t i=0; i<nFlags; ++i) {
307  (*theFlag)[i] = (table) ? table->GetFlag(i) : true;
308  }
309  for(size_t i=nFlags; i<nCouples; ++i) {
310  G4bool yes = (table) ? table->GetFlag(i) : true;
311  theDensityFactor->push_back(1.0);
312  theDensityIdx->push_back(i);
313  theFlag->push_back(yes);
314  }
315  // use base materials
316  for(size_t i=0; i<nCouples; ++i) {
317  // base material is needed only for a couple which is not
318  // initialised and for which tables will be computed
319  auto couple = theCoupleTable->GetMaterialCutsCouple(i);
320  auto pcuts = couple->GetProductionCuts();
321  auto mat = couple->GetMaterial();
322  auto bmat = mat->GetBaseMaterial();
323 
324  // base material exists - find it and check if it can be reused
325  if(bmat) {
326  for(size_t j=0; j<nCouples; ++j) {
327  if(j == i) { continue; }
328  auto bcouple = theCoupleTable->GetMaterialCutsCouple(j);
329 
330  if(bcouple->GetMaterial() == bmat &&
331  bcouple->GetProductionCuts() == pcuts) {
332 
333  // based couple exist in the same region
334  (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity();
335  (*theDensityIdx)[i] = j;
336  (*theFlag)[i] = false;
337 
338  // ensure that there will no double initialisation
339  (*theDensityFactor)[j] = 1.0;
340  (*theDensityIdx)[j] = j;
341  (*theFlag)[j] = true;
342  break;
343  }
344  }
345  }
346  }
347  /*
348  for(size_t i=0; i<nCouples; ++i) {
349  G4cout << "CoupleIdx= " << i << " Flag= " << theFlag[i]
350  << " TableFlag= " << table->GetFlag(i) << " "
351  << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName()
352  << G4endl;
353  }
354  */
355 }
356 
357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
358 
361  G4VEmModel* model,
362  const G4ParticleDefinition* part,
364  G4bool spline)
365 {
366  // check input
368  if(!table) { return table; }
369  if(emin >= emax) {
370  table->clearAndDestroy();
371  delete table;
372  table = nullptr;
373  return table;
374  }
377  G4bool useMB = model->UseBaseMaterials();
378 
379  // Access to materials
380  const G4ProductionCutsTable* theCoupleTable=
382  size_t numOfCouples = theCoupleTable->GetTableSize();
383 
384  G4PhysicsLogVector* aVector = nullptr;
385 
386  for(size_t i=0; i<numOfCouples; ++i) {
387  if ((useMB && GetFlag(i)) || (!useMB && table->GetFlag(i))) {
388 
389  // create physics vector and fill it
390  auto couple = theCoupleTable->GetMaterialCutsCouple(i);
391  delete (*table)[i];
392 
393  // if start from zero then change the scale
394 
395  const G4Material* mat = couple->GetMaterial();
396 
397  G4double tmin = std::max(emin,model->MinPrimaryEnergy(mat,part));
398  if(0.0 >= tmin) { tmin = CLHEP::eV; }
399  G4int n = nbins;
400 
401  if(tmin >= emax) {
402  aVector = nullptr;
403  } else {
404  n *= (G4int)(std::log10(emax/tmin) + 0.5);
405  n = std::max(n, 3);
406  aVector = new G4PhysicsLogVector(tmin, emax, n);
407  }
408 
409  if(aVector) {
410  aVector->SetSpline(spline);
411  //G4cout << part->GetParticleName() << " in " << mat->GetName()
412  // << " tmin= " << tmin << G4endl;
413  for(G4int j=0; j<=n; ++j) {
414  aVector->PutValue(j, model->Value(couple, part,
415  aVector->Energy(j)));
416  }
417  if(spline) { aVector->FillSecondDerivatives(); }
418  }
419  G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
420  }
421  }
422  /*
423  G4cout << "G4LossTableBuilder::BuildTableForModel done for "
424  << part->GetParticleName() << " and "<< model->GetName()
425  << " " << table << G4endl;
426  */
427  return table;
428 }
429 
430 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
431