Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4GSMottCorrection.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 #include "G4EmParameters.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.);
 75   fLogMinEkin     = G4Log(gMinEkin);
 76   fInvLogDelEkin  = (gNumEkin-gNumBeta2)/G4Log(gMidEkin/gMinEkin);
 77   G4double pt2    = gMidEkin*(gMidEkin+2.0*CLHEP::electron_mass_c2);
 78   fMinBeta2       = pt2/(pt2+CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
 79   fInvDelBeta2    = (gNumBeta2-1.)/(gMaxBeta2-fMinBeta2);
 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::GetMottCorrectionFactors(G4double logekin, G4double beta2, G4int matindx, G4double &mcToScr,
 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
123 double G4GSMottCorrection::GetMottRejectionValue(G4double logekin, G4double beta2, G4double q1, G4double cost,
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 
190 void G4GSMottCorrection::Initialise() {
191   // load Mott-correction data for each elements that belongs to materials that are used in the detector
192   InitMCDataPerElement();
193   // clrea Mott-correction data per material
194   ClearMCDataPerMaterial();
195   // initialise Mott-correction data for the materials that are used in the detector
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 are used check the list of elements and load data from file if the
206   // corresponding data has not been loaded yet
207   G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable();
208   G4int numMatCuts = (G4int)thePCTable->GetTableSize();
209   for (G4int 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     std::size_t numElems = elemVect->size();
218     for (std::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 
232 void G4GSMottCorrection::InitMCDataPerMaterials() {
233   // prepare size of the container
234   std::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
239   G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable();
240   G4int numMatCuts = (G4int)thePCTable->GetTableSize();
241   for (G4int 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
255 void G4GSMottCorrection::LoadMCDataElement(const G4Element *elem) {
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()->GetDirLEDATA();
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[izet-1];
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->fDataPerEkin[iek];
278     // 1. get the 3 Mott-correction factors for the current kinetic energy
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->fDataPerDelta[idel];
285       // 2./a.  : first the rejection function values
286       for (G4int iang=0; iang<gNumAngle; ++iang) {
287         infile >> perDelta->fRejFuntion[iang];
288       }
289       // 2./b. : then the 4 spline parameter for the last bin
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 string stream
299 void G4GSMottCorrection::ReadCompressedFile(const std::string& fname, std::istringstream &iss) {
300   std::string *dataString = nullptr;
301   std::string compfilename(fname+".z");
302   // create input stream with binary mode operation and positioning at the end of the file
303   std::ifstream in(compfilename, std::ios::binary | std::ios::ate);
304   if (in.good()) {
305      // get current position in the stream (was set to the end)
306      std::streamoff fileSize = in.tellg();
307      // set current position being the beginning of the stream
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 uncompressed data
315      uLongf complen    = (uLongf)(fileSize*4);
316      Bytef *uncompdata = new Bytef[complen];
317      while (Z_OK!=uncompress(uncompdata, &complen, compdata, fileSize)) {
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 data (will be deallocated by the caller)
326      dataString = new std::string((char*)uncompdata, (long)complen);
327      // delete the uncompressed data buffer
328      delete [] uncompdata;
329   } else {
330     std::string msg = "  Problem while trying to read " + compfilename + " data file.\n";
331     G4Exception("G4GSMottCorrection::ReadCompressedFile","em0006", FatalException,msg.c_str());
332     return;
333   }
334   // create the input string stream from the data string
335   if (dataString) {
336     iss.str(*dataString);
337     in.close();
338     delete dataString;
339   }
340 }
341 
342 
343 void G4GSMottCorrection::InitMCDataMaterial(const G4Material *mat) {
344   constexpr G4double const1   = 7821.6;      // [cm2/g]
345   constexpr G4double const2   = 0.1569;      // [cm2 MeV2 / g]
346   constexpr G4double finstrc2 = 5.325135453E-5; // fine-structure const. square
347 
348   G4double constFactor        = CLHEP::electron_mass_c2*CLHEP::fine_structure_const/0.88534;
349   constFactor                *= constFactor;  // (mc^2)^2\alpha^2/( C_{TF}^2)
350   // allocate memory
351   auto perMat = new DataPerMaterial();
352   AllocateDataPerMaterial(perMat);
353   fMCDataPerMaterial[mat->GetIndex()] = perMat;
354   //
355   const G4ElementVector* elemVect           = mat->GetElementVector();
356   const G4int            numElems           = (G4int)mat->GetNumberOfElements();
357   const G4double*        nbAtomsPerVolVect  = mat->GetVecNbOfAtomsPerVolume();
358   G4double               totNbAtomsPerVol   = mat->GetTotNbOfAtomsPerVolume();
359   //
360   // 1. Compute material dependent part of Moliere's b_c \chi_c^2
361   //    (with \xi=1 (i.e. total sub-threshold scattering power correction)
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]/totNbAtomsPerVol;
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*finstrc2*zet*zet);
380     sa           += ipz*iwa;
381   }
382   G4double density = mat->GetDensity()*CLHEP::cm3/CLHEP::g; // [g/cm3]
383   //
384   G4double z0 = (0.0 == sa) ? 0.0 : zs/sa;
385   G4double z1 = (0.0 == zs) ? 0.0 : (ze - zx)/zs;
386 
387   moliereBc  = const1*density*z0*G4Exp(z1);  //[1/cm]
388   moliereXc2 = const2*density*z0;  // [MeV2/cm]
389   // change to Geant4 internal units of 1/length and energ2/length
390   moliereBc  *= 1.0/CLHEP::cm;
391   moliereXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm;
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 pt2 value
396       G4double ekin          = G4Exp(fLogMinEkin+iek/fInvLogDelEkin);
397       G4double pt2  = ekin*(ekin+2.0*CLHEP::electron_mass_c2);
398       if (ekin>gMidEkin) {
399         G4double b2   = fMinBeta2+(iek-(gNumEkin-gNumBeta2))/fInvDelBeta2;
400         ekin = CLHEP::electron_mass_c2*(1./std::sqrt(1.-b2)-1.);
401         pt2  = ekin*(ekin+2.0*CLHEP::electron_mass_c2);
402       }
403     // 2./b. loop over the elements at the current kinetic energy point
404     for (G4int ielem=0; ielem<numElems; ++ielem) {
405       const G4Element *elem = (*elemVect)[ielem];
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 total sub-threshold scattering power correction
412       G4double nZZPlus1  = nbAtomsPerVolVect[ielem]*zet*(zet+1.0)/totNbAtomsPerVol;
413       G4double Z23       = std::pow(zet,2./3.);
414       //
415       DataPerEkin *perElemPerEkin  = fMCDataPerElement[izet]->fDataPerEkin[iek];
416       DataPerEkin *perMatPerEkin   = perMat->fDataPerEkin[iek];
417       //
418       // 2./b./(i) Add the 3 Mott-correction factors
419       G4double mcScrCF = perElemPerEkin->fMCScreening;     // \kappa_i[1.13+3.76(\alpha Z_i)^2] with \kappa_i=scr_mc/scr_sr
420       // compute the screening parameter correction factor (Z_i contribution to the material)
421       // 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)}
422       // with C = \frac{(mc^2)^\alpha^2} {4(pc)^2 C_{TF}^2} = constFactor/(4*(pc)^2)
423       // 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
424       perMatPerEkin->fMCScreening += nZZPlus1*G4Log(Z23*mcScrCF);
425       // compute the corrected screening parameter for the current Z_i and E_{kin}
426       // 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]
427       mcScrCF *= constFactor*Z23/(4.*pt2);
428       // compute first moment correction factor
429       // 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}
430       // where:
431       // 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
432       // 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)
433       // and \gamma_i = \sigma(Z_i)_{el}^(MC-DCS)/\sigma(Z_i,src(Z_i)_{mc})_{el}^(sr)
434       // 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
435       // A_i x B_i is stored in file per e-/e+, E_{kin} and Z_i
436       // here we compute the \sum_i n_i Z_i(Z_i+1) A_i  B_i part
437       perMatPerEkin->fMCFirstMoment += nZZPlus1*(G4Log(1.+1./mcScrCF)-1./(1.+mcScrCF))*perElemPerEkin->fMCFirstMoment;
438       // compute the second moment correction factor
439       // [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}
440       // 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}}
441       // here we compute the \sum_i n_i Z_i(Z_i+1) A_i part
442       perMatPerEkin->fMCSecondMoment += nZZPlus1*perElemPerEkin->fMCSecondMoment;
443       //
444       // 2./b./(ii) Go for the rejection funtion part
445       // I. loop over delta values
446       for (G4int idel=0; idel<gNumDelta; ++idel) {
447         DataPerDelta *perMatPerDelta  = perMatPerEkin->fDataPerDelta[idel];
448         DataPerDelta *perElemPerDelta = perElemPerEkin->fDataPerDelta[idel];
449         // I./a. loop over angles (i.e. the \sin(0.5\theta) values) and add the rejection function
450         for (G4int iang=0; iang<gNumAngle; ++iang) {
451           perMatPerDelta->fRejFuntion[iang] += nZZPlus1*perElemPerDelta->fRejFuntion[iang];
452         }
453         // I./b. get the last bin spline parameters and add them (a+bx+cx^2+dx^3)
454         perMatPerDelta->fSA += nZZPlus1*perElemPerDelta->fSA;
455         perMatPerDelta->fSB += nZZPlus1*perElemPerDelta->fSB;
456         perMatPerDelta->fSC += nZZPlus1*perElemPerDelta->fSC;
457         perMatPerDelta->fSD += nZZPlus1*perElemPerDelta->fSD;
458       }
459       //
460       // 2./b./(iii) When the last element has been added:
461       if (ielem==numElems-1) {
462         //
463         // 1. the remaining part of the sreening correction and divide the corrected screening par. with Moliere's one:
464         //    (Moliere screening parameter = moliereXc2/(4(pc)^2 moliereBc) )
465         G4double dumScr   = G4Exp(perMatPerEkin->fMCScreening/zs);
466         perMatPerEkin->fMCScreening = constFactor*dumScr*moliereBc/moliereXc2;
467         //
468         // 2. the remaining part of the first moment correction and divide by the one computed by using the corrected
469         //    screening parameter (= (mc^2)^\alpha^2/(4(pc)^2C_{TF}^2) dumScr
470         G4double scrCorTed = constFactor*dumScr/(4.*pt2);
471         G4double dum0      = G4Log(1.+1./scrCorTed);
472         perMatPerEkin->fMCFirstMoment = perMatPerEkin->fMCFirstMoment/(zs*(dum0-1./(1.+scrCorTed)));
473         //
474         // 3. the remaining part of the second moment correction and divide by the one computed by using the corrected
475         //    screening parameter
476         G4double G2PerG1   =  3.*(1.+scrCorTed)*((1.+2.*scrCorTed)*dum0-2.)/((1.+scrCorTed)*dum0-1.);
477         perMatPerEkin->fMCSecondMoment = perMatPerEkin->fMCSecondMoment/(zs*G2PerG1);
478         //
479         // 4. scale the maximum of the rejection function to unity and correct the last bin spline parameters as well
480         // I. loop over delta values
481         for (G4int idel=0; idel<gNumDelta; ++idel) {
482           DataPerDelta *perMatPerDelta  = perMatPerEkin->fDataPerDelta[idel];
483           G4double maxVal = -1.;
484           // II. llop over angles
485           for (G4int iang=0; iang<gNumAngle; ++iang) {
486             if (perMatPerDelta->fRejFuntion[iang]>maxVal)
487               maxVal = perMatPerDelta->fRejFuntion[iang];
488           }
489           for (G4int iang=0; iang<gNumAngle; ++iang) {
490             perMatPerDelta->fRejFuntion[iang] /=maxVal;
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::AllocateDataPerMaterial(DataPerMaterial *data) {
504   data->fDataPerEkin = new DataPerEkin*[gNumEkin]();
505   for (G4int iek=0; iek<gNumEkin; ++iek) {
506     auto perEkin = new DataPerEkin();
507     perEkin->fDataPerDelta = new DataPerDelta*[gNumDelta]();
508     for (G4int idel=0; idel<gNumDelta; ++idel) {
509       auto perDelta                = new DataPerDelta();
510       perDelta->fRejFuntion        = new double[gNumAngle]();
511       perEkin->fDataPerDelta[idel] = perDelta;
512     }
513     data->fDataPerEkin[iek] = perEkin;
514   }
515 }
516 
517 void G4GSMottCorrection::DeAllocateDataPerMaterial(DataPerMaterial *data) {
518   for (G4int iek=0; iek<gNumEkin; ++iek) {
519     DataPerEkin *perEkin = data->fDataPerEkin[iek]; //new DataPerEkin();
520     for (G4int idel=0; idel<gNumDelta; ++idel) {
521       DataPerDelta *perDelta = perEkin->fDataPerDelta[idel];
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.size(); ++i) {
534     if (fMCDataPerElement[i]) {
535       DeAllocateDataPerMaterial(fMCDataPerElement[i]);
536       delete fMCDataPerElement[i];
537     }
538   }
539   fMCDataPerElement.clear();
540 }
541 
542 void G4GSMottCorrection::ClearMCDataPerMaterial() {
543   for (std::size_t i=0; i<fMCDataPerMaterial.size(); ++i) {
544     if (fMCDataPerMaterial[i]) {
545       DeAllocateDataPerMaterial(fMCDataPerMaterial[i]);
546       delete fMCDataPerMaterial[i];
547     }
548   }
549   fMCDataPerMaterial.clear();
550 }
551