| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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
66 fUsedHighEenergy(-1.), fLogMinElEnergy(-1.)
67 {}
68
69 G4SBBremTable::~G4SBBremTable()
70 {
71 ClearSamplingTables();
72 }
73
74 void G4SBBremTable::Initialize(const G4double
75 {
76 fUsedLowEenergy = lowe;
77 fUsedHighEenergy = highe;
78 BuildSamplingTables();
79 InitSamplingTables();
80 // Dump();
81 }
82
83 // run-time method that samples energy transfe
84 double G4SBBremTable::SampleEnergyTransfer(con
85 con
86 con
87 con
88 con
89 con
90 con
91 {
92 static const G4double kAlpha2Pi = CLHEP::two
93 const G4double zet = (G4double)iZet;
94 const G4int izet = std::max(std::min(fMaxZ
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 = fSBSamplingTa
100 // get the gamma cut of this Z that correspo
101 const std::size_t gamCutIndx = stZ->fMatCutI
102 // gcut was not found: should never happen (
103 if (gamCutIndx >= stZ->fNumGammaCuts || stZ-
104 G4String msg = " Gamma cut="+std::to_strin
105 msg += "in case of Z = " + std::to_string(
106 G4Exception("G4SBBremTable::SampleEnergyTr
107 msg.c_str());
108 }
109 const G4double lGCut = stZ->fLogGammaECuts[g
110 // get the random engine
111 CLHEP::HepRandomEngine* rndmEngine = G4Rando
112 // find lower e- energy bin
113 G4bool isCorner = false; // indicate that th
114 G4bool isSimply = false; // simply sampling:
115 G4int elEnergyIndx = stZ->fMaxElEnergyIndx
116 // only if e- ekin is below the maximum valu
117 if (eekin < fElEnergyVect[elEnergyIndx]) {
118 const G4double val = (lElEnergy-fLogMinElE
119 elEnergyIndx = (G4int)val;
120 G4double pIndxH = val-elEnergyIndx;
121 // check if we are at limiting case: lower
122 if (fElEnergyVect[elEnergyIndx]<=gcut) {
123 // recompute the probability of taking t
124 pIndxH = (lElEnergy-lGCut)/(fLElEnergy
125 isCorner = true;
126 }
127 //
128 if (rndmEngine->flat()<pIndxH) {
129 ++elEnergyIndx; // take the table a
130 } else if (isCorner) { // take the table a
131 // special sampling need to be done if l
132 // actually, we "sample" from a table "b
133 isSimply = true;
134 }
135 }
136 // should never happen under normal conditio
137 if (!stZ->fTablesPerEnergy[elEnergyIndx]) {
138 return 0.;
139 }
140 // Do the photon energy sampling:
141 const STable *st = stZ->fTablesPerEnergy[el
142 const std::vector<G4double>& cVect = st->fCu
143 const std::vector<STPoint>& pVect = st->fST
144 const G4double minVal = cVect[gamCutIndx];
145
146 // should never happen under normal conditio
147 if (minVal >= 1.) {
148 return 0.;
149 }
150 // some transfomrmtion variables used in the
151 const G4double lCurKappaC = lGCut - leekin;
152 const G4double lUsedKappaC = lGCut - fLElEne
153 // dielectric (always) and e+ correction sup
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.-minVa
162 // find lower index of the values in the
163 // instead of binary search because it's
164 const G4int cumLIndx = LinSearch(pVect,
165 // const G4int cumLIndx = std::lower_boun
166 // [](const
167 // return p
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]
173 const G4double lKL = fLKappaVect[cumLI
174 const G4double lKH = fLKappaVect[cumLI
175 const G4double dm1 = (cumRV-cumL)/(cum
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
180 // transform lKappa to [log(gcut/eekin),
181 kappa = G4Exp(lKappa*lCurKappaC/lUsedKa
182 } else {
183 // const G4double upLimit = std::min(1.*C
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+[g
191 suppression = 1./(1.+dielSupConst*invEGamm
192 // add positron correction if particle is
193 if (!isElectron) {
194 const G4double e1 = eekin - gcut;
195 const G4double iBeta1 = (e1 + CLHEP::el
196 / std::sqrt(e1*(
197 const G4double e2 = eekin - eGamma;
198 const G4double iBeta2 = (e2 + CLHEP::el
199 / std::sqrt(e2*(
200 const G4double dum = kAlpha2Pi*zet*(i
201 suppression = (dum > -12.) ? suppression
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 struct
215 // loop over all material-cuts and add gamma
216 const G4ProductionCutsTable
217 *thePCTable = G4ProductionCutsTable::GetProd
218 // a temporary vector to store one element
219 std::vector<std::size_t> vtmp(1,0);
220 std::size_t numMatCuts = thePCTable->GetTabl
221 for (G4int imc=0; imc<(G4int)numMatCuts; ++i
222 const G4MaterialCutsCouple *matCut = thePC
223 if (!matCut->IsUsed()) {
224 continue;
225 }
226 const G4Material* mat = matCut->
227 const G4ElementVector* elemVect = mat->Get
228 const G4int indxMC = matCut->
229 const G4double gamCut = (*(thePCTable->Get
230 const std::size_t numElems = elemVect->siz
231 for (std::size_t ielem=0; ielem<numElems;
232 const G4Element *elem = (*elemVect)[iele
233 const G4int izet = std::max(std::min(fMa
234 if (!fSBSamplingTables[izet]) {
235 // create data structure but do not lo
236 // loaded after we know what are the m
237 // will be needed during the simulatio
238 // LoadSamplingTables(izet);
239 fSBSamplingTables[izet] = new Sampling
240 }
241 // add current gamma cut to the list of
242 // cut value is still not tehre)
243 const std::vector<double> &cVect = fSBSa
244 std::size_t indx = std::find(cVect.cbegi
245 if (indx==cVect.size()) {
246 vtmp[0] = imc;
247 fSBSamplingTables[izet]->fGamCutIndxTo
248 fSBSamplingTables[izet]->fGammaECuts.p
249 fSBSamplingTables[izet]->fLogGammaECut
250 ++fSBSamplingTables[izet]->fNumGammaCu
251 } else {
252 fSBSamplingTables[izet]->fGamCutIndxTo
253 }
254 }
255 }
256 }
257
258 void G4SBBremTable::InitSamplingTables() {
259 const std::size_t numMatCuts = G4ProductionC
260 ->GetTableSize();
261 for (G4int iz=1; iz<fMaxZet+1; ++iz) {
262 SamplingTablePerZ* stZ = fSBSamplingTables
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[
271 // 1 indicates that gamma production is
272 stZ->fTablesPerEnergy[iee]->fCumCutValue
273 // sort gamma cuts and other members acc
274 for (std::size_t i=0; i<stZ->fNumGammaCu
275 for (std::size_t j=i+1; j<stZ->fNumGam
276 if (stZ->fGammaECuts[j]<stZ->fGammaE
277 G4double dum0 =
278 G4double dum1 =
279 std::vector<std::size_t> dumv =
280 stZ->fGammaECuts[i] =
281 stZ->fLogGammaECuts[i] =
282 stZ->fGamCutIndxToMatCutIndx[i] =
283 stZ->fGammaECuts[j] =
284 stZ->fLogGammaECuts[j] =
285 stZ->fGamCutIndxToMatCutIndx[j] =
286 }
287 }
288 }
289 // set couple indices to store the corre
290 stZ->fMatCutIndxToGamCutIndx.resize(numM
291 for (std::size_t i=0; i<stZ->fGamCutIndx
292 for (std::size_t j=0; j<stZ->fGamCutIn
293 stZ->fMatCutIndxToGamCutIndx[stZ->fG
294 }
295 }
296 // clear temporary vector
297 for (std::size_t i=0; i<stZ->fGamCutIndx
298 stZ->fGamCutIndxToMatCutIndx[i].clear(
299 }
300 stZ->fGamCutIndxToMatCutIndx.clear();
301 // init for each gamma cut that are bel
302 for (std::size_t ic=0; ic<stZ->fNumGamma
303 const G4double gamCut = stZ->fGammaECu
304 if (elEnergy>gamCut) {
305 // find lower kappa index; compute t
306 // gamCut/elEnergy
307 const G4double cutKappa = std::max(1
308 const std::size_t iKLow = (cutKappa>
309 ? std::lower_bound(fKappaVect.cbegin
310 - fKappaVect.cbegin() -1
311 : 0;
312 const STPoint* stpL = &(stZ->fTables
313 const STPoint* stpH = &(stZ->fTables
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
319 /G4Log(fKappaVe
320 const G4double dum = pA*(alph-1.)-1
321 G4double val = etaL;
322 if (alph==0.) {
323 stZ->fTablesPerEnergy[iee]->fCumCu
324 } else {
325 val = -(dum+std::sqrt(dum*dum-4.*p
326 val = val*(etaH-etaL)+etaL;
327 stZ->fTablesPerEnergy[iee]->fCumCu
328 }
329 }
330 }
331 }
332 }
333 }
334
335 // should be called only from LoadSamplingTabl
336 void G4SBBremTable::LoadSTGrid() {
337 const G4String fname = G4EmParameters::Inst
338 std::ifstream infile(fname,std::ios::in);
339 if (!infile.is_open()) {
340 G4String msgc = "Cannot open file: " + fna
341 G4Exception("G4SBBremTable::LoadSTGrid()",
342 FatalException, msgc.c_str());
343 return;
344 }
345 // get max Z, # electron energies and # kapp
346 infile >> fMaxZet;
347 infile >> fNumElEnergy;
348 infile >> fNumKappa;
349 // allocate space for the data and get them
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[
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 s
367 fSBSamplingTables.resize(fMaxZet+1,nullptr);
368 // init electron energy grid related variabl
369 // fLogMinElEnergy = G4Log(fElEnergyVect[0]
370 // fILDeltaElEnergy = 1./G4Log(fElEnergyVect
371 const G4double elEmin = 100.0*CLHEP::eV; //
372 const G4double elEmax = 10.0*CLHEP::GeV;//
373 fLogMinElEnergy = G4Log(elEmin);
374 fILDeltaElEnergy = 1./(G4Log(elEmax/elEmin)/
375 // reset min/max energies if needed
376 fUsedLowEenergy = std::max(fUsedLowEenergy
377 fUsedHighEenergy = std::min(fUsedHighEenergy
378 //
379 infile.close();
380 }
381
382 void G4SBBremTable::LoadSamplingTables(G4int i
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::Insta
391 + std::to_string(iz);
392 std::istringstream infile(std::ios::in);
393 // read the compressed data file into the st
394 ReadCompressedFile(fname, infile);
395 // the SamplingTablePerZ object was already
396 // then load sampling table data for each el
397 SamplingTablePerZ* zTable = fSBSamplingTable
398 //
399 // Determine min/max elektron kinetic energi
400 const G4double minGammaCut = zTable->fGammaE
401 std::cbegin(zTable->fGammaECu
402 -std::cbegin(zTable->fGammaECu
403 const G4double elEmin = std::max(fUsedLowEen
404 const G4double elEmax = fUsedHighEenergy;
405 // find low/high elecrton energy indices whe
406 // low:
407 zTable->fMinElEnergyIndx = 0;
408 if (elEmin>=fElEnergyVect[fNumElEnergy-1]) {
409 zTable->fMinElEnergyIndx = fNumElEnergy-1;
410 } else {
411 zTable->fMinElEnergyIndx = G4int(std::lowe
412
413 - fElEner
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::lowe
422
423 - fElEner
424 }
425 // protect
426 if (zTable->fMaxElEnergyIndx<=zTable->fMinEl
427 return;
428 }
429 // load sampling tables that are needed: fil
430 zTable->fTablesPerEnergy.resize(fNumElEnergy
431 for (G4int iee=0; iee<fNumElEnergy; ++iee) {
432 // go over data that are not needed
433 if (iee<zTable->fMinElEnergyIndx || iee>zT
434 for (G4int ik=0; ik<fNumKappa; ++ik) {
435 G4double dum;
436 infile >> dum; infile >> dum; infile >
437 }
438 } else { // load data that are needed
439 zTable->fTablesPerEnergy[iee] = new STab
440 zTable->fTablesPerEnergy[iee]->fSTable.r
441 for (G4int ik=0; ik<fNumKappa; ++ik) {
442 STPoint &stP = zTable->fTablesPerEnerg
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 rea
452 void G4SBBremTable::ClearSamplingTables() {
453 for (G4int iz=0; iz<fMaxZet+1; ++iz) {
454 if (fSBSamplingTables[iz]) {
455 for (G4int iee=0; iee<fNumElEnergy; ++ie
456 if (fSBSamplingTables[iz]->fTablesPerE
457 fSBSamplingTables[iz]->fTablesPerEne
458 fSBSamplingTables[iz]->fTablesPerEne
459 }
460 }
461 fSBSamplingTables[iz]->fTablesPerEnergy.
462 fSBSamplingTables[iz]->fGammaECuts.clear
463 fSBSamplingTables[iz]->fLogGammaECuts.cl
464 fSBSamplingTables[iz]->fMatCutIndxToGamC
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" <<
480 // for (G4int iz=0; iz<fMaxZet+1; ++iz) {
481 // if (fSBSamplingTables[iz]) {
482 // G4cerr<< " ----> There are " << fSBS
483 // << " g-cut for Z = " << iz << G4
484 // for (std::size_t ic=0; ic<fSBSamplingT
485 // G4cerr<< " i = " << ic << "
486 // << fSBSamplingTables[iz]->fGam
487 // }
488 // }
489 //}
490
491 // find lower bin index of value: used in acse
492 // while vector elements in [0,1]
493 G4int G4SBBremTable::LinSearch(const std::vect
494 const G4int siz
495 const G4double
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 str
513 void G4SBBremTable::ReadCompressedFile(const G
514 std::is
515 std::string *dataString = nullptr;
516 std::string compfilename(fname+".z");
517 // create input stream with binary mode oper
518 // of the file
519 std::ifstream in(compfilename, std::ios::bin
520 if (in.good()) {
521 // get current position in the stream (wa
522 std::streamoff fileSize = in.tellg();
523 // set current position being the beginni
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 unco
531 uLongf complen = (uLongf)(fileSize*4);
532 Bytef *uncompdata = new Bytef[complen];
533 while (Z_OK!=uncompress(uncompdata, &comp
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
542 dataString = new std::string((char*)uncom
543 // delete the uncompressed data buffer
544 delete [] uncompdata;
545 } else {
546 std::string msg = " Problem while trying
547 + compfilename + " data
548 G4Exception("G4SBBremTable::ReadCompressed
549 FatalException,msg.c_str());
550 return;
551 }
552 // create the input string stream from the d
553 if (dataString) {
554 iss.str(*dataString);
555 in.close();
556 delete dataString;
557 }
558 }
559