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.1.2)


  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 std::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<std::size_t> vtmp(1,0);
220   std::size_t numMatCuts = thePCTable->GetTabl    219   std::size_t numMatCuts = thePCTable->GetTableSize();
221   for (G4int imc=0; imc<(G4int)numMatCuts; ++i    220   for (G4int imc=0; imc<(G4int)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 std::size_t numElems = elemVect->size();
231     for (std::size_t ielem=0; ielem<numElems;     230     for (std::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       std::size_t indx = std::find(cVect.cbegin(), cVect.cend(), gamCut)-cVect.cbegin();
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 std::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 (std::size_t i=0; i<stZ->fNumGammaCuts-1; ++i) {
275         for (std::size_t j=i+1; j<stZ->fNumGam    274         for (std::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<std::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 (std::size_t i=0; i<stZ->fGamCutIndxToMatCutIndx.size(); ++i) {
292         for (std::size_t j=0; j<stZ->fGamCutIn    291         for (std::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 (std::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 (std::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 std::size_t iKLow = (cutKappa>1.e-12)
309           ? std::lower_bound(fKappaVect.cbegin    308           ? std::lower_bound(fKappaVect.cbegin(), fKappaVect.cend(), cutKappa)
310             - fKappaVect.cbegin() -1              309             - fKappaVect.cbegin() -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   const char* path = G4FindDataDir("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   const char* path = G4FindDataDir("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::cbegin(zTable->fGammaECuts),std::cend(zTable->fGammaECuts))
402                 -std::cbegin(zTable->fGammaECu    412                 -std::cbegin(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 = G4int(std::lower_bound(fElEnergyVect.cbegin(),
412                                                   422                                                 fElEnergyVect.cend(), elEmin)
413                                      - fElEner    423                                      - fElEnergyVect.cbegin() -1);
414   }                                               424   }
415   // high:                                        425   // high:
416   zTable->fMaxElEnergyIndx = 0;                   426   zTable->fMaxElEnergyIndx = 0;
417   if (elEmax>=fElEnergyVect[fNumElEnergy-1]) {    427   if (elEmax>=fElEnergyVect[fNumElEnergy-1]) {
418     zTable->fMaxElEnergyIndx = fNumElEnergy-1;    428     zTable->fMaxElEnergyIndx = fNumElEnergy-1;
419   } else {                                        429   } else {
420     // lower + 1                                  430     // lower + 1
421     zTable->fMaxElEnergyIndx = G4int(std::lowe    431     zTable->fMaxElEnergyIndx = G4int(std::lower_bound(fElEnergyVect.cbegin(),
422                                                   432                                                 fElEnergyVect.cend(), elEmax)
423                                      - fElEner    433                                      - fElEnergyVect.cbegin());
424   }                                               434   }
425   // protect                                      435   // protect
426   if (zTable->fMaxElEnergyIndx<=zTable->fMinEl    436   if (zTable->fMaxElEnergyIndx<=zTable->fMinElEnergyIndx) {
427     return;                                       437     return;
428   }                                               438   }
429   // load sampling tables that are needed: fil    439   // load sampling tables that are needed: file is already in the stream
430   zTable->fTablesPerEnergy.resize(fNumElEnergy    440   zTable->fTablesPerEnergy.resize(fNumElEnergy, nullptr);
431   for (G4int iee=0; iee<fNumElEnergy; ++iee) {    441   for (G4int iee=0; iee<fNumElEnergy; ++iee) {
432     // go over data that are not needed           442     // go over data that are not needed
433     if (iee<zTable->fMinElEnergyIndx || iee>zT    443     if (iee<zTable->fMinElEnergyIndx || iee>zTable->fMaxElEnergyIndx) {
434       for (G4int ik=0; ik<fNumKappa; ++ik) {      444       for (G4int ik=0; ik<fNumKappa; ++ik) {
435         G4double dum;                             445         G4double dum;
436         infile >> dum; infile >> dum; infile >    446         infile >> dum; infile >> dum; infile >> dum;
437       }                                           447       }
438     } else { // load data that are needed         448     } else { // load data that are needed
439       zTable->fTablesPerEnergy[iee] = new STab    449       zTable->fTablesPerEnergy[iee] = new STable();
440       zTable->fTablesPerEnergy[iee]->fSTable.r    450       zTable->fTablesPerEnergy[iee]->fSTable.resize(fNumKappa);
441       for (G4int ik=0; ik<fNumKappa; ++ik) {      451       for (G4int ik=0; ik<fNumKappa; ++ik) {
442         STPoint &stP = zTable->fTablesPerEnerg    452         STPoint &stP = zTable->fTablesPerEnergy[iee]->fSTable[ik];
443         infile >> stP.fCum;                       453         infile >> stP.fCum;
444         infile >> stP.fParA;                      454         infile >> stP.fParA;
445         infile >> stP.fParB;                      455         infile >> stP.fParB;
446       }                                           456       }
447     }                                             457     }
448   }                                               458   }
449 }                                                 459 }
450                                                   460 
451 // clean away all sampling tables and make rea    461 // clean away all sampling tables and make ready to re-install
452 void G4SBBremTable::ClearSamplingTables() {       462 void G4SBBremTable::ClearSamplingTables() {
453   for (G4int iz=0; iz<fMaxZet+1; ++iz) {          463   for (G4int iz=0; iz<fMaxZet+1; ++iz) {
454     if (fSBSamplingTables[iz]) {                  464     if (fSBSamplingTables[iz]) {
455       for (G4int iee=0; iee<fNumElEnergy; ++ie    465       for (G4int iee=0; iee<fNumElEnergy; ++iee) {
456         if (fSBSamplingTables[iz]->fTablesPerE    466         if (fSBSamplingTables[iz]->fTablesPerEnergy[iee]) {
457           fSBSamplingTables[iz]->fTablesPerEne    467           fSBSamplingTables[iz]->fTablesPerEnergy[iee]->fSTable.clear();
458           fSBSamplingTables[iz]->fTablesPerEne    468           fSBSamplingTables[iz]->fTablesPerEnergy[iee]->fCumCutValues.clear();
459         }                                         469         }
460       }                                           470       }
461       fSBSamplingTables[iz]->fTablesPerEnergy.    471       fSBSamplingTables[iz]->fTablesPerEnergy.clear();
462       fSBSamplingTables[iz]->fGammaECuts.clear    472       fSBSamplingTables[iz]->fGammaECuts.clear();
463       fSBSamplingTables[iz]->fLogGammaECuts.cl    473       fSBSamplingTables[iz]->fLogGammaECuts.clear();
464       fSBSamplingTables[iz]->fMatCutIndxToGamC    474       fSBSamplingTables[iz]->fMatCutIndxToGamCutIndx.clear();
465       //                                          475       //
466       delete fSBSamplingTables[iz];               476       delete fSBSamplingTables[iz];
467       fSBSamplingTables[iz] = nullptr;            477       fSBSamplingTables[iz] = nullptr;
468     }                                             478     }
469   }                                               479   }
470   fSBSamplingTables.clear();                      480   fSBSamplingTables.clear();
471   fElEnergyVect.clear();                          481   fElEnergyVect.clear();
472   fLElEnergyVect.clear();                         482   fLElEnergyVect.clear();
473   fKappaVect.clear();                             483   fKappaVect.clear();
474   fLKappaVect.clear();                            484   fLKappaVect.clear();
475   fMaxZet = -1;                                   485   fMaxZet = -1;
476 }                                                 486 }
477                                                   487 
478 //void G4SBBremTable::Dump() {                    488 //void G4SBBremTable::Dump() {
479 //  G4cerr<< "\n  =====   Dumping ===== \n" <<    489 //  G4cerr<< "\n  =====   Dumping ===== \n" << G4endl;
480 //  for (G4int iz=0; iz<fMaxZet+1; ++iz) {        490 //  for (G4int iz=0; iz<fMaxZet+1; ++iz) {
481 //    if (fSBSamplingTables[iz]) {                491 //    if (fSBSamplingTables[iz]) {
482 //      G4cerr<< "   ----> There are " << fSBS    492 //      G4cerr<< "   ----> There are " << fSBSamplingTables[iz]->fNumGammaCuts
483 //            << " g-cut for Z = " << iz << G4    493 //            << " g-cut for Z = " << iz << G4endl;
484 //      for (std::size_t ic=0; ic<fSBSamplingT    494 //      for (std::size_t ic=0; ic<fSBSamplingTables[iz]->fGammaECuts.size(); ++ic)
485 //        G4cerr<< "        i = " << ic << "      495 //        G4cerr<< "        i = " << ic << "  "
486 //              << fSBSamplingTables[iz]->fGam    496 //              << fSBSamplingTables[iz]->fGammaECuts[ic] << G4endl;
487 //    }                                           497 //    }
488 //  }                                             498 //  }
489 //}                                               499 //}
490                                                   500 
491 // find lower bin index of value: used in acse    501 // 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]                 502 // while vector elements in [0,1]
493 G4int G4SBBremTable::LinSearch(const std::vect    503 G4int G4SBBremTable::LinSearch(const std::vector<STPoint>& vect,
494                                const G4int siz    504                                const G4int size,
495                                const G4double     505                                const G4double val) {
496   G4int i= 0;                                     506   G4int i= 0;
497   while (i + 3 < size) {                          507   while (i + 3 < size) {
498     if (vect [i + 0].fCum > val) return i + 0;    508     if (vect [i + 0].fCum > val) return i + 0;
499     if (vect [i + 1].fCum > val) return i + 1;    509     if (vect [i + 1].fCum > val) return i + 1;
500     if (vect [i + 2].fCum > val) return i + 2;    510     if (vect [i + 2].fCum > val) return i + 2;
501     if (vect [i + 3].fCum > val) return i + 3;    511     if (vect [i + 3].fCum > val) return i + 3;
502     i += 4;                                       512     i += 4;
503   }                                               513   }
504   while (i < size) {                              514   while (i < size) {
505     if (vect [i].fCum > val)                      515     if (vect [i].fCum > val)
506       break;                                      516       break;
507     ++i;                                          517     ++i;
508   }                                               518   }
509   return i;                                       519   return i;
510 }                                                 520 }
511                                                   521 
512 // uncompress one data file into the input str    522 // uncompress one data file into the input string stream
513 void G4SBBremTable::ReadCompressedFile(const G    523 void G4SBBremTable::ReadCompressedFile(const G4String &fname,
514                                        std::is    524                                        std::istringstream &iss) {
515   std::string *dataString = nullptr;              525   std::string *dataString = nullptr;
516   std::string compfilename(fname+".z");           526   std::string compfilename(fname+".z");
517   // create input stream with binary mode oper    527   // create input stream with binary mode operation and positioning at the end
518   // of the file                                  528   // of the file
519   std::ifstream in(compfilename, std::ios::bin    529   std::ifstream in(compfilename, std::ios::binary | std::ios::ate);
520   if (in.good()) {                                530   if (in.good()) {
521      // get current position in the stream (wa    531      // get current position in the stream (was set to the end)
522      std::streamoff fileSize = in.tellg();        532      std::streamoff fileSize = in.tellg();
523      // set current position being the beginni    533      // set current position being the beginning of the stream
524      in.seekg(0,std::ios::beg);                   534      in.seekg(0,std::ios::beg);
525      // create (zlib) byte buffer for the data    535      // create (zlib) byte buffer for the data
526      Bytef *compdata = new Bytef[fileSize];       536      Bytef *compdata = new Bytef[fileSize];
527      while(in) {                                  537      while(in) {
528         in.read((char*)compdata, fileSize);       538         in.read((char*)compdata, fileSize);
529      }                                            539      }
530      // create (zlib) byte buffer for the unco    540      // create (zlib) byte buffer for the uncompressed data
531      uLongf complen    = (uLongf)(fileSize*4);    541      uLongf complen    = (uLongf)(fileSize*4);
532      Bytef *uncompdata = new Bytef[complen];      542      Bytef *uncompdata = new Bytef[complen];
533      while (Z_OK!=uncompress(uncompdata, &comp    543      while (Z_OK!=uncompress(uncompdata, &complen, compdata, fileSize)) {
534         // increase uncompressed byte buffer      544         // increase uncompressed byte buffer
535         delete[] uncompdata;                      545         delete[] uncompdata;
536         complen   *= 2;                           546         complen   *= 2;
537         uncompdata = new Bytef[complen];          547         uncompdata = new Bytef[complen];
538      }                                            548      }
539      // delete the compressed data buffer         549      // delete the compressed data buffer
540      delete [] compdata;                          550      delete [] compdata;
541      // create a string from the uncompressed     551      // create a string from the uncompressed data (will be deleted by the caller)
542      dataString = new std::string((char*)uncom    552      dataString = new std::string((char*)uncompdata, (long)complen);
543      // delete the uncompressed data buffer       553      // delete the uncompressed data buffer
544      delete [] uncompdata;                        554      delete [] uncompdata;
545   } else {                                        555   } else {
546     std::string msg = "  Problem while trying     556     std::string msg = "  Problem while trying to read "
547                       + compfilename + " data     557                       + compfilename + " data file.\n";
548     G4Exception("G4SBBremTable::ReadCompressed    558     G4Exception("G4SBBremTable::ReadCompressedFile","em0006",
549                 FatalException,msg.c_str());      559                 FatalException,msg.c_str());
550     return;                                       560     return;
551   }                                               561   }
552   // create the input string stream from the d    562   // create the input string stream from the data string
553   if (dataString) {                               563   if (dataString) {
554     iss.str(*dataString);                         564     iss.str(*dataString);
555     in.close();                                   565     in.close();
556     delete dataString;                            566     delete dataString;
557   }                                               567   }
558 }                                                 568 }
559                                                   569