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