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 9.1.p3)


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