| 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 //
30 // File name: G4GSPWACorrections
31 //
32 // Author: Mihaly Novak
33 //
34 // Creation date: 17.10.2017
35 //
36 // Modifications:
37 // 02.02.2018 M.Novak: fixed initialization of
38 //
39 // Class description: see the header file.
40 //
41 // -------------------------------------------
42
43 #include "G4GSPWACorrections.hh"
44
45 #include "G4PhysicalConstants.hh"
46 #include "G4Log.hh"
47 #include "G4Exp.hh"
48
49 #include "G4ProductionCutsTable.hh"
50 #include "G4MaterialCutsCouple.hh"
51 #include "G4Material.hh"
52 #include "G4ElementVector.hh"
53 #include "G4Element.hh"
54 #include "G4EmParameters.hh"
55
56
57 const std::string G4GSPWACorrections::gElemSym
58 "C" ,"N" ,"O" ,"F" ,"Ne","Na","Mg","Al","Si",
59 "Ti","V" ,"Cr","Mn","Fe","Co","Ni","Cu","Zn",
60 "Sr","Y" ,"Zr","Nb","Mo","Tc","Ru","Rh","Pd",
61 "Xe","Cs","Ba","La","Ce","Pr","Nd","Pm","Sm",
62 "Yb","Lu","Hf","Ta","W" ,"Re","Os","Ir","Pt",
63 "Rn","Fr","Ra","Ac","Th","Pa","U" ,"Np","Pu",
64
65 G4GSPWACorrections::G4GSPWACorrections(G4bool
66 // init grids related data member values
67 fMaxEkin = CLHEP::electron_mass_c2*(1
68 fLogMinEkin = G4Log(gMinEkin);
69 fInvLogDelEkin = (gNumEkin-gNumBeta2)/G4Log
70 G4double pt2 = gMidEkin*(gMidEkin+2.0*CLH
71 fMinBeta2 = pt2/(pt2+CLHEP::electron_m
72 fInvDelBeta2 = (gNumBeta2-1.)/(gMaxBeta2-
73 }
74
75
76 G4GSPWACorrections::~G4GSPWACorrections() {
77 ClearDataPerElement();
78 ClearDataPerMaterial();
79 }
80
81
82 void G4GSPWACorrections::GetPWACorrectionFact
83
84 G4int ekinIndxLow = 0;
85 G4double remRfaction = 0.;
86 if (beta2>=gMaxBeta2) {
87 ekinIndxLow = gNumEkin - 1;
88 // remRfaction = -1.
89 } else if (beta2>=fMinBeta2) { // linear in
90 remRfaction = (beta2 - fMinBeta2) * fInv
91 ekinIndxLow = (G4int)remRfaction;
92 remRfaction -= ekinIndxLow;
93 ekinIndxLow += (gNumEkin - gNumBeta2);
94 } else if (logekin>=fLogMinEkin) {
95 remRfaction = (logekin - fLogMinEkin) *
96 ekinIndxLow = (G4int)remRfaction;
97 remRfaction -= ekinIndxLow;
98 } // the defaults otherwise i.e. use the low
99 //
100 DataPerMaterial *data = fDataPerMaterial[mat
101 corToScr = data->fCorScreening[ekinIndx
102 corToQ1 = data->fCorFirstMoment[ekinIn
103 corToG2PerG1 = data->fCorSecondMoment[ekinI
104 if (remRfaction>0.) {
105 corToScr += remRfaction*(data->fCorSc
106 corToQ1 += remRfaction*(data->fCorFi
107 corToG2PerG1 += remRfaction*(data->fCorSe
108 }
109 }
110
111
112 void G4GSPWACorrections::Initialise() {
113 // load PWA correction data for each element
114 InitDataPerElement();
115 // clear PWA correction data per material
116 ClearDataPerMaterial();
117 // initialise PWA correction data for the ma
118 InitDataPerMaterials();
119 }
120
121
122 void G4GSPWACorrections::InitDataPerElement()
123 // do it only once
124 if (fDataPerElement.size()<gMaxZet+1) {
125 fDataPerElement.resize(gMaxZet+1,nullptr);
126 }
127 // loop over all materials, for those that a
128 // corresponding data has not been loaded ye
129 G4ProductionCutsTable *thePCTable = G4Produc
130 G4int numMatCuts = (G4int)thePCTable->GetTab
131 for (G4int imc=0; imc<numMatCuts; ++imc) {
132 const G4MaterialCutsCouple *matCut = theP
133 if (!matCut->IsUsed()) {
134 continue;
135 }
136 const G4Material *mat = matCut->
137 const G4ElementVector *elemVect = mat->Get
138 //
139 std::size_t numElems = elemVect->size();
140 for (std::size_t ielem=0; ielem<numElems;
141 const G4Element *elem = (*elemVect)[iele
142 G4int izet = G4lrint(elem->GetZ());
143 if (izet>gMaxZet) {
144 izet = gMaxZet;
145 }
146 if (!fDataPerElement[izet]) {
147 LoadDataElement(elem);
148 }
149 }
150 }
151 }
152
153
154 void G4GSPWACorrections::InitDataPerMaterials(
155 // prepare size of the container
156 std::size_t numMaterials = G4Material::GetNu
157 if (fDataPerMaterial.size()!=numMaterials) {
158 fDataPerMaterial.resize(numMaterials);
159 }
160 // init. PWA correction data for the Materia
161 G4ProductionCutsTable *thePCTable = G4Produc
162 G4int numMatCuts = (G4int)thePCTable->GetTab
163 for (G4int imc=0; imc<numMatCuts; ++imc) {
164 const G4MaterialCutsCouple *matCut = theP
165 if (!matCut->IsUsed()) {
166 continue;
167 }
168 const G4Material *mat = matCut->GetMateria
169 if (!fDataPerMaterial[mat->GetIndex()]) {
170 InitDataMaterial(mat);
171 }
172 }
173 }
174
175
176 // it's called only if data has not been loade
177 void G4GSPWACorrections::LoadDataElement(const
178 // allocate memory
179 G4int izet = elem->GetZasInt();
180 if (izet>gMaxZet) {
181 izet = gMaxZet;
182 }
183 // load data from file
184 std::string path = G4EmParameters::Instance(
185 if (fIsElectron) {
186 path += "/msc_GS/PWACor/el/";
187 } else {
188 path += "/msc_GS/PWACor/pos/";
189 }
190 std::string fname = path+"cf_"+gElemSymbols
191 std::ifstream infile(fname,std::ios::in);
192 if (!infile.is_open()) {
193 std::string msg = " Problem while trying
194 G4Exception("G4GSPWACorrection::LoadDataEl
195 return;
196 }
197 // allocate data structure
198 auto perElem = new DataPerMaterial();
199 perElem->fCorScreening.resize(gNumEkin,0.0);
200 perElem->fCorFirstMoment.resize(gNumEkin,0.0
201 perElem->fCorSecondMoment.resize(gNumEkin,0.
202 fDataPerElement[izet] = perElem;
203 G4double dum0;
204 for (G4int iek=0; iek<gNumEkin; ++iek) {
205 infile >> dum0;
206 infile >> perElem->fCorScreening[iek];
207 infile >> perElem->fCorFirstMoment[iek];
208 infile >> perElem->fCorSecondMoment[iek];
209 }
210 infile.close();
211 }
212
213
214 void G4GSPWACorrections::InitDataMaterial(cons
215 constexpr G4double const1 = 7821.6; /
216 constexpr G4double const2 = 0.1569; /
217 constexpr G4double finstrc2 = 5.325135453E-5
218
219 G4double constFactor = CLHEP::electro
220 constFactor *= constFactor;
221 // allocate memory
222 auto perMat = new DataPerMaterial();
223 perMat->fCorScreening.resize(gNumEkin,0.0);
224 perMat->fCorFirstMoment.resize(gNumEkin,0.0)
225 perMat->fCorSecondMoment.resize(gNumEkin,0.0
226 fDataPerMaterial[mat->GetIndex()] = perMat;
227 //
228 const G4ElementVector* elemVect =
229 const G4int numElems =
230 const G4double* nbAtomsPerVolVect =
231 G4double totNbAtomsPerVol =
232 //
233 // 1. Compute material dependent part of Mol
234 // (with \xi=1 (i.e. total sub-threshold
235 G4double moliereBc = 0.0;
236 G4double moliereXc2 = 0.0;
237 G4double zs = 0.0;
238 G4double ze = 0.0;
239 G4double zx = 0.0;
240 G4double sa = 0.0;
241 G4double xi = 1.0;
242 for (G4int ielem=0; ielem<numElems; ++ielem)
243 G4double zet = (*elemVect)[ielem]->GetZ();
244 if (zet>gMaxZet) {
245 zet = (G4double)gMaxZet;
246 }
247 G4double iwa = (*elemVect)[ielem]->GetN()
248 G4double ipz = nbAtomsPerVolVect[ielem]/t
249 G4double dum = ipz*zet*(zet+xi);
250 zs += dum;
251 ze += dum*(-2.0/3.0)*G4Log(zet);
252 zx += dum*G4Log(1.0+3.34*finstrc
253 sa += ipz*iwa;
254 }
255 G4double density = mat->GetDensity()*CLHEP::
256 //
257 G4double z0 = (0.0 == sa) ? 0.0 : zs/sa;
258 G4double z1 = (0.0 == zs) ? 0.0 : (ze - zx)/
259 moliereBc = const1*density*z0*G4Exp(z1); /
260 moliereXc2 = const2*density*z0; // [MeV2/cm
261 // change to Geant4 internal units of 1/leng
262 moliereBc *= 1.0/CLHEP::cm;
263 moliereXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::c
264 //
265 // 2. loop over the kinetic energy grid
266 for (G4int iek=0; iek<gNumEkin; ++iek) {
267 // 2./a. set current kinetic energy and pt
268 G4double ekin = G4Exp(fLogMinEk
269 G4double pt2 = ekin*(ekin+2.0*CLHEP::el
270 if (ekin>gMidEkin) {
271 G4double b2 = fMinBeta2+(iek-(gNumEk
272 ekin = CLHEP::electron_mass_c2*(1./std
273 pt2 = ekin*(ekin+2.0*CLHEP::electron_
274 }
275 // 2./b. loop over the elements at the cur
276 for (G4int ielem=0; ielem<numElems; ++iele
277 const G4Element *elem = (*elemVect)[iele
278 G4double zet = elem->GetZ();
279 if (zet>gMaxZet) {
280 zet = (G4double)gMaxZet;
281 }
282 G4int izet = G4lrint(zet);
283 // loaded PWA corrections for the curren
284 DataPerMaterial *perElem = fDataPerElem
285 //
286 // xi should be one i.e. z(z+1) since to
287 G4double nZZPlus1 = nbAtomsPerVolVect[i
288 G4double Z23 = std::pow(zet,2./3.)
289 //
290 // 2./b./(i) Add the 3 PWA correction fa
291 G4double mcScrCF = perElem->fCorScreenin
292 // compute the screening parameter corre
293 // src_{mc} = C \exp\left[ \frac{ \sum_i
294 // with C = \frac{(mc^2)^\alpha^2} {4(pc
295 // here we compute the \sum_i n_i Z_i(Z_
296 perMat->fCorScreening[iek] += nZZPlus1*G
297 // compute the corrected screening param
298 // src(Z_i)_{mc} = \frac{(mc^2)^\alpha^2
299 mcScrCF *= constFactor*Z23/(4.*pt2);
300 // compute first moment correction facto
301 // q1_{mc} = \frac{ \sum_i n_i Z_i(Z_i+1
302 // where:
303 // A_i(src(Z_i)_{mc}) = [\ln(1+1/src(Z_i
304 // B_i = \beta_i \gamma_i with beta_i(Z_
305 // and \gamma_i = \sigma(Z_i)_{el}^(MC-D
306 // C(src_{mc}) = [\ln(1+1/src_{mc}) - 1/
307 // A_i x B_i is stored in file per e-/e+
308 // here we compute the \sum_i n_i Z_i(Z_
309 perMat->fCorFirstMoment[iek] += nZZPlus1
310 // compute the second moment correction
311 // [G2/G1]_{mc} = \frac{ \sum_i n_i Z_i(
312 // with A_i(Z_i) = G2(Z_i)^{PWA}/G1(Z_i)
313 // here we compute the \sum_i n_i Z_i(Z_
314 perMat->fCorSecondMoment[iek] += nZZPlus
315 //
316 // 2./b./(ii) When the last element has
317 if (ielem==numElems-1) {
318 //
319 // 1. the remaining part of the sreeni
320 // (Moliere screening parameter = m
321 G4double dumScr = G4Exp(perMat->fCor
322 perMat->fCorScreening[iek] = constFact
323 //
324 // 2. the remaining part of the first
325 // screening parameter (= (mc^2)^\a
326 G4double scrCorTed = constFactor*dumSc
327 G4double dum0 = G4Log(1.+1./scrCo
328 perMat->fCorFirstMoment[iek] = perMat-
329 //
330 // 3. the remaining part of the second
331 // screening parameter
332 G4double G2PerG1 = 3.*(1.+scrCorTed
333 perMat->fCorSecondMoment[iek] = perMat
334 }
335 }
336 }
337 }
338
339
340
341 void G4GSPWACorrections::ClearDataPerElement()
342 for (std::size_t i=0; i<fDataPerElement.size
343 if (fDataPerElement[i]) {
344 fDataPerElement[i]->fCorScreening.clear(
345 fDataPerElement[i]->fCorFirstMoment.clea
346 fDataPerElement[i]->fCorSecondMoment.cle
347 delete fDataPerElement[i];
348 }
349 }
350 fDataPerElement.clear();
351 }
352
353
354 void G4GSPWACorrections::ClearDataPerMaterial(
355 for (std::size_t i=0; i<fDataPerMaterial.siz
356 if (fDataPerMaterial[i]) {
357 fDataPerMaterial[i]->fCorScreening.clear
358 fDataPerMaterial[i]->fCorFirstMoment.cle
359 fDataPerMaterial[i]->fCorSecondMoment.cl
360 delete fDataPerMaterial[i];
361 }
362 }
363 fDataPerMaterial.clear();
364 }
365