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.4.p4)


  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.32 2009-08-11 17:24:53 vnivanch Exp $
                                                   >>  27 // GEANT4 tag $Name: not supported by cvs2svn $
                                                   >>  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     49 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
 47 //                                                 50 //
 48 // Class Description:                              51 // Class Description:
 49 //                                                 52 //
 50 // -------------------------------------------     53 // -------------------------------------------------------------------
 51 //                                                 54 //
 52 //....oooOO0OOooo........oooOO0OOooo........oo     55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 53 //....oooOO0OOooo........oooOO0OOooo........oo     56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 54                                                    57 
 55 #include "G4LossTableBuilder.hh"                   58 #include "G4LossTableBuilder.hh"
 56 #include "G4SystemOfUnits.hh"                  << 
 57 #include "G4PhysicsTable.hh"                       59 #include "G4PhysicsTable.hh"
 58 #include "G4PhysicsLogVector.hh"                   60 #include "G4PhysicsLogVector.hh"
 59 #include "G4PhysicsTableHelper.hh"                 61 #include "G4PhysicsTableHelper.hh"
 60 #include "G4PhysicsFreeVector.hh"              <<  62 #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                                                    63 
 74 //....oooOO0OOooo........oooOO0OOooo........oo     64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 75                                                    65 
 76 G4LossTableBuilder::G4LossTableBuilder(G4bool  <<  66 G4LossTableBuilder::G4LossTableBuilder() 
 77   : isInitializer(master)                      << 
 78 {                                                  67 {
 79   theParameters = G4EmParameters::Instance();  <<  68   splineFlag = true;
 80   if (nullptr == theFlag) {                    << 
 81     theDensityFactor = new std::vector<G4doubl << 
 82     theDensityIdx = new std::vector<G4int>;    << 
 83     theFlag = new std::vector<G4bool>;         << 
 84   }                                            << 
 85 }                                                  69 }
 86                                                    70 
 87 //....oooOO0OOooo........oooOO0OOooo........oo     71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 88                                                    72 
 89 G4LossTableBuilder::~G4LossTableBuilder()          73 G4LossTableBuilder::~G4LossTableBuilder() 
 90 {                                              <<  74 {}
 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                                                    75 
129 //....oooOO0OOooo........oooOO0OOooo........oo     76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
130                                                    77 
131 void                                               78 void 
132 G4LossTableBuilder::BuildDEDXTable(G4PhysicsTa     79 G4LossTableBuilder::BuildDEDXTable(G4PhysicsTable* dedxTable,
133                                    const std:: <<  80            const std::vector<G4PhysicsTable*>& list)
134 {                                                  81 {
135   InitialiseBaseMaterials(dedxTable);          <<  82   size_t n_processes = list.size();
136   std::size_t n_processes = list.size();       <<  83   if(1 >= n_processes) return;
137   if(1 >= n_processes) { return; }             <<  84 
138                                                <<  85   size_t n_vectors = (list[0])->length();
139   std::size_t nCouples = dedxTable->size();    <<  86   if(0 >= n_vectors) return;
140   //G4cout << "Nproc= " << n_processes << " nC <<  87 
141   //   << dedxTable->size() << G4endl;         <<  88   G4PhysicsLogVector* pv0 = static_cast<G4PhysicsLogVector*>((*(list[0]))[0]);
142   if(0 >= nCouples) { return; }                <<  89   size_t npoints = pv0->GetVectorLength();
143                                                <<  90   for (size_t i=0; i<n_vectors; i++) {
144   for (std::size_t i=0; i<nCouples; ++i) {     <<  91 
145     auto pv0 = static_cast<G4PhysicsLogVector* <<  92     G4PhysicsLogVector* pv = new G4PhysicsLogVector(*pv0);
146     //if (0 == i) G4cout << i << ". pv0=" << p <<  93     //    pv = new G4PhysicsLogVector(elow, ehigh, npoints-1);
147     if(pv0 == nullptr) { continue; }           <<  94     pv->SetSpline(splineFlag);
148     std::size_t npoints = pv0->GetVectorLength <<  95     for (size_t j=0; j<npoints; j++) {
149     auto pv = new G4PhysicsLogVector(*pv0);    << 
150     for (std::size_t j=0; j<npoints; ++j) {    << 
151       G4double dedx = 0.0;                         96       G4double dedx = 0.0;
152       for (std::size_t k=0; k<n_processes; ++k <<  97       for (size_t k=0; k<n_processes; k++) {
153   const G4PhysicsVector* pv1 = (*(list[k]))[i] <<  98         G4PhysicsVector* pv1   = (*(list[k]))[i];
154   //if (0 == i) G4cout << "     " << k << ". p << 
155   dedx += (*pv1)[j];                               99   dedx += (*pv1)[j];
156       }                                           100       }
157       pv->PutValue(j, dedx);                      101       pv->PutValue(j, dedx);
158     }                                             102     }
159     if(splineFlag) { pv->FillSecondDerivatives << 103     if(splineFlag) pv->FillSecondDerivatives();
160     G4PhysicsTableHelper::SetPhysicsVector(ded    104     G4PhysicsTableHelper::SetPhysicsVector(dedxTable, i, pv);
161   }                                               105   }
162   //G4cout << "### G4LossTableBuilder::BuildDE << 
163   //G4cout << *dedxTable << G4endl;            << 
164 }                                                 106 }
165                                                   107 
166 //....oooOO0OOooo........oooOO0OOooo........oo    108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
167                                                   109 
168 void G4LossTableBuilder::BuildRangeTable(const    110 void G4LossTableBuilder::BuildRangeTable(const G4PhysicsTable* dedxTable,
169                                          G4Phy << 111            G4PhysicsTable* rangeTable,
                                                   >> 112            G4bool isIonisation)
170 // Build range table from the energy loss tabl    113 // Build range table from the energy loss table
171 {                                                 114 {
172   //G4cout << "### G4LossTableBuilder::BuildRa << 115   size_t n_vectors = dedxTable->length();
173   //G4cout << *const_cast<G4PhysicsTable*>(ded << 116   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                                                   117 
200     // initialisation of a new vector          << 118   size_t n = 100;
201     if(npoints < 3) { npoints = 3; }           << 119   G4double del = 1.0/(G4double)n;
202                                                   120 
203     delete (*rangeTable)[i];                   << 121   for (size_t i=0; i<n_vectors; 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                                                   122 
281 void G4LossTableBuilder::InitialiseBaseMateria << 123     if (rangeTable->GetFlag(i) || !isIonisation) {
282 {                                              << 124       G4PhysicsLogVector* pv = 
283   if(!isInitializer) { return; }               << 125   static_cast<G4PhysicsLogVector*>((*dedxTable)[i]);
284   const G4ProductionCutsTable* theCoupleTable= << 126       size_t npoints = pv->GetVectorLength();
285     G4ProductionCutsTable::GetProductionCutsTa << 127       size_t bin0    = 0;
286   std::size_t nCouples = theCoupleTable->GetTa << 128       G4double elow  = pv->Energy(0);
287   std::size_t nFlags = theFlag->size();        << 129       G4double ehigh = pv->Energy(npoints-1);
288   /*                                           << 130       G4double dedx1 = pv->Value(elow);
289   G4cout << "### InitialiseBaseMaterials: nCou << 131 
290    << " nFlags=" << nFlags << " isInit:" << is << 132       //G4cout << "nbins= " << nbins << " dedx1= " << dedx1 << G4endl;
291    << " baseMat:" << baseMatFlag << G4endl;    << 133 
292   */                                           << 134       // protection for specific cases dedx=0
293   // define base material flag                 << 135       if(dedx1 == 0.0) {
294   if(isBaseMatActive && !baseMatFlag) {        << 136         for (size_t k=1; k<npoints; k++) {
295     for(G4int i=0; i<(G4int)nCouples; ++i) {   << 137           bin0++;
296       if(nullptr != theCoupleTable->GetMateria << 138           elow  = pv->Energy(k);
297   baseMatFlag = true;                          << 139           dedx1 = (*pv)[k];
298         isInitialized = false;                 << 140           if(dedx1 > 0.0) break;
299   break;                                       << 141         }
                                                   >> 142         npoints -= bin0;
300       }                                           143       }
301     }                                          << 144       // G4cout<<"New Range vector" << G4endl;
302   }                                            << 145       // G4cout<<"nbins= "<<npoints-1<<" elow= "<<elow<<" ehigh= "<<ehigh<<G4endl;
303                                                << 146       // initialisation of a new vector
304   if(nFlags != nCouples) { isInitialized = fal << 147       if(npoints < 2) npoints = 2;
305   if(isInitialized) { return; }                << 148       G4PhysicsLogVector* v;
306                                                << 149       if(0 == bin0) { v = new G4PhysicsLogVector(*pv); }
307   // reserve memory                            << 150       else          { v = new G4PhysicsLogVector(elow, ehigh, npoints-1); }
308   theFlag->resize(nCouples, true);             << 151       // dedx is exect zero
309   theDensityFactor->resize(nCouples,1.0);      << 152       if(2 == npoints) {
310   theDensityIdx->resize(nCouples, 0);          << 153   v->PutValue(0,1000.);
311                                                << 154   v->PutValue(1,2000.);
312   // define default flag and index of used mat << 155   G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
313   for (G4int i=0; i<(G4int)nCouples; ++i) {    << 156   return;
314     (*theFlag)[i] = (nullptr == table) ? true  << 157       }
315     (*theDensityIdx)[i] = i;                   << 158       v->SetSpline(splineFlag);
316   }                                            << 159 
317   isInitialized = true;                        << 160       // assumed dedx proportional to beta
318   if (!baseMatFlag) { return; }                << 161       G4double range  = 2.*elow/dedx1;
319                                                << 162       v->PutValue(0,range);
320   // use base materials                        << 163       G4double energy1 = elow;
321   for (G4int i=0; i<(G4int)nCouples; ++i) {    << 164 
322     // base material is needed only for a coup << 165       for (size_t j=1; j<npoints; j++) {
323     // initialised and for which tables will b << 166 
324     auto couple = theCoupleTable->GetMaterialC << 167         G4double energy2 = pv->Energy(j+bin0);
325     auto pcuts = couple->GetProductionCuts();  << 168         G4double de      = (energy2 - energy1) * del;
326     auto mat  = couple->GetMaterial();         << 169         G4double energy  = energy2 + de*0.5;
327     auto bmat = mat->GetBaseMaterial();        << 170         G4double sum = 0.0;
328                                                << 171         for (size_t k=0; k<n; k++) {
329     // base material exists - find it and chec << 172           energy -= de;
330     if(nullptr != bmat) {                      << 173           dedx1 = pv->Value(energy);
331       for(G4int j=0; j<(G4int)nCouples; ++j) { << 174           if(dedx1 > 0.0) sum += de/dedx1;
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   }                                               175   }
                                                   >> 176         range += sum;
                                                   >> 177   //  G4cout << "Range i= " <<i << " j= " << j << G4endl;
                                                   >> 178         v->PutValue(j,range);
                                                   >> 179         energy1 = energy2;
349       }                                           180       }
                                                   >> 181       if(splineFlag) v->FillSecondDerivatives();
                                                   >> 182       G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
350     }                                             183     }
351   }                                               184   }
352 }                                                 185 }
353                                                   186 
354 //....oooOO0OOooo........oooOO0OOooo........oo    187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
355                                                   188 
356 G4PhysicsTable*                                << 189 void G4LossTableBuilder::BuildInverseRangeTable(const G4PhysicsTable* rangeTable,
357 G4LossTableBuilder::BuildTableForModel(G4Physi << 190             G4PhysicsTable* invRangeTable,
358                                        G4VEmMo << 191             G4bool isIonisation)
359                                        const G << 192 // Build inverse range table from the energy loss table
360                                        G4doubl << 
361                                        G4bool  << 
362 {                                                 193 {
363   // check input                               << 194   size_t n_vectors = rangeTable->length();
364   G4PhysicsTable* table = G4PhysicsTableHelper << 195   if(!n_vectors) return;
365   if (nullptr == table) { return table; }      << 
366   if (aTable != nullptr && aTable != table) {  << 
367     aTable->clearAndDestroy();                 << 
368     delete aTable;                             << 
369   }                                            << 
370                                                   196 
371   InitialiseBaseMaterials(table);              << 197   for (size_t i=0; i<n_vectors; i++) {
372   G4int nbins = theParameters->NumberOfBinsPer << 
373                                                   198 
374   // Access to materials                       << 199     if (invRangeTable->GetFlag(i) || !isIonisation) {
375   const G4ProductionCutsTable* theCoupleTable= << 200       G4PhysicsVector* pv = (*rangeTable)[i];
376         G4ProductionCutsTable::GetProductionCu << 201       size_t npoints = pv->GetVectorLength();
377   std::size_t numOfCouples = theCoupleTable->G << 202       G4double rlow  = (*pv)[0];
378   /*                                           << 203       G4double rhigh = (*pv)[npoints-1];
379   G4cout << "   G4LossTableBuilder::BuildTable << 204       
380    << " isMaster=" << isInitializer << " model << 205       G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(npoints,rlow,rhigh);
381    << "  " << part->GetParticleName() << G4end << 206       v->SetSpline(splineFlag);
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                                                   207 
408       if(nullptr != aVector) {                 << 208       for (size_t j=0; j<npoints; j++) {
409         //G4cout << part->GetParticleName() << << 209   G4double e  = pv->Energy(j);
410   //     << " emin= " << tmin << " emax=" << e << 210   G4double r  = (*pv)[j];
411         for(G4int j=0; j<=n; ++j) {            << 211         v->PutValues(j,r,e);
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       }                                           212       }
419       G4PhysicsTableHelper::SetPhysicsVector(t << 213       if(splineFlag) v->FillSecondDerivatives();
                                                   >> 214 
                                                   >> 215       G4PhysicsTableHelper::SetPhysicsVector(invRangeTable, i, v);
420     }                                             216     }
421   }                                               217   }
422   /*                                           << 
423   G4cout << "G4LossTableBuilder::BuildTableFor << 
424          << part->GetParticleName() << " and " << 
425          << "  " << table << G4endl;           << 
426   */                                           << 
427   //G4cout << *table << G4endl;                << 
428   return table;                                << 
429 }                                                 218 }
430                                                   219 
431 //....oooOO0OOooo........oooOO0OOooo........oo    220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 221 
432                                                   222