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