Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/utils/src/G4LossTableBuilder.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/utils/src/G4LossTableBuilder.cc (Version 11.3.0) and /processes/electromagnetic/utils/src/G4LossTableBuilder.cc (Version 7.1.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
                                                   >>  23 // $Id: G4LossTableBuilder.cc,v 1.17 2005/04/12 18:13:04 vnivanch Exp $
                                                   >>  24 // GEANT4 tag $Name: geant4-07-01-patch-01 $
                                                   >>  25 //
 26 // -------------------------------------------     26 // -------------------------------------------------------------------
 27 //                                                 27 //
 28 // GEANT4 Class file                               28 // GEANT4 Class file
 29 //                                                 29 //
 30 //                                                 30 //
 31 // File name:     G4LossTableBuilder               31 // File name:     G4LossTableBuilder
 32 //                                                 32 //
 33 // Author:        Vladimir Ivanchenko              33 // Author:        Vladimir Ivanchenko
 34 //                                                 34 //
 35 // Creation date: 03.01.2002                       35 // Creation date: 03.01.2002
 36 //                                                 36 //
 37 // Modifications:                                  37 // Modifications:
 38 //                                                 38 //
 39 // 23-01-03 V.Ivanchenko Cut per region            39 // 23-01-03 V.Ivanchenko Cut per region
 40 // 21-07-04 V.Ivanchenko Fix problem of range      40 // 21-07-04 V.Ivanchenko Fix problem of range for dedx=0
 41 // 08-11-04 Migration to new interface of Stor <<  41 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivantchenko)
 42 // 07-12-04 Fix of BuildDEDX table (V.Ivanchen <<  42 // 07-12-04 Fix of BuildDEDX table (V.Ivantchenko)
 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 //                                                 43 //
 48 // Class Description:                              44 // Class Description:
 49 //                                                 45 //
 50 // -------------------------------------------     46 // -------------------------------------------------------------------
 51 //                                                 47 //
 52 //....oooOO0OOooo........oooOO0OOooo........oo     48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 53 //....oooOO0OOooo........oooOO0OOooo........oo     49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 54                                                    50 
 55 #include "G4LossTableBuilder.hh"                   51 #include "G4LossTableBuilder.hh"
 56 #include "G4SystemOfUnits.hh"                  << 
 57 #include "G4PhysicsTable.hh"                       52 #include "G4PhysicsTable.hh"
 58 #include "G4PhysicsLogVector.hh"                   53 #include "G4PhysicsLogVector.hh"
 59 #include "G4PhysicsTableHelper.hh"                 54 #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                                                    55 
 74 //....oooOO0OOooo........oooOO0OOooo........oo     56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 75                                                    57 
 76 G4LossTableBuilder::G4LossTableBuilder(G4bool  <<  58 void G4LossTableBuilder::BuildDEDXTable(G4PhysicsTable* dedxTable,
 77   : isInitializer(master)                      <<  59                                   const std::vector<G4PhysicsTable*>& list)
 78 {                                                  60 {
 79   theParameters = G4EmParameters::Instance();  <<  61   size_t n_processes = list.size();
 80   if (nullptr == theFlag) {                    <<  62   if(1 >= n_processes) return;
 81     theDensityFactor = new std::vector<G4doubl << 
 82     theDensityIdx = new std::vector<G4int>;    << 
 83     theFlag = new std::vector<G4bool>;         << 
 84   }                                            << 
 85 }                                              << 
 86                                                    63 
 87 //....oooOO0OOooo........oooOO0OOooo........oo <<  64   size_t n_vectors = dedxTable->length();
                                                   >>  65   if(0 >= n_vectors) return;
 88                                                    66 
 89 G4LossTableBuilder::~G4LossTableBuilder()      <<  67   G4bool b;
 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                                                    68 
101 //....oooOO0OOooo........oooOO0OOooo........oo <<  69   for (size_t i=0; i<n_vectors; i++) {
102                                                    70 
103 const std::vector<G4int>* G4LossTableBuilder:: <<  71     G4PhysicsVector* pv = (*dedxTable)[i];
104 {                                              <<  72     size_t nbins = pv->GetVectorLength();
105   return theDensityIdx;                        << 
106 }                                              << 
107                                                << 
108 //....oooOO0OOooo........oooOO0OOooo........oo << 
109                                                << 
110 const std::vector<G4double>* G4LossTableBuilde << 
111 {                                              << 
112   return theDensityFactor;                     << 
113 }                                              << 
114                                                    73 
115 //....oooOO0OOooo........oooOO0OOooo........oo <<  74     for (size_t j=0; j<nbins; j++) {
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;                         75       G4double dedx = 0.0;
152       for (std::size_t k=0; k<n_processes; ++k <<  76       G4double energy = pv->GetLowEdgeEnergy(j);
153   const G4PhysicsVector* pv1 = (*(list[k]))[i] <<  77 
154   //if (0 == i) G4cout << "     " << k << ". p <<  78       for (size_t k=0; k<n_processes; k++) {
155   dedx += (*pv1)[j];                           <<  79   dedx += ((*(list[k]))[i])->GetValue(energy, b);
156       }                                            80       }
157       pv->PutValue(j, dedx);                       81       pv->PutValue(j, dedx);
158     }                                              82     }
159     if(splineFlag) { pv->FillSecondDerivatives << 
160     G4PhysicsTableHelper::SetPhysicsVector(ded << 
161   }                                                83   }
162   //G4cout << "### G4LossTableBuilder::BuildDE << 
163   //G4cout << *dedxTable << G4endl;            << 
164 }                                                  84 }
165                                                    85 
166 //....oooOO0OOooo........oooOO0OOooo........oo     86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
167                                                    87 
168 void G4LossTableBuilder::BuildRangeTable(const     88 void G4LossTableBuilder::BuildRangeTable(const G4PhysicsTable* dedxTable,
169                                          G4Phy <<  89                                                G4PhysicsTable* rangeTable)
170 // Build range table from the energy loss tabl     90 // Build range table from the energy loss table
171 {                                                  91 {
172   //G4cout << "### G4LossTableBuilder::BuildRa <<  92   size_t n_vectors = dedxTable->length();
173   //G4cout << *const_cast<G4PhysicsTable*>(ded <<  93   if(!n_vectors) return;
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                                                    94 
200     // initialisation of a new vector          <<  95   G4bool b;
201     if(npoints < 3) { npoints = 3; }           <<  96   size_t n = 100;
202                                                <<  97   G4double del = 1.0/(G4double)n;
203     delete (*rangeTable)[i];                   <<  98 
204     G4PhysicsLogVector* v;                     <<  99   for (size_t i=0; i<n_vectors; i++) {
205     if(0 == bin0) { v = new G4PhysicsLogVector << 100 
206     else { v = new G4PhysicsLogVector(elow, eh << 101     if (rangeTable->GetFlag(i)) {
207                                                << 102       G4PhysicsVector* pv = (*dedxTable)[i];
208     // assumed dedx proportional to beta       << 103       size_t nbins = pv->GetVectorLength();
209     G4double energy1 = v->Energy(0);           << 104       size_t bin0  = 0;
210     G4double range   = 2.*energy1/dedx1;       << 105       G4double elow = pv->GetLowEdgeEnergy(0);
211     /*                                         << 106       G4double ehigh = pv->GetLowEdgeEnergy(nbins);
212     G4cout << "New Range vector Npoints=" << v << 107       G4double dedx1  = pv->GetValue(elow, b);
213      << " coupleIdx=" << i << " spline=" << v- << 108 
214      << " Elow=" << v->GetMinEnergy() <<" Ehig << 109       if(dedx1 == 0.0) {
215      << " DEDX(Elow)=" << dedx1 << " R(Elow)=" << 110         for (size_t k=1; k<nbins; k++) {
216     */                                         << 111           bin0++;
217     v->PutValue(0,range);                      << 112           elow  = pv->GetLowEdgeEnergy(k);
218                                                << 113           dedx1 = pv->GetValue(elow, b);
219     for (std::size_t j=1; j<npoints; ++j) {    << 114           if(dedx1 > 0.0) break;
220                                                << 115         }
221       G4double energy2 = v->Energy(j);         << 116         nbins -= bin0;
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       }                                           117       }
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                                                   118 
247 //....oooOO0OOooo........oooOO0OOooo........oo << 119       G4PhysicsLogVector* v = new G4PhysicsLogVector(elow, ehigh, nbins);
248                                                   120 
249 void                                           << 121       G4double range  = 2.*elow/dedx1;
250 G4LossTableBuilder::BuildInverseRangeTable(con << 122       //G4double range  = elow/dedx1;
251                                            G4P << 123       v->PutValue(0,range);
252 // Build inverse range table from the energy l << 124       G4double energy1 = elow;
253 {                                              << 
254   std::size_t nCouples = rangeTable->size();   << 
255   if(0 >= nCouples) { return; }                << 
256                                                   125 
257   for (std::size_t i=0; i<nCouples; ++i) {     << 126       for (size_t j=1; j<nbins; j++) {
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                                                   127 
279 //....oooOO0OOooo........oooOO0OOooo........oo << 128         G4double energy2 = pv->GetLowEdgeEnergy(j+bin0);
                                                   >> 129         G4double dedx2   = pv->GetValue(energy2, b);
                                                   >> 130         G4double de      = (energy2 - energy1) * del;
                                                   >> 131         G4double energy  = energy1 - de*0.5;
280                                                   132 
281 void G4LossTableBuilder::InitialiseBaseMateria << 133         G4bool   yes     = true;
282 {                                              << 134         if(dedx1 < DBL_MIN || dedx2 < DBL_MIN) yes = false;   
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                                                   135 
304   if(nFlags != nCouples) { isInitialized = fal << 136         G4double fac, f;
305   if(isInitialized) { return; }                << 
306                                                   137 
307   // reserve memory                            << 138         if(yes) fac = std::log(dedx2/dedx1)/std::log(energy2/energy1);
308   theFlag->resize(nCouples, true);             << 139         else    fac = (dedx2 - dedx1)/(energy2 - energy1);
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                                                   140 
320   // use base materials                        << 141         for (size_t k=0; k<n; k++) {
321   for (G4int i=0; i<(G4int)nCouples; ++i) {    << 142           energy += de;
322     // base material is needed only for a coup << 143           if(yes) f = dedx1*std::exp(fac*std::log(energy/energy1));
323     // initialised and for which tables will b << 144           else    f = dedx1 + fac*(energy - energy1);
324     auto couple = theCoupleTable->GetMaterialC << 145           if(f > DBL_MIN) range  += de/f;
325     auto pcuts = couple->GetProductionCuts();  << 146         }
326     auto mat  = couple->GetMaterial();         << 147         v->PutValue(j,range);
327     auto bmat = mat->GetBaseMaterial();        << 148         energy1 = energy2;
328                                                << 149         dedx1   = dedx2;
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       }                                           150       }
                                                   >> 151       G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
350     }                                             152     }
351   }                                               153   }
352 }                                                 154 }
353                                                   155 
354 //....oooOO0OOooo........oooOO0OOooo........oo    156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
355                                                   157 
356 G4PhysicsTable*                                << 158 void G4LossTableBuilder::BuildInverseRangeTable(const G4PhysicsTable* rangeTable,
357 G4LossTableBuilder::BuildTableForModel(G4Physi << 159                                                       G4PhysicsTable* invRangeTable)
358                                        G4VEmMo << 160 // Build inverse range table from the energy loss table
359                                        const G << 
360                                        G4doubl << 
361                                        G4bool  << 
362 {                                                 161 {
363   // check input                               << 162   size_t n_vectors = rangeTable->length();
364   G4PhysicsTable* table = G4PhysicsTableHelper << 163   if(!n_vectors) return;
365   if (nullptr == table) { return table; }      << 164   G4bool b;
366   if (aTable != nullptr && aTable != table) {  << 165 
367     aTable->clearAndDestroy();                 << 166   for (size_t i=0; i<n_vectors; i++) {
368     delete aTable;                             << 167 
369   }                                            << 168     if (invRangeTable->GetFlag(i)) {
                                                   >> 169       G4PhysicsVector* pv = (*rangeTable)[i];
                                                   >> 170       size_t nbins   = pv->GetVectorLength();
                                                   >> 171       G4double elow  = pv->GetLowEdgeEnergy(0);
                                                   >> 172       G4double ehigh = pv->GetLowEdgeEnergy(nbins-1);
                                                   >> 173       G4double rlow  = pv->GetValue(elow, b);
                                                   >> 174       G4double rhigh = pv->GetValue(ehigh, b);
                                                   >> 175 
                                                   >> 176       rhigh *= std::exp(std::log(rhigh/rlow)/((G4double)(nbins-1)));
                                                   >> 177 
                                                   >> 178       G4PhysicsLogVector* v = new G4PhysicsLogVector(rlow, rhigh, nbins);
                                                   >> 179 
                                                   >> 180       v->PutValue(0,elow);
                                                   >> 181       G4double energy1 = elow;
                                                   >> 182       G4double range1  = rlow;
                                                   >> 183       G4double energy2 = elow;
                                                   >> 184       G4double range2  = rlow;
                                                   >> 185       size_t ilow      = 0;
                                                   >> 186       size_t ihigh;
                                                   >> 187 
                                                   >> 188       for (size_t j=1; j<nbins; j++) {
                                                   >> 189 
                                                   >> 190         G4double range = v->GetLowEdgeEnergy(j);
                                                   >> 191 
                                                   >> 192         for (ihigh=ilow+1; ihigh<nbins; ihigh++) {
                                                   >> 193           energy2 = pv->GetLowEdgeEnergy(ihigh);
                                                   >> 194           range2  = pv->GetValue(energy2, b);
                                                   >> 195           if(range2 >= range || ihigh == nbins-1) {
                                                   >> 196             ilow = ihigh - 1;
                                                   >> 197             energy1 = pv->GetLowEdgeEnergy(ilow);
                                                   >> 198             range1  = pv->GetValue(energy1, b);
                                                   >> 199             break;
                                                   >> 200     }
                                                   >> 201         }
370                                                   202 
371   InitialiseBaseMaterials(table);              << 203         G4double e = std::log(energy1) + 
372   G4int nbins = theParameters->NumberOfBinsPer << 204                      std::log(energy2/energy1)*std::log(range/range1)/std::log(range2/range1);
373                                                   205 
374   // Access to materials                       << 206         v->PutValue(j,std::exp(e));
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       }                                           207       }
407                                                << 208       G4PhysicsTableHelper::SetPhysicsVector(invRangeTable, i, v);
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     }                                             209     }
421   }                                               210   }
422   /*                                           << 
423   G4cout << "G4LossTableBuilder::BuildTableFor << 
424          << part->GetParticleName() << " and " << 
425          << "  " << table << G4endl;           << 
426   */                                           << 
427   //G4cout << *table << G4endl;                << 
428   return table;                                << 
429 }                                                 211 }
430                                                   212 
431 //....oooOO0OOooo........oooOO0OOooo........oo    213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 214 
432                                                   215