Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4SBBremTable.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/standard/src/G4SBBremTable.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4SBBremTable.cc (Version 11.0.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 //                                                 26 //
 27 // -------------------------------------------     27 // -------------------------------------------------------------------
 28 //                                                 28 //
 29 // GEANT4 Class file                               29 // GEANT4 Class file
 30 //                                                 30 //
 31 //                                                 31 //
 32 // File name:     G4SBBremTable                    32 // File name:     G4SBBremTable
 33 //                                                 33 //
 34 // Author:        Mihaly Novak                     34 // Author:        Mihaly Novak
 35 //                                                 35 //
 36 // Creation date: 15.07.2018                       36 // Creation date: 15.07.2018
 37 //                                                 37 //
 38 // Modifications:                                  38 // Modifications:
 39 //                                                 39 //
 40 // -------------------------------------------     40 // -------------------------------------------------------------------
 41 //                                                 41 //
 42 #include "G4SBBremTable.hh"                        42 #include "G4SBBremTable.hh"
 43                                                    43 
 44 #include "G4SystemOfUnits.hh"                      44 #include "G4SystemOfUnits.hh"
 45                                                    45 
 46 #include "G4Material.hh"                           46 #include "G4Material.hh"
 47 #include "G4ProductionCutsTable.hh"                47 #include "G4ProductionCutsTable.hh"
 48 #include "G4MaterialCutsCouple.hh"                 48 #include "G4MaterialCutsCouple.hh"
 49 #include "Randomize.hh"                            49 #include "Randomize.hh"
 50 #include "G4EmParameters.hh"                   << 
 51                                                    50 
 52 #include "G4String.hh"                             51 #include "G4String.hh"
 53                                                    52 
 54 #include "G4Log.hh"                                53 #include "G4Log.hh"
 55 #include "G4Exp.hh"                                54 #include "G4Exp.hh"
 56                                                    55 
 57 #include "zlib.h"                                  56 #include "zlib.h"
 58                                                    57 
 59 #include <iostream>                                58 #include <iostream>
 60 #include <fstream>                                 59 #include <fstream>
 61 #include <sstream>                                 60 #include <sstream>
 62 #include <algorithm>                               61 #include <algorithm>
 63                                                    62 
 64 G4SBBremTable::G4SBBremTable()                     63 G4SBBremTable::G4SBBremTable()
 65  : fMaxZet(-1), fNumElEnergy(-1), fNumKappa(-1     64  : fMaxZet(-1), fNumElEnergy(-1), fNumKappa(-1), fUsedLowEenergy(-1.),
 66    fUsedHighEenergy(-1.), fLogMinElEnergy(-1.)     65    fUsedHighEenergy(-1.), fLogMinElEnergy(-1.), fILDeltaElEnergy(-1.)
 67 {}                                                 66 {}
 68                                                    67 
 69 G4SBBremTable::~G4SBBremTable()                    68 G4SBBremTable::~G4SBBremTable()
 70 {                                                  69 {
 71   ClearSamplingTables();                           70   ClearSamplingTables();
 72 }                                                  71 }
 73                                                    72 
 74 void G4SBBremTable::Initialize(const G4double      73 void G4SBBremTable::Initialize(const G4double lowe, const G4double highe)
 75 {                                                  74 {
 76   fUsedLowEenergy  = lowe;                         75   fUsedLowEenergy  = lowe;
 77   fUsedHighEenergy = highe;                        76   fUsedHighEenergy = highe;
 78   BuildSamplingTables();                           77   BuildSamplingTables();
 79   InitSamplingTables();                            78   InitSamplingTables();
 80 //  Dump();                                        79 //  Dump();
 81 }                                                  80 }
 82                                                    81 
 83 // run-time method that samples energy transfe     82 // run-time method that samples energy transferred to the emitted gamma photon
 84 double G4SBBremTable::SampleEnergyTransfer(con     83 double G4SBBremTable::SampleEnergyTransfer(const G4double eekin,
 85                                            con     84                                            const G4double leekin,
 86                                            con     85                                            const G4double gcut,
 87                                            con     86                                            const G4double dielSupConst,
 88                                            con     87                                            const G4int    iZet,
 89                                            con     88                                            const G4int    matCutIndx,
 90                                            con     89                                            const G4bool   isElectron)
 91 {                                                  90 {
 92   static const G4double kAlpha2Pi = CLHEP::two     91   static const G4double kAlpha2Pi = CLHEP::twopi*CLHEP::fine_structure_const;
 93   const G4double zet = (G4double)iZet;             92   const G4double zet = (G4double)iZet;
 94   const G4int   izet = std::max(std::min(fMaxZ     93   const G4int   izet = std::max(std::min(fMaxZet, iZet),1);
 95   G4double eGamma    = 0.;                         94   G4double eGamma    = 0.;
 96   // this should be checked in the caller          95   // this should be checked in the caller
 97   // if (eekin<=gcut) return kappa;                96   // if (eekin<=gcut) return kappa;
 98   const G4double lElEnergy     = leekin;           97   const G4double lElEnergy     = leekin;
 99   const SamplingTablePerZ* stZ = fSBSamplingTa     98   const SamplingTablePerZ* stZ = fSBSamplingTables[izet];
100   // get the gamma cut of this Z that correspo     99   // get the gamma cut of this Z that corresponds to the current mat-cuts
101   const std::size_t gamCutIndx = stZ->fMatCutI << 100   const size_t gamCutIndx = stZ->fMatCutIndxToGamCutIndx[matCutIndx];
102   // gcut was not found: should never happen (    101   // gcut was not found: should never happen (only in verbose mode)
103   if (gamCutIndx >= stZ->fNumGammaCuts || stZ-    102   if (gamCutIndx >= stZ->fNumGammaCuts || stZ->fGammaECuts[gamCutIndx]!=gcut) {
104     G4String msg = " Gamma cut="+std::to_strin    103     G4String msg = " Gamma cut="+std::to_string(gcut) + " [MeV] was not found ";
105     msg += "in case of Z = " + std::to_string(    104     msg += "in case of Z = " + std::to_string(iZet) + ". ";
106     G4Exception("G4SBBremTable::SampleEnergyTr    105     G4Exception("G4SBBremTable::SampleEnergyTransfer()","em0X",FatalException,
107                 msg.c_str());                     106                 msg.c_str());
108   }                                               107   }
109   const G4double lGCut = stZ->fLogGammaECuts[g    108   const G4double lGCut = stZ->fLogGammaECuts[gamCutIndx];
110   // get the random engine                        109   // get the random engine
111   CLHEP::HepRandomEngine* rndmEngine = G4Rando    110   CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
112   // find lower e- energy bin                     111   // find lower e- energy bin
113   G4bool isCorner = false; // indicate that th    112   G4bool isCorner = false; // indicate that the lower edge e- energy < gam-gut
114   G4bool isSimply = false; // simply sampling:    113   G4bool isSimply = false; // simply sampling: isCorner+lower egde is selected
115   G4int elEnergyIndx   = stZ->fMaxElEnergyIndx    114   G4int elEnergyIndx   = stZ->fMaxElEnergyIndx;
116   // only if e- ekin is below the maximum valu    115   // only if e- ekin is below the maximum value(use table at maximum otherwise)
117   if (eekin < fElEnergyVect[elEnergyIndx]) {      116   if (eekin < fElEnergyVect[elEnergyIndx]) {
118     const G4double val = (lElEnergy-fLogMinElE    117     const G4double val = (lElEnergy-fLogMinElEnergy)*fILDeltaElEnergy;
119     elEnergyIndx       = (G4int)val;              118     elEnergyIndx       = (G4int)val;
120     G4double pIndxH    = val-elEnergyIndx;        119     G4double pIndxH    = val-elEnergyIndx;
121     // check if we are at limiting case: lower    120     // check if we are at limiting case: lower edge e- energy < gam-gut
122     if (fElEnergyVect[elEnergyIndx]<=gcut) {      121     if (fElEnergyVect[elEnergyIndx]<=gcut) {
123       // recompute the probability of taking t    122       // recompute the probability of taking the higher e- energy bin()
124       pIndxH   = (lElEnergy-lGCut)/(fLElEnergy    123       pIndxH   = (lElEnergy-lGCut)/(fLElEnergyVect[elEnergyIndx+1]-lGCut);
125       isCorner = true;                            124       isCorner = true;
126     }                                             125     }
127     //                                            126     //
128     if (rndmEngine->flat()<pIndxH) {              127     if (rndmEngine->flat()<pIndxH) {
129       ++elEnergyIndx;      // take the table a    128       ++elEnergyIndx;      // take the table at the higher e- energy bin
130     } else if (isCorner) { // take the table a    129     } else if (isCorner) { // take the table at the lower  e- energy bin
131       // special sampling need to be done if l    130       // special sampling need to be done if lower edge e- energy < gam-gut:
132       // actually, we "sample" from a table "b    131       // actually, we "sample" from a table "built" at the gamm-cut i.e. delta
133       isSimply = true;                            132       isSimply = true;
134     }                                             133     }
135   }                                               134   }
136   // should never happen under normal conditio    135   // should never happen under normal conditions but add protection
137   if (!stZ->fTablesPerEnergy[elEnergyIndx]) {     136   if (!stZ->fTablesPerEnergy[elEnergyIndx]) {
138     return 0.;                                    137     return 0.;
139   }                                               138   }
140   // Do the photon energy sampling:               139   // Do the photon energy sampling:
141   const STable *st =  stZ->fTablesPerEnergy[el    140   const STable *st =  stZ->fTablesPerEnergy[elEnergyIndx];
142   const std::vector<G4double>& cVect = st->fCu    141   const std::vector<G4double>& cVect = st->fCumCutValues;
143   const std::vector<STPoint>&  pVect = st->fST    142   const std::vector<STPoint>&  pVect = st->fSTable;
144   const G4double minVal = cVect[gamCutIndx];      143   const G4double minVal = cVect[gamCutIndx];
145                                                   144 
146   // should never happen under normal conditio    145   // should never happen under normal conditions but add protection
147   if (minVal >= 1.) {                             146   if (minVal >= 1.) {
148     return 0.;                                    147     return 0.;
149   }                                               148   }
150   // some transfomrmtion variables used in the    149   // some transfomrmtion variables used in the looop
151   const G4double lCurKappaC  = lGCut - leekin;    150   const G4double lCurKappaC  = lGCut - leekin;
152   const G4double lUsedKappaC = lGCut - fLElEne    151   const G4double lUsedKappaC = lGCut - fLElEnergyVect[elEnergyIndx];
153   // dielectric (always) and e+ correction sup    152   // dielectric (always) and e+ correction suppressions (if the primary is e+)
154   G4double suppression = 1.;                      153   G4double suppression = 1.;
155   G4double rndm[2];                               154   G4double rndm[2];
156   // rejection loop starts here                   155   // rejection loop starts here
157   do {                                            156   do {
158     rndmEngine->flatArray(2, rndm);               157     rndmEngine->flatArray(2, rndm);
159     G4double kappa = 1.0;                         158     G4double kappa = 1.0;
160     if (!isSimply) {                              159     if (!isSimply) {
161       const G4double cumRV = rndm[0]*(1.-minVa    160       const G4double cumRV = rndm[0]*(1.-minVal)+minVal;
162       // find lower index of the values in the    161       // find lower index of the values in the Cumulative Function: use linear
163       // instead of binary search because it's    162       // instead of binary search because it's faster in our case
164       const G4int cumLIndx = LinSearch(pVect,     163       const G4int cumLIndx = LinSearch(pVect, fNumKappa, cumRV)-1;
165 //      const G4int cumLIndx = std::lower_boun    164 //      const G4int cumLIndx = std::lower_bound( pVect.begin(), pVect.end(), cumRV,
166 //                                    [](const    165 //                                    [](const STPoint& p, const double& cumV) {
167 //                                    return p    166 //                                    return p.fCum<cumV; } ) - pVect.begin() -1;
168       const STPoint& stPL  = pVect[cumLIndx];     167       const STPoint& stPL  = pVect[cumLIndx];
169       const G4double pA    = stPL.fParA;          168       const G4double pA    = stPL.fParA;
170       const G4double pB    = stPL.fParB;          169       const G4double pB    = stPL.fParB;
171       const G4double cumL  = stPL.fCum;           170       const G4double cumL  = stPL.fCum;
172       const G4double cumH  = pVect[cumLIndx+1]    171       const G4double cumH  = pVect[cumLIndx+1].fCum;
173       const G4double lKL   = fLKappaVect[cumLI    172       const G4double lKL   = fLKappaVect[cumLIndx];
174       const G4double lKH   = fLKappaVect[cumLI    173       const G4double lKH   = fLKappaVect[cumLIndx+1];
175       const G4double dm1   = (cumRV-cumL)/(cum    174       const G4double dm1   = (cumRV-cumL)/(cumH-cumL);
176       const G4double dm2   = (1.+pA+pB)*dm1;      175       const G4double dm2   = (1.+pA+pB)*dm1;
177       const G4double dm3   = 1.+dm1*(pA+pB*dm1    176       const G4double dm3   = 1.+dm1*(pA+pB*dm1);
178       // kappa sampled at E_i e- energy           177       // kappa sampled at E_i e- energy
179       const G4double lKappa = lKL+dm2/dm3*(lKH    178       const G4double lKappa = lKL+dm2/dm3*(lKH-lKL);
180       // transform lKappa to [log(gcut/eekin),    179       // transform lKappa to [log(gcut/eekin),0] form [log(gcut/E_i),0]
181       kappa  = G4Exp(lKappa*lCurKappaC/lUsedKa    180       kappa  = G4Exp(lKappa*lCurKappaC/lUsedKappaC);
182     } else {                                      181     } else {
183 //      const G4double upLimit = std::min(1.*C    182 //      const G4double upLimit = std::min(1.*CLHEP::eV,eekin-gcut);
184 //      kappa = (gcut+rndm[0]*upLimit)/eekin;     183 //      kappa = (gcut+rndm[0]*upLimit)/eekin;
185       kappa = 1.-rndm[0]*(1.-gcut/eekin);         184       kappa = 1.-rndm[0]*(1.-gcut/eekin);
186     }                                             185     }
187     // compute the emitted photon energy: k       186     // compute the emitted photon energy: k
188     eGamma = kappa*eekin;                         187     eGamma = kappa*eekin;
189     const G4double invEGamma = 1./eGamma;         188     const G4double invEGamma = 1./eGamma;
190     // compute dielectric suppression: 1/(1+[g    189     // compute dielectric suppression: 1/(1+[gk_p/k]^2)
191     suppression = 1./(1.+dielSupConst*invEGamm    190     suppression = 1./(1.+dielSupConst*invEGamma*invEGamma);
192     // add positron correction if particle is     191     // add positron correction if particle is e+
193     if (!isElectron) {                            192     if (!isElectron) {
194       const G4double e1     = eekin - gcut;       193       const G4double e1     = eekin - gcut;
195       const G4double iBeta1 =  (e1 + CLHEP::el    194       const G4double iBeta1 =  (e1 + CLHEP::electron_mass_c2)
196                               / std::sqrt(e1*(    195                               / std::sqrt(e1*(e1 + 2.*CLHEP::electron_mass_c2));
197       const G4double e2     = eekin - eGamma;     196       const G4double e2     = eekin - eGamma;
198       const G4double iBeta2 =  (e2 + CLHEP::el    197       const G4double iBeta2 =  (e2 + CLHEP::electron_mass_c2)
199                               / std::sqrt(e2*(    198                               / std::sqrt(e2*(e2 + 2.*CLHEP::electron_mass_c2));
200       const G4double dum    = kAlpha2Pi*zet*(i    199       const G4double dum    = kAlpha2Pi*zet*(iBeta1 - iBeta2);
201       suppression = (dum > -12.) ? suppression    200       suppression = (dum > -12.) ? suppression*G4Exp(dum) : 0.;
202     }                                             201     }
203   } while (rndm[1] > suppression);                202   } while (rndm[1] > suppression);
204   // end of rejection loop                        203   // end of rejection loop
205   // return the sampled photon energy value k     204   // return the sampled photon energy value k
206   return eGamma;                                  205   return eGamma;
207 }                                                 206 }
208                                                   207 
209                                                   208 
210 void G4SBBremTable::BuildSamplingTables() {       209 void G4SBBremTable::BuildSamplingTables() {
211   // claer                                        210   // claer
212   ClearSamplingTables();                          211   ClearSamplingTables();
213   LoadSTGrid();                                   212   LoadSTGrid();
214   // First elements and gamma cuts data struct    213   // First elements and gamma cuts data structures need to be built:
215   // loop over all material-cuts and add gamma    214   // loop over all material-cuts and add gamma cut to the list of elements
216   const G4ProductionCutsTable                     215   const G4ProductionCutsTable
217   *thePCTable = G4ProductionCutsTable::GetProd    216   *thePCTable = G4ProductionCutsTable::GetProductionCutsTable();
218   // a temporary vector to store one element      217   // a temporary vector to store one element
219   std::vector<std::size_t> vtmp(1,0);          << 218   std::vector<size_t> vtmp(1,0);
220   std::size_t numMatCuts = thePCTable->GetTabl << 219   size_t numMatCuts = thePCTable->GetTableSize();
221   for (G4int imc=0; imc<(G4int)numMatCuts; ++i << 220   for (size_t imc=0; imc<numMatCuts; ++imc) {
222     const G4MaterialCutsCouple *matCut = thePC    221     const G4MaterialCutsCouple *matCut = thePCTable->GetMaterialCutsCouple(imc);
223     if (!matCut->IsUsed()) {                      222     if (!matCut->IsUsed()) {
224       continue;                                   223       continue;
225     }                                             224     }
226     const G4Material*           mat = matCut->    225     const G4Material*           mat = matCut->GetMaterial();
227     const G4ElementVector* elemVect = mat->Get    226     const G4ElementVector* elemVect = mat->GetElementVector();
228     const G4int              indxMC = matCut->    227     const G4int              indxMC = matCut->GetIndex();
229     const G4double gamCut = (*(thePCTable->Get    228     const G4double gamCut = (*(thePCTable->GetEnergyCutsVector(0)))[indxMC];
230     const std::size_t numElems = elemVect->siz << 229     const size_t numElems = elemVect->size();
231     for (std::size_t ielem=0; ielem<numElems;  << 230     for (size_t ielem=0; ielem<numElems; ++ielem) {
232       const G4Element *elem = (*elemVect)[iele    231       const G4Element *elem = (*elemVect)[ielem];
233       const G4int izet = std::max(std::min(fMa    232       const G4int izet = std::max(std::min(fMaxZet, elem->GetZasInt()),1);
234       if (!fSBSamplingTables[izet]) {             233       if (!fSBSamplingTables[izet]) {
235         // create data structure but do not lo    234         // create data structure but do not load sampling tables yet: will be
236         // loaded after we know what are the m    235         // loaded after we know what are the min/max e- energies where sampling
237         // will be needed during the simulatio    236         // will be needed during the simulation for this Z
238         // LoadSamplingTables(izet);              237         // LoadSamplingTables(izet);
239         fSBSamplingTables[izet] = new Sampling    238         fSBSamplingTables[izet] = new SamplingTablePerZ();
240       }                                           239       }
241       // add current gamma cut to the list of     240       // add current gamma cut to the list of this element data (only if this
242       // cut value is still not tehre)            241       // cut value is still not tehre)
243       const std::vector<double> &cVect = fSBSa    242       const std::vector<double> &cVect = fSBSamplingTables[izet]->fGammaECuts;
244       std::size_t indx = std::find(cVect.cbegi << 243       size_t indx = std::find(cVect.begin(), cVect.end(), gamCut)-cVect.begin();
245       if (indx==cVect.size()) {                   244       if (indx==cVect.size()) {
246         vtmp[0] = imc;                            245         vtmp[0] = imc;
247         fSBSamplingTables[izet]->fGamCutIndxTo    246         fSBSamplingTables[izet]->fGamCutIndxToMatCutIndx.push_back(vtmp);
248         fSBSamplingTables[izet]->fGammaECuts.p    247         fSBSamplingTables[izet]->fGammaECuts.push_back(gamCut);
249         fSBSamplingTables[izet]->fLogGammaECut    248         fSBSamplingTables[izet]->fLogGammaECuts.push_back(G4Log(gamCut));
250         ++fSBSamplingTables[izet]->fNumGammaCu    249         ++fSBSamplingTables[izet]->fNumGammaCuts;
251       } else {                                    250       } else {
252         fSBSamplingTables[izet]->fGamCutIndxTo    251         fSBSamplingTables[izet]->fGamCutIndxToMatCutIndx[indx].push_back(imc);
253       }                                           252       }
254     }                                             253     }
255   }                                               254   }
256 }                                                 255 }
257                                                   256 
258 void G4SBBremTable::InitSamplingTables() {        257 void G4SBBremTable::InitSamplingTables() {
259   const std::size_t numMatCuts = G4ProductionC << 258   const size_t numMatCuts = G4ProductionCutsTable::GetProductionCutsTable()
260                             ->GetTableSize();     259                             ->GetTableSize();
261   for (G4int iz=1; iz<fMaxZet+1; ++iz) {          260   for (G4int iz=1; iz<fMaxZet+1; ++iz) {
262     SamplingTablePerZ* stZ = fSBSamplingTables    261     SamplingTablePerZ* stZ = fSBSamplingTables[iz];
263     if (!stZ) continue;                           262     if (!stZ) continue;
264     // Load-in sampling table data:               263     // Load-in sampling table data:
265     LoadSamplingTables(iz);                       264     LoadSamplingTables(iz);
266     // init data                                  265     // init data
267     for (G4int iee=0; iee<fNumElEnergy; ++iee)    266     for (G4int iee=0; iee<fNumElEnergy; ++iee) {
268       if (!stZ->fTablesPerEnergy[iee])            267       if (!stZ->fTablesPerEnergy[iee])
269         continue;                                 268         continue;
270       const G4double elEnergy = fElEnergyVect[    269       const G4double elEnergy = fElEnergyVect[iee];
271       // 1 indicates that gamma production is     270       // 1 indicates that gamma production is not possible at this e- energy
272       stZ->fTablesPerEnergy[iee]->fCumCutValue    271       stZ->fTablesPerEnergy[iee]->fCumCutValues.resize(stZ->fNumGammaCuts,1.);
273       // sort gamma cuts and other members acc    272       // sort gamma cuts and other members accordingly
274       for (std::size_t i=0; i<stZ->fNumGammaCu << 273       for (size_t i=0; i<stZ->fNumGammaCuts-1; ++i) {
275         for (std::size_t j=i+1; j<stZ->fNumGam << 274         for (size_t j=i+1; j<stZ->fNumGammaCuts; ++j) {
276           if (stZ->fGammaECuts[j]<stZ->fGammaE    275           if (stZ->fGammaECuts[j]<stZ->fGammaECuts[i]) {
277             G4double dum0                   =     276             G4double dum0                   = stZ->fGammaECuts[i];
278             G4double dum1                   =     277             G4double dum1                   = stZ->fLogGammaECuts[i];
279             std::vector<std::size_t>   dumv =  << 278             std::vector<size_t>   dumv      = stZ->fGamCutIndxToMatCutIndx[i];
280             stZ->fGammaECuts[i]             =     279             stZ->fGammaECuts[i]             = stZ->fGammaECuts[j];
281             stZ->fLogGammaECuts[i]          =     280             stZ->fLogGammaECuts[i]          = stZ->fLogGammaECuts[j];
282             stZ->fGamCutIndxToMatCutIndx[i] =     281             stZ->fGamCutIndxToMatCutIndx[i] = stZ->fGamCutIndxToMatCutIndx[j];
283             stZ->fGammaECuts[j]             =     282             stZ->fGammaECuts[j]             = dum0;
284             stZ->fLogGammaECuts[j]          =     283             stZ->fLogGammaECuts[j]          = dum1;
285             stZ->fGamCutIndxToMatCutIndx[j] =  << 284             stZ->fGamCutIndxToMatCutIndx[j] = dumv;
286           }                                       285           }
287         }                                         286         }
288       }                                           287       }
289       // set couple indices to store the corre    288       // set couple indices to store the corresponding gamma cut index
290       stZ->fMatCutIndxToGamCutIndx.resize(numM    289       stZ->fMatCutIndxToGamCutIndx.resize(numMatCuts,-1);
291       for (std::size_t i=0; i<stZ->fGamCutIndx << 290       for (size_t i=0; i<stZ->fGamCutIndxToMatCutIndx.size(); ++i) {
292         for (std::size_t j=0; j<stZ->fGamCutIn << 291         for (size_t j=0; j<stZ->fGamCutIndxToMatCutIndx[i].size(); ++j) {
293           stZ->fMatCutIndxToGamCutIndx[stZ->fG    292           stZ->fMatCutIndxToGamCutIndx[stZ->fGamCutIndxToMatCutIndx[i][j]] = i;
294         }                                         293         }
295       }                                           294       }
296       // clear temporary vector                   295       // clear temporary vector
297       for (std::size_t i=0; i<stZ->fGamCutIndx << 296       for (size_t i=0; i<stZ->fGamCutIndxToMatCutIndx.size(); ++i) {
298         stZ->fGamCutIndxToMatCutIndx[i].clear(    297         stZ->fGamCutIndxToMatCutIndx[i].clear();
299       }                                           298       }
300       stZ->fGamCutIndxToMatCutIndx.clear();       299       stZ->fGamCutIndxToMatCutIndx.clear();
301       //  init for each gamma cut that are bel    300       //  init for each gamma cut that are below the e- energy
302       for (std::size_t ic=0; ic<stZ->fNumGamma << 301       for (size_t ic=0; ic<stZ->fNumGammaCuts; ++ic) {
303         const G4double gamCut = stZ->fGammaECu    302         const G4double gamCut = stZ->fGammaECuts[ic];
304         if (elEnergy>gamCut) {                    303         if (elEnergy>gamCut) {
305           // find lower kappa index; compute t    304           // find lower kappa index; compute the 'xi' i.e. cummulative value for
306           // gamCut/elEnergy                      305           // gamCut/elEnergy
307           const G4double cutKappa = std::max(1    306           const G4double cutKappa = std::max(1.e-12, gamCut/elEnergy);
308           const std::size_t iKLow = (cutKappa> << 307           const G4int       iKLow = (cutKappa>1.e-12)
309           ? std::lower_bound(fKappaVect.cbegin << 308           ? std::lower_bound(fKappaVect.begin(), fKappaVect.end(), cutKappa)
310             - fKappaVect.cbegin() -1           << 309             - fKappaVect.begin() -1
311           : 0;                                    310           : 0;
312           const STPoint* stpL = &(stZ->fTables    311           const STPoint* stpL = &(stZ->fTablesPerEnergy[iee]->fSTable[iKLow]);
313           const STPoint* stpH = &(stZ->fTables    312           const STPoint* stpH = &(stZ->fTablesPerEnergy[iee]->fSTable[iKLow+1]);
314           const G4double pA   = stpL->fParA;      313           const G4double pA   = stpL->fParA;
315           const G4double pB   = stpL->fParB;      314           const G4double pB   = stpL->fParB;
316           const G4double etaL = stpL->fCum;       315           const G4double etaL = stpL->fCum;
317           const G4double etaH = stpH->fCum;       316           const G4double etaH = stpH->fCum;
318           const G4double alph = G4Log(cutKappa    317           const G4double alph = G4Log(cutKappa/fKappaVect[iKLow])
319                                /G4Log(fKappaVe    318                                /G4Log(fKappaVect[iKLow+1]/fKappaVect[iKLow]);
320           const G4double dum  = pA*(alph-1.)-1    319           const G4double dum  = pA*(alph-1.)-1.-pB;
321           G4double val = etaL;                    320           G4double val = etaL;
322           if (alph==0.) {                         321           if (alph==0.) {
323             stZ->fTablesPerEnergy[iee]->fCumCu    322             stZ->fTablesPerEnergy[iee]->fCumCutValues[ic] = val;
324           } else {                                323           } else {
325             val = -(dum+std::sqrt(dum*dum-4.*p    324             val = -(dum+std::sqrt(dum*dum-4.*pB*alph*alph))/(2.*pB*alph);
326             val = val*(etaH-etaL)+etaL;           325             val = val*(etaH-etaL)+etaL;
327             stZ->fTablesPerEnergy[iee]->fCumCu    326             stZ->fTablesPerEnergy[iee]->fCumCutValues[ic] = val;
328           }                                       327           }
329         }                                         328         }
330       }                                           329       }
331     }                                             330     }
332   }                                               331   }
333 }                                                 332 }
334                                                   333 
335 // should be called only from LoadSamplingTabl    334 // should be called only from LoadSamplingTables(G4int) and once
336 void G4SBBremTable::LoadSTGrid() {                335 void G4SBBremTable::LoadSTGrid() {
337   const G4String fname =  G4EmParameters::Inst << 336   char* path = std::getenv("G4LEDATA");
                                                   >> 337   if (!path) {
                                                   >> 338     G4Exception("G4SBBremTable::LoadSTGrid()","em0006",
                                                   >> 339                 FatalException, "Environment variable G4LEDATA not defined");
                                                   >> 340     return;
                                                   >> 341   }
                                                   >> 342   const G4String fname =  G4String(path) + "/brem_SB/SBTables/grid";
338   std::ifstream infile(fname,std::ios::in);       343   std::ifstream infile(fname,std::ios::in);
339   if (!infile.is_open()) {                        344   if (!infile.is_open()) {
340     G4String msgc = "Cannot open file: " + fna    345     G4String msgc = "Cannot open file: " + fname;
341     G4Exception("G4SBBremTable::LoadSTGrid()",    346     G4Exception("G4SBBremTable::LoadSTGrid()","em0006",
342                 FatalException, msgc.c_str());    347                 FatalException, msgc.c_str());
343     return;                                       348     return;
344   }                                               349   }
345   // get max Z, # electron energies and # kapp    350   // get max Z, # electron energies and # kappa values
346   infile >> fMaxZet;                              351   infile >> fMaxZet;
347   infile >> fNumElEnergy;                         352   infile >> fNumElEnergy;
348   infile >> fNumKappa;                            353   infile >> fNumKappa;
349   // allocate space for the data and get them     354   // allocate space for the data and get them in:
350   // (1.) first eletron energy grid               355   // (1.) first eletron energy grid
351   fElEnergyVect.resize(fNumElEnergy);             356   fElEnergyVect.resize(fNumElEnergy);
352   fLElEnergyVect.resize(fNumElEnergy);            357   fLElEnergyVect.resize(fNumElEnergy);
353   for (G4int iee=0; iee<fNumElEnergy; ++iee) {    358   for (G4int iee=0; iee<fNumElEnergy; ++iee) {
354     G4double  dum;                                359     G4double  dum;
355     infile >> dum;                                360     infile >> dum;
356     fElEnergyVect[iee]  = dum*CLHEP::MeV;         361     fElEnergyVect[iee]  = dum*CLHEP::MeV;
357     fLElEnergyVect[iee] = G4Log(fElEnergyVect[    362     fLElEnergyVect[iee] = G4Log(fElEnergyVect[iee]);
358   }                                               363   }
359   // (2.) then the kappa grid                     364   // (2.) then the kappa grid
360   fKappaVect.resize(fNumKappa);                   365   fKappaVect.resize(fNumKappa);
361   fLKappaVect.resize(fNumKappa);                  366   fLKappaVect.resize(fNumKappa);
362   for (G4int ik=0; ik<fNumKappa; ++ik) {          367   for (G4int ik=0; ik<fNumKappa; ++ik) {
363     infile >> fKappaVect[ik];                     368     infile >> fKappaVect[ik];
364     fLKappaVect[ik] = G4Log(fKappaVect[ik]);      369     fLKappaVect[ik] = G4Log(fKappaVect[ik]);
365   }                                               370   }
366   // (3.) set size of the main container for s    371   // (3.) set size of the main container for sampling tables
367   fSBSamplingTables.resize(fMaxZet+1,nullptr);    372   fSBSamplingTables.resize(fMaxZet+1,nullptr);
368   // init electron energy grid related variabl    373   // init electron energy grid related variables: use accurate values !!!
369 //  fLogMinElEnergy   = G4Log(fElEnergyVect[0]    374 //  fLogMinElEnergy   = G4Log(fElEnergyVect[0]);
370 //  fILDeltaElEnergy  = 1./G4Log(fElEnergyVect    375 //  fILDeltaElEnergy  = 1./G4Log(fElEnergyVect[1]/fElEnergyVect[0]);
371   const G4double elEmin  = 100.0*CLHEP::eV; //    376   const G4double elEmin  = 100.0*CLHEP::eV; //fElEnergyVect[0];
372   const G4double elEmax  =  10.0*CLHEP::GeV;//    377   const G4double elEmax  =  10.0*CLHEP::GeV;//fElEnergyVect[fNumElEnergy-1];
373   fLogMinElEnergy  = G4Log(elEmin);               378   fLogMinElEnergy  = G4Log(elEmin);
374   fILDeltaElEnergy = 1./(G4Log(elEmax/elEmin)/    379   fILDeltaElEnergy = 1./(G4Log(elEmax/elEmin)/(fNumElEnergy-1.0));
375   // reset min/max energies if needed             380   // reset min/max energies if needed
376   fUsedLowEenergy  = std::max(fUsedLowEenergy     381   fUsedLowEenergy  = std::max(fUsedLowEenergy ,elEmin);
377   fUsedHighEenergy = std::min(fUsedHighEenergy    382   fUsedHighEenergy = std::min(fUsedHighEenergy,elEmax);
378   //                                              383   //
379   infile.close();                                 384   infile.close();
380 }                                                 385 }
381                                                   386 
382 void G4SBBremTable::LoadSamplingTables(G4int i    387 void G4SBBremTable::LoadSamplingTables(G4int iz) {
383   // check if grid needs to be loaded first       388   // check if grid needs to be loaded first
384   if (fMaxZet<0) {                                389   if (fMaxZet<0) {
385     LoadSTGrid();                                 390     LoadSTGrid();
386   }                                               391   }
387   // load data for a given Z only once            392   // load data for a given Z only once
388   iz = std::max(std::min(fMaxZet, iz),1);         393   iz = std::max(std::min(fMaxZet, iz),1);
389                                                << 394   char* path = std::getenv("G4LEDATA");
390   const G4String fname = G4EmParameters::Insta << 395   if (!path) {
                                                   >> 396     G4Exception("G4SBBremTable::LoadSamplingTables()","em0006",
                                                   >> 397                 FatalException, "Environment variable G4LEDATA not defined");
                                                   >> 398     return;
                                                   >> 399   }
                                                   >> 400   const G4String fname =  G4String(path) + "/brem_SB/SBTables/sTableSB_"
391                         + std::to_string(iz);     401                         + std::to_string(iz);
392   std::istringstream infile(std::ios::in);        402   std::istringstream infile(std::ios::in);
393   // read the compressed data file into the st    403   // read the compressed data file into the stream
394   ReadCompressedFile(fname, infile);              404   ReadCompressedFile(fname, infile);
395   // the SamplingTablePerZ object was already     405   // the SamplingTablePerZ object was already created, set size of containers
396   // then load sampling table data for each el    406   // then load sampling table data for each electron energies
397   SamplingTablePerZ* zTable = fSBSamplingTable    407   SamplingTablePerZ* zTable = fSBSamplingTables[iz];
398   //                                              408   //
399   // Determine min/max elektron kinetic energi    409   // Determine min/max elektron kinetic energies and indices
400   const G4double minGammaCut = zTable->fGammaE    410   const G4double minGammaCut = zTable->fGammaECuts[ std::min_element(
401                  std::cbegin(zTable->fGammaECu << 411                  std::begin(zTable->fGammaECuts),std::end(zTable->fGammaECuts))
402                 -std::cbegin(zTable->fGammaECu << 412                 -std::begin(zTable->fGammaECuts)];
403   const G4double elEmin = std::max(fUsedLowEen    413   const G4double elEmin = std::max(fUsedLowEenergy, minGammaCut);
404   const G4double elEmax = fUsedHighEenergy;       414   const G4double elEmax = fUsedHighEenergy;
405   // find low/high elecrton energy indices whe    415   // find low/high elecrton energy indices where tables will be needed
406   // low:                                         416   // low:
407   zTable->fMinElEnergyIndx = 0;                   417   zTable->fMinElEnergyIndx = 0;
408   if (elEmin>=fElEnergyVect[fNumElEnergy-1]) {    418   if (elEmin>=fElEnergyVect[fNumElEnergy-1]) {
409     zTable->fMinElEnergyIndx = fNumElEnergy-1;    419     zTable->fMinElEnergyIndx = fNumElEnergy-1;
410   } else {                                        420   } else {
411     zTable->fMinElEnergyIndx = G4int(std::lowe << 421     zTable->fMinElEnergyIndx = std::lower_bound(fElEnergyVect.begin(),
412                                                << 422                         fElEnergyVect.end(), elEmin) - fElEnergyVect.begin() -1;
413                                      - fElEner << 
414   }                                               423   }
415   // high:                                        424   // high:
416   zTable->fMaxElEnergyIndx = 0;                   425   zTable->fMaxElEnergyIndx = 0;
417   if (elEmax>=fElEnergyVect[fNumElEnergy-1]) {    426   if (elEmax>=fElEnergyVect[fNumElEnergy-1]) {
418     zTable->fMaxElEnergyIndx = fNumElEnergy-1;    427     zTable->fMaxElEnergyIndx = fNumElEnergy-1;
419   } else {                                        428   } else {
420     // lower + 1                                  429     // lower + 1
421     zTable->fMaxElEnergyIndx = G4int(std::lowe << 430     zTable->fMaxElEnergyIndx = std::lower_bound(fElEnergyVect.begin(),
422                                                << 431                            fElEnergyVect.end(), elEmax) - fElEnergyVect.begin();
423                                      - fElEner << 
424   }                                               432   }
425   // protect                                      433   // protect
426   if (zTable->fMaxElEnergyIndx<=zTable->fMinEl    434   if (zTable->fMaxElEnergyIndx<=zTable->fMinElEnergyIndx) {
427     return;                                       435     return;
428   }                                               436   }
429   // load sampling tables that are needed: fil    437   // load sampling tables that are needed: file is already in the stream
430   zTable->fTablesPerEnergy.resize(fNumElEnergy    438   zTable->fTablesPerEnergy.resize(fNumElEnergy, nullptr);
431   for (G4int iee=0; iee<fNumElEnergy; ++iee) {    439   for (G4int iee=0; iee<fNumElEnergy; ++iee) {
432     // go over data that are not needed           440     // go over data that are not needed
433     if (iee<zTable->fMinElEnergyIndx || iee>zT    441     if (iee<zTable->fMinElEnergyIndx || iee>zTable->fMaxElEnergyIndx) {
434       for (G4int ik=0; ik<fNumKappa; ++ik) {      442       for (G4int ik=0; ik<fNumKappa; ++ik) {
435         G4double dum;                             443         G4double dum;
436         infile >> dum; infile >> dum; infile >    444         infile >> dum; infile >> dum; infile >> dum;
437       }                                           445       }
438     } else { // load data that are needed         446     } else { // load data that are needed
439       zTable->fTablesPerEnergy[iee] = new STab    447       zTable->fTablesPerEnergy[iee] = new STable();
440       zTable->fTablesPerEnergy[iee]->fSTable.r    448       zTable->fTablesPerEnergy[iee]->fSTable.resize(fNumKappa);
441       for (G4int ik=0; ik<fNumKappa; ++ik) {      449       for (G4int ik=0; ik<fNumKappa; ++ik) {
442         STPoint &stP = zTable->fTablesPerEnerg    450         STPoint &stP = zTable->fTablesPerEnergy[iee]->fSTable[ik];
443         infile >> stP.fCum;                       451         infile >> stP.fCum;
444         infile >> stP.fParA;                      452         infile >> stP.fParA;
445         infile >> stP.fParB;                      453         infile >> stP.fParB;
446       }                                           454       }
447     }                                             455     }
448   }                                               456   }
449 }                                                 457 }
450                                                   458 
451 // clean away all sampling tables and make rea    459 // clean away all sampling tables and make ready to re-install
452 void G4SBBremTable::ClearSamplingTables() {       460 void G4SBBremTable::ClearSamplingTables() {
453   for (G4int iz=0; iz<fMaxZet+1; ++iz) {          461   for (G4int iz=0; iz<fMaxZet+1; ++iz) {
454     if (fSBSamplingTables[iz]) {                  462     if (fSBSamplingTables[iz]) {
455       for (G4int iee=0; iee<fNumElEnergy; ++ie    463       for (G4int iee=0; iee<fNumElEnergy; ++iee) {
456         if (fSBSamplingTables[iz]->fTablesPerE    464         if (fSBSamplingTables[iz]->fTablesPerEnergy[iee]) {
457           fSBSamplingTables[iz]->fTablesPerEne    465           fSBSamplingTables[iz]->fTablesPerEnergy[iee]->fSTable.clear();
458           fSBSamplingTables[iz]->fTablesPerEne    466           fSBSamplingTables[iz]->fTablesPerEnergy[iee]->fCumCutValues.clear();
459         }                                         467         }
460       }                                           468       }
461       fSBSamplingTables[iz]->fTablesPerEnergy.    469       fSBSamplingTables[iz]->fTablesPerEnergy.clear();
462       fSBSamplingTables[iz]->fGammaECuts.clear    470       fSBSamplingTables[iz]->fGammaECuts.clear();
463       fSBSamplingTables[iz]->fLogGammaECuts.cl    471       fSBSamplingTables[iz]->fLogGammaECuts.clear();
464       fSBSamplingTables[iz]->fMatCutIndxToGamC    472       fSBSamplingTables[iz]->fMatCutIndxToGamCutIndx.clear();
465       //                                          473       //
466       delete fSBSamplingTables[iz];               474       delete fSBSamplingTables[iz];
467       fSBSamplingTables[iz] = nullptr;            475       fSBSamplingTables[iz] = nullptr;
468     }                                             476     }
469   }                                               477   }
470   fSBSamplingTables.clear();                      478   fSBSamplingTables.clear();
471   fElEnergyVect.clear();                          479   fElEnergyVect.clear();
472   fLElEnergyVect.clear();                         480   fLElEnergyVect.clear();
473   fKappaVect.clear();                             481   fKappaVect.clear();
474   fLKappaVect.clear();                            482   fLKappaVect.clear();
475   fMaxZet = -1;                                   483   fMaxZet = -1;
476 }                                                 484 }
477                                                   485 
478 //void G4SBBremTable::Dump() {                    486 //void G4SBBremTable::Dump() {
479 //  G4cerr<< "\n  =====   Dumping ===== \n" <<    487 //  G4cerr<< "\n  =====   Dumping ===== \n" << G4endl;
480 //  for (G4int iz=0; iz<fMaxZet+1; ++iz) {        488 //  for (G4int iz=0; iz<fMaxZet+1; ++iz) {
481 //    if (fSBSamplingTables[iz]) {                489 //    if (fSBSamplingTables[iz]) {
482 //      G4cerr<< "   ----> There are " << fSBS    490 //      G4cerr<< "   ----> There are " << fSBSamplingTables[iz]->fNumGammaCuts
483 //            << " g-cut for Z = " << iz << G4    491 //            << " g-cut for Z = " << iz << G4endl;
484 //      for (std::size_t ic=0; ic<fSBSamplingT << 492 //      for (size_t ic=0; ic<fSBSamplingTables[iz]->fGammaECuts.size(); ++ic)
485 //        G4cerr<< "        i = " << ic << "      493 //        G4cerr<< "        i = " << ic << "  "
486 //              << fSBSamplingTables[iz]->fGam    494 //              << fSBSamplingTables[iz]->fGammaECuts[ic] << G4endl;
487 //    }                                           495 //    }
488 //  }                                             496 //  }
489 //}                                               497 //}
490                                                   498 
491 // find lower bin index of value: used in acse    499 // find lower bin index of value: used in acse of CDF values i.e. val in [0,1)
492 // while vector elements in [0,1]                 500 // while vector elements in [0,1]
493 G4int G4SBBremTable::LinSearch(const std::vect    501 G4int G4SBBremTable::LinSearch(const std::vector<STPoint>& vect,
494                                const G4int siz    502                                const G4int size,
495                                const G4double     503                                const G4double val) {
496   G4int i= 0;                                     504   G4int i= 0;
497   while (i + 3 < size) {                          505   while (i + 3 < size) {
498     if (vect [i + 0].fCum > val) return i + 0;    506     if (vect [i + 0].fCum > val) return i + 0;
499     if (vect [i + 1].fCum > val) return i + 1;    507     if (vect [i + 1].fCum > val) return i + 1;
500     if (vect [i + 2].fCum > val) return i + 2;    508     if (vect [i + 2].fCum > val) return i + 2;
501     if (vect [i + 3].fCum > val) return i + 3;    509     if (vect [i + 3].fCum > val) return i + 3;
502     i += 4;                                       510     i += 4;
503   }                                               511   }
504   while (i < size) {                              512   while (i < size) {
505     if (vect [i].fCum > val)                      513     if (vect [i].fCum > val)
506       break;                                      514       break;
507     ++i;                                          515     ++i;
508   }                                               516   }
509   return i;                                       517   return i;
510 }                                                 518 }
511                                                   519 
512 // uncompress one data file into the input str    520 // uncompress one data file into the input string stream
513 void G4SBBremTable::ReadCompressedFile(const G    521 void G4SBBremTable::ReadCompressedFile(const G4String &fname,
514                                        std::is    522                                        std::istringstream &iss) {
515   std::string *dataString = nullptr;              523   std::string *dataString = nullptr;
516   std::string compfilename(fname+".z");           524   std::string compfilename(fname+".z");
517   // create input stream with binary mode oper    525   // create input stream with binary mode operation and positioning at the end
518   // of the file                                  526   // of the file
519   std::ifstream in(compfilename, std::ios::bin    527   std::ifstream in(compfilename, std::ios::binary | std::ios::ate);
520   if (in.good()) {                                528   if (in.good()) {
521      // get current position in the stream (wa    529      // get current position in the stream (was set to the end)
522      std::streamoff fileSize = in.tellg();     << 530      int fileSize = in.tellg();
523      // set current position being the beginni    531      // set current position being the beginning of the stream
524      in.seekg(0,std::ios::beg);                   532      in.seekg(0,std::ios::beg);
525      // create (zlib) byte buffer for the data    533      // create (zlib) byte buffer for the data
526      Bytef *compdata = new Bytef[fileSize];       534      Bytef *compdata = new Bytef[fileSize];
527      while(in) {                                  535      while(in) {
528         in.read((char*)compdata, fileSize);       536         in.read((char*)compdata, fileSize);
529      }                                            537      }
530      // create (zlib) byte buffer for the unco    538      // create (zlib) byte buffer for the uncompressed data
531      uLongf complen    = (uLongf)(fileSize*4);    539      uLongf complen    = (uLongf)(fileSize*4);
532      Bytef *uncompdata = new Bytef[complen];      540      Bytef *uncompdata = new Bytef[complen];
533      while (Z_OK!=uncompress(uncompdata, &comp    541      while (Z_OK!=uncompress(uncompdata, &complen, compdata, fileSize)) {
534         // increase uncompressed byte buffer      542         // increase uncompressed byte buffer
535         delete[] uncompdata;                      543         delete[] uncompdata;
536         complen   *= 2;                           544         complen   *= 2;
537         uncompdata = new Bytef[complen];          545         uncompdata = new Bytef[complen];
538      }                                            546      }
539      // delete the compressed data buffer         547      // delete the compressed data buffer
540      delete [] compdata;                          548      delete [] compdata;
541      // create a string from the uncompressed     549      // create a string from the uncompressed data (will be deleted by the caller)
542      dataString = new std::string((char*)uncom    550      dataString = new std::string((char*)uncompdata, (long)complen);
543      // delete the uncompressed data buffer       551      // delete the uncompressed data buffer
544      delete [] uncompdata;                        552      delete [] uncompdata;
545   } else {                                        553   } else {
546     std::string msg = "  Problem while trying     554     std::string msg = "  Problem while trying to read "
547                       + compfilename + " data     555                       + compfilename + " data file.\n";
548     G4Exception("G4SBBremTable::ReadCompressed    556     G4Exception("G4SBBremTable::ReadCompressedFile","em0006",
549                 FatalException,msg.c_str());      557                 FatalException,msg.c_str());
550     return;                                       558     return;
551   }                                               559   }
552   // create the input string stream from the d    560   // create the input string stream from the data string
553   if (dataString) {                               561   if (dataString) {
554     iss.str(*dataString);                         562     iss.str(*dataString);
555     in.close();                                   563     in.close();
556     delete dataString;                            564     delete dataString;
557   }                                               565   }
558 }                                                 566 }
559                                                   567