| 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: G4GSMottCorrection
31 //
32 // Author: Mihaly Novak
33 //
34 // Creation date: 23.08.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 "G4GSMottCorrection.hh"
44
45 #include "G4PhysicalConstants.hh"
46 #include "zlib.h"
47 #include "Randomize.hh"
48 #include "G4Log.hh"
49 #include "G4Exp.hh"
50
51 #include "G4ProductionCutsTable.hh"
52 #include "G4MaterialCutsCouple.hh"
53 #include "G4Material.hh"
54 #include "G4ElementVector.hh"
55 #include "G4Element.hh"
56 #include "G4EmParameters.hh"
57
58 #include <iostream>
59 #include <fstream>
60 #include <cmath>
61 #include <algorithm>
62
63
64 const std::string G4GSMottCorrection::gElemSym
65 "C" ,"N" ,"O" ,"F" ,"Ne","Na","Mg","Al","Si",
66 "Ti","V" ,"Cr","Mn","Fe","Co","Ni","Cu","Zn",
67 "Sr","Y" ,"Zr","Nb","Mo","Tc","Ru","Rh","Pd",
68 "Xe","Cs","Ba","La","Ce","Pr","Nd","Pm","Sm",
69 "Yb","Lu","Hf","Ta","W" ,"Re","Os","Ir","Pt",
70 "Rn","Fr","Ra","Ac","Th","Pa","U" ,"Np","Pu",
71
72 G4GSMottCorrection::G4GSMottCorrection(G4bool
73 // init grids related data member values
74 fMaxEkin = CLHEP::electron_mass_c2*(1
75 fLogMinEkin = G4Log(gMinEkin);
76 fInvLogDelEkin = (gNumEkin-gNumBeta2)/G4Log
77 G4double pt2 = gMidEkin*(gMidEkin+2.0*CLH
78 fMinBeta2 = pt2/(pt2+CLHEP::electron_m
79 fInvDelBeta2 = (gNumBeta2-1.)/(gMaxBeta2-
80 fInvDelDelta = (gNumDelta-1.)/gMaxDelta;
81 fInvDelAngle = gNumAngle-1.;
82 }
83
84
85 G4GSMottCorrection::~G4GSMottCorrection() {
86 ClearMCDataPerElement();
87 ClearMCDataPerMaterial();
88 }
89
90
91 void G4GSMottCorrection::GetMottCorrectionFact
92
93 G4int ekinIndxLow = 0;
94 G4double remRfaction = 0.;
95 if (beta2>=gMaxBeta2) {
96 ekinIndxLow = gNumEkin - 1;
97 // remRfaction = -1.
98 } else if (beta2>=fMinBeta2) { // linear in
99 remRfaction = (beta2 - fMinBeta2) * fInv
100 ekinIndxLow = (G4int)remRfaction;
101 remRfaction -= ekinIndxLow;
102 ekinIndxLow += (gNumEkin - gNumBeta2);
103 } else if (logekin>=fLogMinEkin) {
104 remRfaction = (logekin - fLogMinEkin) *
105 ekinIndxLow = (G4int)remRfaction;
106 remRfaction -= ekinIndxLow;
107 } // the defaults otherwise i.e. use the low
108 //
109 DataPerEkin *perEkinLow = fMCDataPerMateria
110 mcToScr = perEkinLow->fMCScreening;
111 mcToQ1 = perEkinLow->fMCFirstMoment;
112 mcToG2PerG1 = perEkinLow->fMCSecondMoment;
113 if (remRfaction>0.) {
114 DataPerEkin *perEkinHigh = fMCDataPerMater
115 mcToScr += remRfaction*(perEkinHigh->
116 mcToQ1 += remRfaction*(perEkinHigh->
117 mcToG2PerG1 += remRfaction*(perEkinHigh->
118 }
119 }
120
121
122 // accept cost if rndm [0,1] < return value
123 double G4GSMottCorrection::GetMottRejectionVal
124
125 G4double val = 1.0;
126 G4double delta = q1/(0.5+q1);
127 // check if converged to 1 for all angles =>
128 if (delta>=gMaxDelta) {
129 return val;
130 }
131 //
132 // check if kinetic energy index needs to be
133 if (ekindx<0) {
134 G4int ekinIndxLow = 0;
135 G4double probIndxHigh = 0.; // will be th
136 if (beta2>gMaxBeta2) {
137 ekinIndxLow = gNumEkin - 1;
138 // probIndxHigh = -1.
139 } else if (beta2>=fMinBeta2) { // linea
140 probIndxHigh = (beta2 - fMinBeta2) * fI
141 ekinIndxLow = (G4int)probIndxHigh;
142 probIndxHigh -= ekinIndxLow;
143 ekinIndxLow += (gNumEkin - gNumBeta2);
144 } else if (logekin>fLogMinEkin) { // linea
145 probIndxHigh = (logekin - fLogMinEkin)
146 ekinIndxLow = (G4int)probIndxHigh;
147 probIndxHigh -= ekinIndxLow;
148 } // the defaults otherwise i.e. use the l
149 //
150 // check if need to take the higher ekin i
151 if (G4UniformRand()<probIndxHigh) {
152 ++ekinIndxLow;
153 }
154 // set kinetic energy grid index
155 ekindx = ekinIndxLow;
156 }
157 // check if delta value index needs to be de
158 // by the caller but the ekindx will be -1:
159 if (deltindx<0) {
160 // note: delta is for sure < gMaxDelta at
161 G4double probIndxHigh = delta*fInvDelDelta
162 G4int deltIndxLow = (G4int)probIndxHig
163 probIndxHigh -= deltIndxLow;
164 // check if need to take the higher delta
165 if (G4UniformRand()<probIndxHigh) {
166 ++deltIndxLow;
167 }
168 // set the delta value grid index
169 deltindx = deltIndxLow;
170 }
171 //
172 // get the corresponding distribution
173 DataPerDelta *perDelta = fMCDataPerMaterial
174 //
175 // determine lower index of the angular bin
176 G4double ang = std::sqrt(0.5*(1.-cos
177 G4double remRfaction = ang*fInvDelAngle;
178 G4int angIndx = (G4int)remRfaction;
179 remRfaction -= angIndx;
180 if (angIndx<gNumAngle-2) { // normal case: l
181 val = remRfaction*(perDelta->fRej
182 } else { // last bin
183 G4double dum = ang-1.+1./fInvDelAngle;
184 val = perDelta->fSA + dum*(perDel
185 }
186 return val;
187 }
188
189
190 void G4GSMottCorrection::Initialise() {
191 // load Mott-correction data for each elemen
192 InitMCDataPerElement();
193 // clrea Mott-correction data per material
194 ClearMCDataPerMaterial();
195 // initialise Mott-correction data for the m
196 InitMCDataPerMaterials();
197 }
198
199
200 void G4GSMottCorrection::InitMCDataPerElement(
201 // do it only once
202 if (fMCDataPerElement.size()<gMaxZet+1) {
203 fMCDataPerElement.resize(gMaxZet+1,nullptr
204 }
205 // loop over all materials, for those that a
206 // corresponding data has not been loaded ye
207 G4ProductionCutsTable *thePCTable = G4Produc
208 G4int numMatCuts = (G4int)thePCTable->GetTab
209 for (G4int imc=0; imc<numMatCuts; ++imc) {
210 const G4MaterialCutsCouple *matCut = thePC
211 if (!matCut->IsUsed()) {
212 continue;
213 }
214 const G4Material *mat = matCut->
215 const G4ElementVector *elemVect = mat->Get
216 //
217 std::size_t numElems = elemVect->size();
218 for (std::size_t ielem=0; ielem<numElems;
219 const G4Element *elem = (*elemVect)[iele
220 G4int izet = G4lrint(elem->GetZ());
221 if (izet>gMaxZet) {
222 izet = gMaxZet;
223 }
224 if (!fMCDataPerElement[izet]) {
225 LoadMCDataElement(elem);
226 }
227 }
228 }
229 }
230
231
232 void G4GSMottCorrection::InitMCDataPerMaterial
233 // prepare size of the container
234 std::size_t numMaterials = G4Material::GetNu
235 if (fMCDataPerMaterial.size()!=numMaterials)
236 fMCDataPerMaterial.resize(numMaterials);
237 }
238 // init. Mott-correction data for the Materi
239 G4ProductionCutsTable *thePCTable = G4Produc
240 G4int numMatCuts = (G4int)thePCTable->GetTab
241 for (G4int imc=0; imc<numMatCuts; ++imc) {
242 const G4MaterialCutsCouple *matCut = thePC
243 if (!matCut->IsUsed()) {
244 continue;
245 }
246 const G4Material *mat = matCut->GetMateria
247 if (!fMCDataPerMaterial[mat->GetIndex()])
248 InitMCDataMaterial(mat);
249 }
250 }
251 }
252
253
254 // it's called only if data has not been loade
255 void G4GSMottCorrection::LoadMCDataElement(con
256 // allocate memory
257 G4int izet = elem->GetZasInt();
258 if (izet>gMaxZet) {
259 izet = gMaxZet;
260 }
261 auto perElem = new DataPerMaterial();
262 AllocateDataPerMaterial(perElem);
263 fMCDataPerElement[izet] = perElem;
264 //
265 // load data from file
266 std::string path = G4EmParameters::Instance(
267 if (fIsElectron) {
268 path += "/msc_GS/MottCor/el/";
269 } else {
270 path += "/msc_GS/MottCor/pos/";
271 }
272 std::string fname = path+"rej_"+gElemSymbols
273 std::istringstream infile(std::ios::in);
274 ReadCompressedFile(fname, infile);
275 // check if file is open !!!
276 for (G4int iek=0; iek<gNumEkin; ++iek) {
277 DataPerEkin *perEkin = perElem->fDataPerEk
278 // 1. get the 3 Mott-correction factors fo
279 infile >> perEkin->fMCScreening;
280 infile >> perEkin->fMCFirstMoment;
281 infile >> perEkin->fMCSecondMoment;
282 // 2. load each data per delta:
283 for (G4int idel=0; idel<gNumDelta; ++idel)
284 DataPerDelta *perDelta = perEkin->fDataP
285 // 2./a. : first the rejection function
286 for (G4int iang=0; iang<gNumAngle; ++ian
287 infile >> perDelta->fRejFuntion[iang];
288 }
289 // 2./b. : then the 4 spline parameter f
290 infile >> perDelta->fSA;
291 infile >> perDelta->fSB;
292 infile >> perDelta->fSC;
293 infile >> perDelta->fSD;
294 }
295 }
296 }
297
298 // uncompress one data file into the input str
299 void G4GSMottCorrection::ReadCompressedFile(co
300 std::string *dataString = nullptr;
301 std::string compfilename(fname+".z");
302 // create input stream with binary mode oper
303 std::ifstream in(compfilename, std::ios::bin
304 if (in.good()) {
305 // get current position in the stream (wa
306 std::streamoff fileSize = in.tellg();
307 // set current position being the beginni
308 in.seekg(0,std::ios::beg);
309 // create (zlib) byte buffer for the data
310 Bytef *compdata = new Bytef[fileSize];
311 while(in) {
312 in.read((char*)compdata, fileSize);
313 }
314 // create (zlib) byte buffer for the unco
315 uLongf complen = (uLongf)(fileSize*4);
316 Bytef *uncompdata = new Bytef[complen];
317 while (Z_OK!=uncompress(uncompdata, &comp
318 // increase uncompressed byte buffer
319 delete[] uncompdata;
320 complen *= 2;
321 uncompdata = new Bytef[complen];
322 }
323 // delete the compressed data buffer
324 delete [] compdata;
325 // create a string from the uncompressed
326 dataString = new std::string((char*)uncom
327 // delete the uncompressed data buffer
328 delete [] uncompdata;
329 } else {
330 std::string msg = " Problem while trying
331 G4Exception("G4GSMottCorrection::ReadCompr
332 return;
333 }
334 // create the input string stream from the d
335 if (dataString) {
336 iss.str(*dataString);
337 in.close();
338 delete dataString;
339 }
340 }
341
342
343 void G4GSMottCorrection::InitMCDataMaterial(co
344 constexpr G4double const1 = 7821.6; /
345 constexpr G4double const2 = 0.1569; /
346 constexpr G4double finstrc2 = 5.325135453E-5
347
348 G4double constFactor = CLHEP::electro
349 constFactor *= constFactor;
350 // allocate memory
351 auto perMat = new DataPerMaterial();
352 AllocateDataPerMaterial(perMat);
353 fMCDataPerMaterial[mat->GetIndex()] = perMat
354 //
355 const G4ElementVector* elemVect =
356 const G4int numElems =
357 const G4double* nbAtomsPerVolVect =
358 G4double totNbAtomsPerVol =
359 //
360 // 1. Compute material dependent part of Mol
361 // (with \xi=1 (i.e. total sub-threshold
362 G4double moliereBc = 0.0;
363 G4double moliereXc2 = 0.0;
364 G4double zs = 0.0;
365 G4double ze = 0.0;
366 G4double zx = 0.0;
367 G4double sa = 0.0;
368 G4double xi = 1.0;
369 for (G4int ielem=0; ielem<numElems; ++ielem)
370 G4double zet = (*elemVect)[ielem]->GetZ();
371 if (zet>gMaxZet) {
372 zet = (G4double)gMaxZet;
373 }
374 G4double iwa = (*elemVect)[ielem]->GetN()
375 G4double ipz = nbAtomsPerVolVect[ielem]/t
376 G4double dum = ipz*zet*(zet+xi);
377 zs += dum;
378 ze += dum*(-2.0/3.0)*G4Log(zet);
379 zx += dum*G4Log(1.0+3.34*finstrc
380 sa += ipz*iwa;
381 }
382 G4double density = mat->GetDensity()*CLHEP::
383 //
384 G4double z0 = (0.0 == sa) ? 0.0 : zs/sa;
385 G4double z1 = (0.0 == zs) ? 0.0 : (ze - zx)/
386
387 moliereBc = const1*density*z0*G4Exp(z1); /
388 moliereXc2 = const2*density*z0; // [MeV2/cm
389 // change to Geant4 internal units of 1/leng
390 moliereBc *= 1.0/CLHEP::cm;
391 moliereXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::c
392 //
393 // 2. loop over the kinetic energy grid
394 for (G4int iek=0; iek<gNumEkin; ++iek) {
395 // 2./a. set current kinetic energy and pt
396 G4double ekin = G4Exp(fLogMinEk
397 G4double pt2 = ekin*(ekin+2.0*CLHEP::el
398 if (ekin>gMidEkin) {
399 G4double b2 = fMinBeta2+(iek-(gNumEk
400 ekin = CLHEP::electron_mass_c2*(1./std
401 pt2 = ekin*(ekin+2.0*CLHEP::electron_
402 }
403 // 2./b. loop over the elements at the cur
404 for (G4int ielem=0; ielem<numElems; ++iele
405 const G4Element *elem = (*elemVect)[iele
406 G4double zet = elem->GetZ();
407 if (zet>gMaxZet) {
408 zet = (G4double)gMaxZet;
409 }
410 G4int izet = G4lrint(zet);
411 // xi should be one i.e. z(z+1) since to
412 G4double nZZPlus1 = nbAtomsPerVolVect[i
413 G4double Z23 = std::pow(zet,2./3.)
414 //
415 DataPerEkin *perElemPerEkin = fMCDataPe
416 DataPerEkin *perMatPerEkin = perMat->f
417 //
418 // 2./b./(i) Add the 3 Mott-correction f
419 G4double mcScrCF = perElemPerEkin->fMCSc
420 // compute the screening parameter corre
421 // src_{mc} = C \exp\left[ \frac{ \sum_i
422 // with C = \frac{(mc^2)^\alpha^2} {4(pc
423 // here we compute the \sum_i n_i Z_i(Z_
424 perMatPerEkin->fMCScreening += nZZPlus1*
425 // compute the corrected screening param
426 // src(Z_i)_{mc} = \frac{(mc^2)^\alpha^2
427 mcScrCF *= constFactor*Z23/(4.*pt2);
428 // compute first moment correction facto
429 // q1_{mc} = \frac{ \sum_i n_i Z_i(Z_i+1
430 // where:
431 // A_i(src(Z_i)_{mc}) = [\ln(1+1/src(Z_i
432 // B_i = \beta_i \gamma_i with beta_i(Z_
433 // and \gamma_i = \sigma(Z_i)_{el}^(MC-D
434 // C(src_{mc}) = [\ln(1+1/src_{mc}) - 1/
435 // A_i x B_i is stored in file per e-/e+
436 // here we compute the \sum_i n_i Z_i(Z_
437 perMatPerEkin->fMCFirstMoment += nZZPlus
438 // compute the second moment correction
439 // [G2/G1]_{mc} = \frac{ \sum_i n_i Z_i(
440 // with A_i(Z_i) = G2(Z_i)^{PWA}/G1(Z_i)
441 // here we compute the \sum_i n_i Z_i(Z_
442 perMatPerEkin->fMCSecondMoment += nZZPlu
443 //
444 // 2./b./(ii) Go for the rejection funti
445 // I. loop over delta values
446 for (G4int idel=0; idel<gNumDelta; ++ide
447 DataPerDelta *perMatPerDelta = perMat
448 DataPerDelta *perElemPerDelta = perEle
449 // I./a. loop over angles (i.e. the \s
450 for (G4int iang=0; iang<gNumAngle; ++i
451 perMatPerDelta->fRejFuntion[iang] +=
452 }
453 // I./b. get the last bin spline param
454 perMatPerDelta->fSA += nZZPlus1*perEle
455 perMatPerDelta->fSB += nZZPlus1*perEle
456 perMatPerDelta->fSC += nZZPlus1*perEle
457 perMatPerDelta->fSD += nZZPlus1*perEle
458 }
459 //
460 // 2./b./(iii) When the last element has
461 if (ielem==numElems-1) {
462 //
463 // 1. the remaining part of the sreeni
464 // (Moliere screening parameter = m
465 G4double dumScr = G4Exp(perMatPerEki
466 perMatPerEkin->fMCScreening = constFac
467 //
468 // 2. the remaining part of the first
469 // screening parameter (= (mc^2)^\a
470 G4double scrCorTed = constFactor*dumSc
471 G4double dum0 = G4Log(1.+1./scrCo
472 perMatPerEkin->fMCFirstMoment = perMat
473 //
474 // 3. the remaining part of the second
475 // screening parameter
476 G4double G2PerG1 = 3.*(1.+scrCorTed
477 perMatPerEkin->fMCSecondMoment = perMa
478 //
479 // 4. scale the maximum of the rejecti
480 // I. loop over delta values
481 for (G4int idel=0; idel<gNumDelta; ++i
482 DataPerDelta *perMatPerDelta = perM
483 G4double maxVal = -1.;
484 // II. llop over angles
485 for (G4int iang=0; iang<gNumAngle; +
486 if (perMatPerDelta->fRejFuntion[ia
487 maxVal = perMatPerDelta->fRejFun
488 }
489 for (G4int iang=0; iang<gNumAngle; +
490 perMatPerDelta->fRejFuntion[iang]
491 }
492 perMatPerDelta->fSA /= maxVal;
493 perMatPerDelta->fSB /= maxVal;
494 perMatPerDelta->fSC /= maxVal;
495 perMatPerDelta->fSD /= maxVal;
496 }
497 }
498 }
499 }
500 }
501
502
503 void G4GSMottCorrection::AllocateDataPerMateri
504 data->fDataPerEkin = new DataPerEkin*[gNumEk
505 for (G4int iek=0; iek<gNumEkin; ++iek) {
506 auto perEkin = new DataPerEkin();
507 perEkin->fDataPerDelta = new DataPerDelta*
508 for (G4int idel=0; idel<gNumDelta; ++idel)
509 auto perDelta = new DataP
510 perDelta->fRejFuntion = new doubl
511 perEkin->fDataPerDelta[idel] = perDelta;
512 }
513 data->fDataPerEkin[iek] = perEkin;
514 }
515 }
516
517 void G4GSMottCorrection::DeAllocateDataPerMate
518 for (G4int iek=0; iek<gNumEkin; ++iek) {
519 DataPerEkin *perEkin = data->fDataPerEkin[
520 for (G4int idel=0; idel<gNumDelta; ++idel)
521 DataPerDelta *perDelta = perEkin->fDataP
522 delete [] perDelta->fRejFuntion;
523 delete perDelta;
524 }
525 delete [] perEkin->fDataPerDelta;
526 delete perEkin;
527 }
528 delete [] data->fDataPerEkin;
529 }
530
531
532 void G4GSMottCorrection::ClearMCDataPerElement
533 for (std::size_t i=0; i<fMCDataPerElement.si
534 if (fMCDataPerElement[i]) {
535 DeAllocateDataPerMaterial(fMCDataPerElem
536 delete fMCDataPerElement[i];
537 }
538 }
539 fMCDataPerElement.clear();
540 }
541
542 void G4GSMottCorrection::ClearMCDataPerMateria
543 for (std::size_t i=0; i<fMCDataPerMaterial.s
544 if (fMCDataPerMaterial[i]) {
545 DeAllocateDataPerMaterial(fMCDataPerMate
546 delete fMCDataPerMaterial[i];
547 }
548 }
549 fMCDataPerMaterial.clear();
550 }
551