Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4GSMottCorrection.cc
이 파일의 문서화 페이지로 가기
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 // $Id: $
27 //
28 // ----------------------------------------------------------------------------
29 //
30 //
31 // File name: G4GSMottCorrection
32 //
33 // Author: Mihaly Novak
34 //
35 // Creation date: 23.08.2017
36 //
37 // Modifications:
38 // 02.02.2018 M.Novak: fixed initialization of first moment correction.
39 //
40 // Class description: see the header file.
41 //
42 // -----------------------------------------------------------------------------
43 
44 #include "G4GSMottCorrection.hh"
45 
46 #include "G4PhysicalConstants.hh"
47 #include "zlib.h"
48 #include "Randomize.hh"
49 #include "G4Log.hh"
50 #include "G4Exp.hh"
51 
52 #include "G4ProductionCutsTable.hh"
53 #include "G4MaterialCutsCouple.hh"
54 #include "G4Material.hh"
55 #include "G4ElementVector.hh"
56 #include "G4Element.hh"
57 
58 #include <iostream>
59 #include <fstream>
60 #include <cmath>
61 #include <algorithm>
62 
63 
64 const std::string G4GSMottCorrection::gElemSymbols[] = {"H","He","Li","Be","B" ,
65  "C" ,"N" ,"O" ,"F" ,"Ne","Na","Mg","Al","Si","P" , "S","Cl","Ar","K" ,"Ca","Sc",
66  "Ti","V" ,"Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr","Rb",
67  "Sr","Y" ,"Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I" ,
68  "Xe","Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm",
69  "Yb","Lu","Hf","Ta","W" ,"Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At",
70  "Rn","Fr","Ra","Ac","Th","Pa","U" ,"Np","Pu","Am","Cm","Bk","Cf"};
71 
72 G4GSMottCorrection::G4GSMottCorrection(G4bool iselectron) : fIsElectron(iselectron) {
73  // init grids related data member values
74  fMaxEkin = CLHEP::electron_mass_c2*(1./std::sqrt(1.-gMaxBeta2)-1.);
82 }
83 
84 
88 }
89 
90 
92  G4double &mcToQ1, G4double &mcToG2PerG1) {
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 interpolation on \beta^2
99  remRfaction = (beta2 - fMinBeta2) * fInvDelBeta2;
100  ekinIndxLow = (G4int)remRfaction;
101  remRfaction -= ekinIndxLow;
102  ekinIndxLow += (gNumEkin - gNumBeta2);
103  } else if (logekin>=fLogMinEkin) {
104  remRfaction = (logekin - fLogMinEkin) * fInvLogDelEkin;
105  ekinIndxLow = (G4int)remRfaction;
106  remRfaction -= ekinIndxLow;
107  } // the defaults otherwise i.e. use the lowest energy values when ekin is smaller than the minum ekin
108  //
109  DataPerEkin *perEkinLow = fMCDataPerMaterial[matindx]->fDataPerEkin[ekinIndxLow];
110  mcToScr = perEkinLow->fMCScreening;
111  mcToQ1 = perEkinLow->fMCFirstMoment;
112  mcToG2PerG1 = perEkinLow->fMCSecondMoment;
113  if (remRfaction>0.) {
114  DataPerEkin *perEkinHigh = fMCDataPerMaterial[matindx]->fDataPerEkin[ekinIndxLow+1];
115  mcToScr += remRfaction*(perEkinHigh->fMCScreening - perEkinLow->fMCScreening);
116  mcToQ1 += remRfaction*(perEkinHigh->fMCFirstMoment - perEkinLow->fMCFirstMoment);
117  mcToG2PerG1 += remRfaction*(perEkinHigh->fMCSecondMoment - perEkinLow->fMCSecondMoment);
118  }
119 }
120 
121 
122 // accept cost if rndm [0,1] < return value
124  G4int matindx, G4int &ekindx, G4int &deltindx) {
125  G4double val = 1.0;
126  G4double delta = q1/(0.5+q1);
127  // check if converged to 1 for all angles => accept cost
128  if (delta>=gMaxDelta) {
129  return val;
130  }
131  //
132  // check if kinetic energy index needs to be determined
133  if (ekindx<0) {
134  G4int ekinIndxLow = 0;
135  G4double probIndxHigh = 0.; // will be the prob. of taking the ekinIndxLow+1 bin
136  if (beta2>gMaxBeta2) {
137  ekinIndxLow = gNumEkin - 1;
138  // probIndxHigh = -1.
139  } else if (beta2>=fMinBeta2) { // linear interpolation on \beta^2
140  probIndxHigh = (beta2 - fMinBeta2) * fInvDelBeta2;
141  ekinIndxLow = (G4int)probIndxHigh;
142  probIndxHigh -= ekinIndxLow;
143  ekinIndxLow += (gNumEkin - gNumBeta2);
144  } else if (logekin>fLogMinEkin) { // linear interpolation on \ln(E_{kin})
145  probIndxHigh = (logekin - fLogMinEkin) * fInvLogDelEkin;
146  ekinIndxLow = (G4int)probIndxHigh;
147  probIndxHigh -= ekinIndxLow;
148  } // the defaults otherwise i.e. use the lowest energy values when ekin is smaller than the minum ekin
149  //
150  // check if need to take the higher ekin index
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 determined (note: in case of single scattering deltindx will be set to 0 by
158  // by the caller but the ekindx will be -1: kinetic energy index is not known but the delta index is known)
159  if (deltindx<0) {
160  // note: delta is for sure < gMaxDelta at this point ( and minimum delta value is 0)
161  G4double probIndxHigh = delta*fInvDelDelta; // will be the prob. of taking the deltIndxLow+1 bin
162  G4int deltIndxLow = (G4int)probIndxHigh;
163  probIndxHigh -= deltIndxLow;
164  // check if need to take the higher delta index
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[matindx]->fDataPerEkin[ekindx]->fDataPerDelta[deltindx];
174  //
175  // determine lower index of the angular bin
176  G4double ang = std::sqrt(0.5*(1.-cost)); // sin(0.5\theta) in [0,1]
177  G4double remRfaction = ang*fInvDelAngle;
178  G4int angIndx = (G4int)remRfaction;
179  remRfaction -= angIndx;
180  if (angIndx<gNumAngle-2) { // normal case: linear interpolation
181  val = remRfaction*(perDelta->fRejFuntion[angIndx+1]-perDelta->fRejFuntion[angIndx]) + perDelta->fRejFuntion[angIndx];
182  } else { // last bin
183  G4double dum = ang-1.+1./fInvDelAngle;
184  val = perDelta->fSA + dum*(perDelta->fSB + dum*(perDelta->fSC + dum*perDelta->fSD));
185  }
186  return val;
187 }
188 
189 
191  // load Mott-correction data for each elements that belongs to materials that are used in the detector
193  // clrea Mott-correction data per material
195  // initialise Mott-correction data for the materials that are used in the detector
197 }
198 
199 
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 are used check the list of elements and load data from file if the
206  // corresponding data has not been loaded yet
208  size_t numMatCuts = thePCTable->GetTableSize();
209  for (size_t imc=0; imc<numMatCuts; ++imc) {
210  const G4MaterialCutsCouple *matCut = thePCTable->GetMaterialCutsCouple(imc);
211  if (!matCut->IsUsed()) {
212  continue;
213  }
214  const G4Material *mat = matCut->GetMaterial();
215  const G4ElementVector *elemVect = mat->GetElementVector();
216  //
217  size_t numElems = elemVect->size();
218  for (size_t ielem=0; ielem<numElems; ++ielem) {
219  const G4Element *elem = (*elemVect)[ielem];
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 
233  // prepare size of the container
234  size_t numMaterials = G4Material::GetNumberOfMaterials();
235  if (fMCDataPerMaterial.size()!=numMaterials) {
236  fMCDataPerMaterial.resize(numMaterials);
237  }
238  // init. Mott-correction data for the Materials that are used in the geometry
240  size_t numMatCuts = thePCTable->GetTableSize();
241  for (size_t imc=0; imc<numMatCuts; ++imc) {
242  const G4MaterialCutsCouple *matCut = thePCTable->GetMaterialCutsCouple(imc);
243  if (!matCut->IsUsed()) {
244  continue;
245  }
246  const G4Material *mat = matCut->GetMaterial();
247  if (!fMCDataPerMaterial[mat->GetIndex()]) {
248  InitMCDataMaterial(mat);
249  }
250  }
251 }
252 
253 
254 // it's called only if data has not been loaded for this element yet
256  // allocate memory
257  G4int izet = G4lrint(elem->GetZ());
258  if (izet>gMaxZet) {
259  izet = gMaxZet;
260  }
261  DataPerMaterial *perElem = new DataPerMaterial();
262  AllocateDataPerMaterial(perElem);
263  fMCDataPerElement[izet] = perElem;
264  //
265  // load data from file
266  char* tmppath = getenv("G4LEDATA");
267  if (!tmppath) {
268  G4Exception("G4GSMottCorrection::LoadMCDataElement()","em0006",
270  "Environment variable G4LEDATA not defined");
271  return;
272  }
273  std::string path(tmppath);
274  if (fIsElectron) {
275  path += "/msc_GS/MottCor/el/";
276  } else {
277  path += "/msc_GS/MottCor/pos/";
278  }
279  std::string fname = path+"rej_"+gElemSymbols[izet-1];
280  std::istringstream infile(std::ios::in);
281  ReadCompressedFile(fname, infile);
282  // check if file is open !!!
283  for (G4int iek=0; iek<gNumEkin; ++iek) {
284  DataPerEkin *perEkin = perElem->fDataPerEkin[iek];
285  // 1. get the 3 Mott-correction factors for the current kinetic energy
286  infile >> perEkin->fMCScreening;
287  infile >> perEkin->fMCFirstMoment;
288  infile >> perEkin->fMCSecondMoment;
289  // 2. load each data per delta:
290  for (G4int idel=0; idel<gNumDelta; ++idel) {
291  DataPerDelta *perDelta = perEkin->fDataPerDelta[idel];
292  // 2./a. : first the rejection function values
293  for (G4int iang=0; iang<gNumAngle; ++iang) {
294  infile >> perDelta->fRejFuntion[iang];
295  }
296  // 2./b. : then the 4 spline parameter for the last bin
297  infile >> perDelta->fSA;
298  infile >> perDelta->fSB;
299  infile >> perDelta->fSC;
300  infile >> perDelta->fSD;
301  }
302  }
303 }
304 
305 // uncompress one data file into the input string stream
306 void G4GSMottCorrection::ReadCompressedFile(std::string fname, std::istringstream &iss) {
307  std::string *dataString = nullptr;
308  std::string compfilename(fname+".z");
309  // create input stream with binary mode operation and positioning at the end of the file
310  std::ifstream in(compfilename, std::ios::binary | std::ios::ate);
311  if (in.good()) {
312  // get current position in the stream (was set to the end)
313  int fileSize = in.tellg();
314  // set current position being the beginning of the stream
315  in.seekg(0,std::ios::beg);
316  // create (zlib) byte buffer for the data
317  Bytef *compdata = new Bytef[fileSize];
318  while(in) {
319  in.read((char*)compdata, fileSize);
320  }
321  // create (zlib) byte buffer for the uncompressed data
322  uLongf complen = (uLongf)(fileSize*4);
323  Bytef *uncompdata = new Bytef[complen];
324  while (Z_OK!=uncompress(uncompdata, &complen, compdata, fileSize)) {
325  // increase uncompressed byte buffer
326  delete[] uncompdata;
327  complen *= 2;
328  uncompdata = new Bytef[complen];
329  }
330  // delete the compressed data buffer
331  delete [] compdata;
332  // create a string from the uncompressed data (will be deallocated by the caller)
333  dataString = new std::string((char*)uncompdata, (long)complen);
334  // delete the uncompressed data buffer
335  delete [] uncompdata;
336  } else {
337  std::string msg = " Problem while trying to read " + compfilename + " data file.\n";
338  G4Exception("G4GSMottCorrection::ReadCompressedFile","em0006", FatalException,msg.c_str());
339  return;
340  }
341  // create the input string stream from the data string
342  if (dataString) {
343  iss.str(*dataString);
344  in.close();
345  delete dataString;
346  }
347 }
348 
349 
351  constexpr G4double const1 = 7821.6; // [cm2/g]
352  constexpr G4double const2 = 0.1569; // [cm2 MeV2 / g]
353  constexpr G4double finstrc2 = 5.325135453E-5; // fine-structure const. square
354 
356  constFactor *= constFactor; // (mc^2)^2\alpha^2/( C_{TF}^2)
357  // allocate memory
358  DataPerMaterial *perMat = new DataPerMaterial();
359  AllocateDataPerMaterial(perMat);
360  fMCDataPerMaterial[mat->GetIndex()] = perMat;
361  //
362  const G4ElementVector* elemVect = mat->GetElementVector();
363  const G4int numElems = mat->GetNumberOfElements();
364  const G4double* nbAtomsPerVolVect = mat->GetVecNbOfAtomsPerVolume();
365  G4double totNbAtomsPerVol = mat->GetTotNbOfAtomsPerVolume();
366  //
367  // 1. Compute material dependent part of Moliere's b_c \chi_c^2
368  // (with \xi=1 (i.e. total sub-threshold scattering power correction)
369  G4double moliereBc = 0.0;
370  G4double moliereXc2 = 0.0;
371  G4double zs = 0.0;
372  G4double ze = 0.0;
373  G4double zx = 0.0;
374  G4double sa = 0.0;
375  G4double xi = 1.0;
376  for (G4int ielem=0; ielem<numElems; ++ielem) {
377  G4double zet = (*elemVect)[ielem]->GetZ();
378  if (zet>gMaxZet) {
379  zet = (G4double)gMaxZet;
380  }
381  G4double iwa = (*elemVect)[ielem]->GetN();
382  G4double ipz = nbAtomsPerVolVect[ielem]/totNbAtomsPerVol;
383  G4double dum = ipz*zet*(zet+xi);
384  zs += dum;
385  ze += dum*(-2.0/3.0)*G4Log(zet);
386  zx += dum*G4Log(1.0+3.34*finstrc2*zet*zet);
387  sa += ipz*iwa;
388  }
389  G4double density = mat->GetDensity()*CLHEP::cm3/CLHEP::g; // [g/cm3]
390  //
391  moliereBc = const1*density*zs/sa*G4Exp(ze/zs)/G4Exp(zx/zs); //[1/cm]
392  moliereXc2 = const2*density*zs/sa; // [MeV2/cm]
393  // change to Geant4 internal units of 1/length and energ2/length
394  moliereBc *= 1.0/CLHEP::cm;
395  moliereXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm;
396  //
397  // 2. loop over the kinetic energy grid
398  for (G4int iek=0; iek<gNumEkin; ++iek) {
399  // 2./a. set current kinetic energy and pt2 value
401  G4double pt2 = ekin*(ekin+2.0*CLHEP::electron_mass_c2);
402  if (ekin>gMidEkin) {
403  G4double b2 = fMinBeta2+(iek-(gNumEkin-gNumBeta2))/fInvDelBeta2;
404  ekin = CLHEP::electron_mass_c2*(1./std::sqrt(1.-b2)-1.);
405  pt2 = ekin*(ekin+2.0*CLHEP::electron_mass_c2);
406  }
407  // 2./b. loop over the elements at the current kinetic energy point
408  for (G4int ielem=0; ielem<numElems; ++ielem) {
409  const G4Element *elem = (*elemVect)[ielem];
410  G4double zet = elem->GetZ();
411  if (zet>gMaxZet) {
412  zet = (G4double)gMaxZet;
413  }
414  G4int izet = G4lrint(zet);
415  // xi should be one i.e. z(z+1) since total sub-threshold scattering power correction
416  G4double nZZPlus1 = nbAtomsPerVolVect[ielem]*zet*(zet+1.0)/totNbAtomsPerVol;
417  G4double Z23 = std::pow(zet,2./3.);
418  //
419  DataPerEkin *perElemPerEkin = fMCDataPerElement[izet]->fDataPerEkin[iek];
420  DataPerEkin *perMatPerEkin = perMat->fDataPerEkin[iek];
421  //
422  // 2./b./(i) Add the 3 Mott-correction factors
423  G4double mcScrCF = perElemPerEkin->fMCScreening; // \kappa_i[1.13+3.76(\alpha Z_i)^2] with \kappa_i=scr_mc/scr_sr
424  // compute the screening parameter correction factor (Z_i contribution to the material)
425  // src_{mc} = C \exp\left[ \frac{ \sum_i n_i Z_i(Z_i+1)\ln[Z_{i}^{2/3}\kappa_i(1.13+3.76(\alpha Z_i)^2)] } {\sum_i n_i Z_i(Z_i+1)}
426  // with C = \frac{(mc^2)^\alpha^2} {4(pc)^2 C_{TF}^2} = constFactor/(4*(pc)^2)
427  // here we compute the \sum_i n_i Z_i(Z_i+1)\ln[Z_{i}^{2/3}\kappa_i(1.13+3.76(\alpha Z_i)^2)] part
428  perMatPerEkin->fMCScreening += nZZPlus1*G4Log(Z23*mcScrCF);
429  // compute the corrected screening parameter for the current Z_i and E_{kin}
430  // src(Z_i)_{mc} = \frac{(mc^2)^\alpha^2 Z_i^{2/3}} {4(pc)^2 C_{TF}^2} \kappa_i[1.13+3.76(\alpha Z_i)^2]
431  mcScrCF *= constFactor*Z23/(4.*pt2);
432  // compute first moment correction factor
433  // q1_{mc} = \frac{ \sum_i n_i Z_i(Z_i+1) A_i B_i } {\sum_i n_i Z_i(Z_i+1)} \frac{1}{C}
434  // where:
435  // A_i(src(Z_i)_{mc}) = [\ln(1+1/src(Z_i)_{mc}) - 1/(1+src(Z_i)_{mc})]; where \sigma(Z_i)_{tr1}^(sr) = A_i(src(Z_i)_{mc}) [2\pi r_0 Z_i mc^2/(pc)\beta]^2
436  // B_i = \beta_i \gamma_i with beta_i(Z_i) = \sigma(Z_i)_{tr1}^(PWA)/\sigma(Z_i,src(Z_i)_{mc})_{tr1}^(sr)
437  // and \gamma_i = \sigma(Z_i)_{el}^(MC-DCS)/\sigma(Z_i,src(Z_i)_{mc})_{el}^(sr)
438  // C(src_{mc}) = [\ln(1+1/src_{mc}) - 1/(1+src_{mc})]; where \sigma_{tr1}^(sr) = C(src_{mc}) [2\pi r_0 Z_i mc^2/(pc)\beta]^2
439  // A_i x B_i is stored in file per e-/e+, E_{kin} and Z_i
440  // here we compute the \sum_i n_i Z_i(Z_i+1) A_i B_i part
441  perMatPerEkin->fMCFirstMoment += nZZPlus1*(G4Log(1.+1./mcScrCF)-1./(1.+mcScrCF))*perElemPerEkin->fMCFirstMoment;
442  // compute the second moment correction factor
443  // [G2/G1]_{mc} = \frac{ \sum_i n_i Z_i(Z_i+1) A_i } {\sum_i n_i Z_i(Z_i+1)} \frac{1}{C}
444  // with A_i(Z_i) = G2(Z_i)^{PWA}/G1(Z_i)^{PWA} and C=G2(Z_i,scr_{mc})^{sr}/G1(Z_i,scr_{mc})^{sr}}
445  // here we compute the \sum_i n_i Z_i(Z_i+1) A_i part
446  perMatPerEkin->fMCSecondMoment += nZZPlus1*perElemPerEkin->fMCSecondMoment;
447  //
448  // 2./b./(ii) Go for the rejection funtion part
449  // I. loop over delta values
450  for (G4int idel=0; idel<gNumDelta; ++idel) {
451  DataPerDelta *perMatPerDelta = perMatPerEkin->fDataPerDelta[idel];
452  DataPerDelta *perElemPerDelta = perElemPerEkin->fDataPerDelta[idel];
453  // I./a. loop over angles (i.e. the \sin(0.5\theta) values) and add the rejection function
454  for (G4int iang=0; iang<gNumAngle; ++iang) {
455  perMatPerDelta->fRejFuntion[iang] += nZZPlus1*perElemPerDelta->fRejFuntion[iang];
456  }
457  // I./b. get the last bin spline parameters and add them (a+bx+cx^2+dx^3)
458  perMatPerDelta->fSA += nZZPlus1*perElemPerDelta->fSA;
459  perMatPerDelta->fSB += nZZPlus1*perElemPerDelta->fSB;
460  perMatPerDelta->fSC += nZZPlus1*perElemPerDelta->fSC;
461  perMatPerDelta->fSD += nZZPlus1*perElemPerDelta->fSD;
462  }
463  //
464  // 2./b./(iii) When the last element has been added:
465  if (ielem==numElems-1) {
466  //
467  // 1. the remaining part of the sreening correction and divide the corrected screening par. with Moliere's one:
468  // (Moliere screening parameter = moliereXc2/(4(pc)^2 moliereBc) )
469  G4double dumScr = G4Exp(perMatPerEkin->fMCScreening/zs);
470  perMatPerEkin->fMCScreening = constFactor*dumScr*moliereBc/moliereXc2;
471  //
472  // 2. the remaining part of the first moment correction and divide by the one computed by using the corrected
473  // screening parameter (= (mc^2)^\alpha^2/(4(pc)^2C_{TF}^2) dumScr
474  G4double scrCorTed = constFactor*dumScr/(4.*pt2);
475  G4double dum0 = G4Log(1.+1./scrCorTed);
476  perMatPerEkin->fMCFirstMoment = perMatPerEkin->fMCFirstMoment/(zs*(dum0-1./(1.+scrCorTed)));
477  //
478  // 3. the remaining part of the second moment correction and divide by the one computed by using the corrected
479  // screening parameter
480  G4double G2PerG1 = 3.*(1.+scrCorTed)*((1.+2.*scrCorTed)*dum0-2.)/((1.+scrCorTed)*dum0-1.);
481  perMatPerEkin->fMCSecondMoment = perMatPerEkin->fMCSecondMoment/(zs*G2PerG1);
482  //
483  // 4. scale the maximum of the rejection function to unity and correct the last bin spline parameters as well
484  // I. loop over delta values
485  for (G4int idel=0; idel<gNumDelta; ++idel) {
486  DataPerDelta *perMatPerDelta = perMatPerEkin->fDataPerDelta[idel];
487  G4double maxVal = -1.;
488  // II. llop over angles
489  for (G4int iang=0; iang<gNumAngle; ++iang) {
490  if (perMatPerDelta->fRejFuntion[iang]>maxVal)
491  maxVal = perMatPerDelta->fRejFuntion[iang];
492  }
493  for (G4int iang=0; iang<gNumAngle; ++iang) {
494  perMatPerDelta->fRejFuntion[iang] /=maxVal;
495  }
496  perMatPerDelta->fSA /= maxVal;
497  perMatPerDelta->fSB /= maxVal;
498  perMatPerDelta->fSC /= maxVal;
499  perMatPerDelta->fSD /= maxVal;
500  }
501  }
502  }
503  }
504 }
505 
506 
508  data->fDataPerEkin = new DataPerEkin*[gNumEkin]();
509  for (G4int iek=0; iek<gNumEkin; ++iek) {
510  DataPerEkin *perEkin = new DataPerEkin();
511  perEkin->fDataPerDelta = new DataPerDelta*[gNumDelta]();
512  for (G4int idel=0; idel<gNumDelta; ++idel) {
513  DataPerDelta *perDelta = new DataPerDelta();
514  perDelta->fRejFuntion = new double[gNumAngle]();
515  perEkin->fDataPerDelta[idel] = perDelta;
516  }
517  data->fDataPerEkin[iek] = perEkin;
518  }
519 }
520 
522  for (G4int iek=0; iek<gNumEkin; ++iek) {
523  DataPerEkin *perEkin = data->fDataPerEkin[iek]; //new DataPerEkin();
524  for (G4int idel=0; idel<gNumDelta; ++idel) {
525  DataPerDelta *perDelta = perEkin->fDataPerDelta[idel];
526  delete [] perDelta->fRejFuntion;
527  delete perDelta;
528  }
529  delete [] perEkin->fDataPerDelta;
530  delete perEkin;
531  }
532  delete [] data->fDataPerEkin;
533 }
534 
535 
537  for (size_t i=0; i<fMCDataPerElement.size(); ++i) {
538  if (fMCDataPerElement[i]) {
540  delete fMCDataPerElement[i];
541  }
542  }
543  fMCDataPerElement.clear();
544 }
545 
547  for (size_t i=0; i<fMCDataPerMaterial.size(); ++i) {
548  if (fMCDataPerMaterial[i]) {
550  delete fMCDataPerMaterial[i];
551  }
552  }
553  fMCDataPerMaterial.clear();
554 }
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static constexpr G4int gNumDelta
std::vector< DataPerMaterial * > fMCDataPerElement
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:600
static constexpr double cm3
static constexpr double g
static constexpr G4int gNumBeta2
size_t GetIndex() const
Definition: G4Material.hh:262
static constexpr G4double gMaxBeta2
static constexpr G4int gNumEkin
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:191
void LoadMCDataElement(const G4Element *)
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:207
static constexpr G4double gMaxDelta
static constexpr G4int gNumAngle
G4double G4Log(G4double x)
Definition: G4Log.hh:230
const XML_Char const XML_Char * data
Definition: expat.h:268
static constexpr double cm
Definition: SystemOfUnits.h:99
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
static constexpr G4int gMaxZet
static constexpr G4double gMinEkin
void GetMottCorrectionFactors(G4double logekin, G4double beta2, G4int matindx, G4double &mcToScr, G4double &mcToQ1, G4double &mcToG2PerG1)
static constexpr double electron_mass_c2
void InitMCDataMaterial(const G4Material *)
static constexpr double MeV
#define G4UniformRand()
Definition: Randomize.hh:53
static const std::string gElemSymbols[]
G4double GetMottRejectionValue(G4double logekin, G4double G4beta2, G4double q1, G4double cost, G4int matindx, G4int &ekindx, G4int &deltindx)
void AllocateDataPerMaterial(DataPerMaterial *)
Float_t mat
static G4ProductionCutsTable * GetProductionCutsTable()
static constexpr G4double gMidEkin
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
int G4lrint(double ad)
Definition: templates.hh:151
std::vector< G4Element * > G4ElementVector
int G4int
Definition: G4Types.hh:78
ifstream in
Definition: comparison.C:7
void DeAllocateDataPerMaterial(DataPerMaterial *)
void ReadCompressedFile(std::string fname, std::istringstream &iss)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
#define Z_OK
Definition: zlib.h:177
G4double GetZ() const
Definition: G4Element.hh:131
G4GSMottCorrection(G4bool iselectron=true)
const G4Material * GetMaterial() const
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:210
static constexpr double fine_structure_const
std::vector< DataPerMaterial * > fMCDataPerMaterial
size_t GetNumberOfElements() const
Definition: G4Material.hh:187
G4double GetDensity() const
Definition: G4Material.hh:181