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