ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GSPWACorrections.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4GSPWACorrections.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: 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 first moment correction.
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 
55 
56 const std::string G4GSPWACorrections::gElemSymbols[] = {"H","He","Li","Be","B" ,
57  "C" ,"N" ,"O" ,"F" ,"Ne","Na","Mg","Al","Si","P" , "S","Cl","Ar","K" ,"Ca","Sc",
58  "Ti","V" ,"Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr","Rb",
59  "Sr","Y" ,"Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I" ,
60  "Xe","Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm",
61  "Yb","Lu","Hf","Ta","W" ,"Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At",
62  "Rn","Fr","Ra","Ac","Th","Pa","U" ,"Np","Pu","Am","Cm","Bk","Cf"};
63 
64 G4GSPWACorrections::G4GSPWACorrections(G4bool iselectron) : fIsElectron(iselectron) {
65  // init grids related data member values
66  fMaxEkin = CLHEP::electron_mass_c2*(1./std::sqrt(1.-gMaxBeta2)-1.);
67  fLogMinEkin = G4Log(gMinEkin);
68  fInvLogDelEkin = (gNumEkin-gNumBeta2)/G4Log(gMidEkin/gMinEkin);
69  G4double pt2 = gMidEkin*(gMidEkin+2.0*CLHEP::electron_mass_c2);
70  fMinBeta2 = pt2/(pt2+CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
71  fInvDelBeta2 = (gNumBeta2-1.)/(gMaxBeta2-fMinBeta2);
72 }
73 
74 
75 G4GSPWACorrections::~G4GSPWACorrections() {
76  ClearDataPerElement();
77  ClearDataPerMaterial();
78 }
79 
80 
81 void G4GSPWACorrections::GetPWACorrectionFactors(G4double logekin, G4double beta2, G4int matindx,
82  G4double &corToScr, G4double &corToQ1, G4double &corToG2PerG1) {
83  G4int ekinIndxLow = 0;
84  G4double remRfaction = 0.;
85  if (beta2>=gMaxBeta2) {
86  ekinIndxLow = gNumEkin - 1;
87  // remRfaction = -1.
88  } else if (beta2>=fMinBeta2) { // linear interpolation on \beta^2
89  remRfaction = (beta2 - fMinBeta2) * fInvDelBeta2;
90  ekinIndxLow = (G4int)remRfaction;
91  remRfaction -= ekinIndxLow;
92  ekinIndxLow += (gNumEkin - gNumBeta2);
93  } else if (logekin>=fLogMinEkin) {
94  remRfaction = (logekin - fLogMinEkin) * fInvLogDelEkin;
95  ekinIndxLow = (G4int)remRfaction;
96  remRfaction -= ekinIndxLow;
97  } // the defaults otherwise i.e. use the lowest energy values when ekin is smaller than the minum ekin
98  //
99  DataPerMaterial *data = fDataPerMaterial[matindx];
100  corToScr = data->fCorScreening[ekinIndxLow];
101  corToQ1 = data->fCorFirstMoment[ekinIndxLow];
102  corToG2PerG1 = data->fCorSecondMoment[ekinIndxLow];
103  if (remRfaction>0.) {
104  corToScr += remRfaction*(data->fCorScreening[ekinIndxLow+1] - data->fCorScreening[ekinIndxLow]);
105  corToQ1 += remRfaction*(data->fCorFirstMoment[ekinIndxLow+1] - data->fCorFirstMoment[ekinIndxLow]);
106  corToG2PerG1 += remRfaction*(data->fCorSecondMoment[ekinIndxLow+1] - data->fCorSecondMoment[ekinIndxLow]);
107  }
108 }
109 
110 
111 void G4GSPWACorrections::Initialise() {
112  // load PWA correction data for each elements that belongs to materials that are used in the detector
113  InitDataPerElement();
114  // clear PWA correction data per material
115  ClearDataPerMaterial();
116  // initialise PWA correction data for the materials that are used in the detector
117  InitDataPerMaterials();
118 }
119 
120 
121 void G4GSPWACorrections::InitDataPerElement() {
122  // do it only once
123  if (fDataPerElement.size()<gMaxZet+1) {
124  fDataPerElement.resize(gMaxZet+1,nullptr);
125  }
126  // loop over all materials, for those that are used check the list of elements and load data from file if the
127  // corresponding data has not been loaded yet
128  G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable();
129  size_t numMatCuts = thePCTable->GetTableSize();
130  for (size_t imc=0; imc<numMatCuts; ++imc) {
131  const G4MaterialCutsCouple *matCut = thePCTable->GetMaterialCutsCouple(imc);
132  if (!matCut->IsUsed()) {
133  continue;
134  }
135  const G4Material *mat = matCut->GetMaterial();
136  const G4ElementVector *elemVect = mat->GetElementVector();
137  //
138  size_t numElems = elemVect->size();
139  for (size_t ielem=0; ielem<numElems; ++ielem) {
140  const G4Element *elem = (*elemVect)[ielem];
141  G4int izet = G4lrint(elem->GetZ());
142  if (izet>gMaxZet) {
143  izet = gMaxZet;
144  }
145  if (!fDataPerElement[izet]) {
146  LoadDataElement(elem);
147  }
148  }
149  }
150 }
151 
152 
153 void G4GSPWACorrections::InitDataPerMaterials() {
154  // prepare size of the container
155  size_t numMaterials = G4Material::GetNumberOfMaterials();
156  if (fDataPerMaterial.size()!=numMaterials) {
157  fDataPerMaterial.resize(numMaterials);
158  }
159  // init. PWA correction data for the Materials that are used in the geometry
160  G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable();
161  size_t numMatCuts = thePCTable->GetTableSize();
162  for (size_t imc=0; imc<numMatCuts; ++imc) {
163  const G4MaterialCutsCouple *matCut = thePCTable->GetMaterialCutsCouple(imc);
164  if (!matCut->IsUsed()) {
165  continue;
166  }
167  const G4Material *mat = matCut->GetMaterial();
168  if (!fDataPerMaterial[mat->GetIndex()]) {
169  InitDataMaterial(mat);
170  }
171  }
172 }
173 
174 
175 // it's called only if data has not been loaded for this element yet
176 void G4GSPWACorrections::LoadDataElement(const G4Element *elem) {
177  // allocate memory
178  G4int izet = elem->GetZasInt();
179  if (izet>gMaxZet) {
180  izet = gMaxZet;
181  }
182  // load data from file
183  char* tmppath = std::getenv("G4LEDATA");
184  if (!tmppath) {
185  G4Exception("G4GSPWACorrection::LoadDataElement()","em0006",
186  FatalException,
187  "Environment variable G4LEDATA not defined");
188  return;
189  }
190  std::string path(tmppath);
191  if (fIsElectron) {
192  path += "/msc_GS/PWACor/el/";
193  } else {
194  path += "/msc_GS/PWACor/pos/";
195  }
196  std::string fname = path+"cf_"+gElemSymbols[izet-1];
197  std::ifstream infile(fname,std::ios::in);
198  if (!infile.is_open()) {
199  std::string msg = " Problem while trying to read " + fname + " data file.\n";
200  G4Exception("G4GSPWACorrection::LoadDataElement","em0006", FatalException,msg.c_str());
201  return;
202  }
203  // allocate data structure
204  DataPerMaterial *perElem = new DataPerMaterial();
205  perElem->fCorScreening.resize(gNumEkin,0.0);
206  perElem->fCorFirstMoment.resize(gNumEkin,0.0);
207  perElem->fCorSecondMoment.resize(gNumEkin,0.0);
208  fDataPerElement[izet] = perElem;
209  G4double dum0;
210  for (G4int iek=0; iek<gNumEkin; ++iek) {
211  infile >> dum0;
212  infile >> perElem->fCorScreening[iek];
213  infile >> perElem->fCorFirstMoment[iek];
214  infile >> perElem->fCorSecondMoment[iek];
215  }
216  infile.close();
217 }
218 
219 
220 void G4GSPWACorrections::InitDataMaterial(const G4Material *mat) {
221  constexpr G4double const1 = 7821.6; // [cm2/g]
222  constexpr G4double const2 = 0.1569; // [cm2 MeV2 / g]
223  constexpr G4double finstrc2 = 5.325135453E-5; // fine-structure const. square
224 
225  G4double constFactor = CLHEP::electron_mass_c2*CLHEP::fine_structure_const/0.88534;
226  constFactor *= constFactor; // (mc^2)^2\alpha^2/( C_{TF}^2)
227  // allocate memory
228  DataPerMaterial *perMat = new DataPerMaterial();
229  perMat->fCorScreening.resize(gNumEkin,0.0);
230  perMat->fCorFirstMoment.resize(gNumEkin,0.0);
231  perMat->fCorSecondMoment.resize(gNumEkin,0.0);
232  fDataPerMaterial[mat->GetIndex()] = perMat;
233  //
234  const G4ElementVector* elemVect = mat->GetElementVector();
235  const G4int numElems = mat->GetNumberOfElements();
236  const G4double* nbAtomsPerVolVect = mat->GetVecNbOfAtomsPerVolume();
237  G4double totNbAtomsPerVol = mat->GetTotNbOfAtomsPerVolume();
238  //
239  // 1. Compute material dependent part of Moliere's b_c \chi_c^2
240  // (with \xi=1 (i.e. total sub-threshold scattering power correction)
241  G4double moliereBc = 0.0;
242  G4double moliereXc2 = 0.0;
243  G4double zs = 0.0;
244  G4double ze = 0.0;
245  G4double zx = 0.0;
246  G4double sa = 0.0;
247  G4double xi = 1.0;
248  for (G4int ielem=0; ielem<numElems; ++ielem) {
249  G4double zet = (*elemVect)[ielem]->GetZ();
250  if (zet>gMaxZet) {
251  zet = (G4double)gMaxZet;
252  }
253  G4double iwa = (*elemVect)[ielem]->GetN();
254  G4double ipz = nbAtomsPerVolVect[ielem]/totNbAtomsPerVol;
255  G4double dum = ipz*zet*(zet+xi);
256  zs += dum;
257  ze += dum*(-2.0/3.0)*G4Log(zet);
258  zx += dum*G4Log(1.0+3.34*finstrc2*zet*zet);
259  sa += ipz*iwa;
260  }
261  G4double density = mat->GetDensity()*CLHEP::cm3/CLHEP::g; // [g/cm3]
262  //
263  moliereBc = const1*density*zs/sa*G4Exp(ze/zs)/G4Exp(zx/zs); //[1/cm]
264  moliereXc2 = const2*density*zs/sa; // [MeV2/cm]
265  // change to Geant4 internal units of 1/length and energ2/length
266  moliereBc *= 1.0/CLHEP::cm;
267  moliereXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm;
268  //
269  // 2. loop over the kinetic energy grid
270  for (G4int iek=0; iek<gNumEkin; ++iek) {
271  // 2./a. set current kinetic energy and pt2 value
272  G4double ekin = G4Exp(fLogMinEkin+iek/fInvLogDelEkin);
273  G4double pt2 = ekin*(ekin+2.0*CLHEP::electron_mass_c2);
274  if (ekin>gMidEkin) {
275  G4double b2 = fMinBeta2+(iek-(gNumEkin-gNumBeta2))/fInvDelBeta2;
276  ekin = CLHEP::electron_mass_c2*(1./std::sqrt(1.-b2)-1.);
277  pt2 = ekin*(ekin+2.0*CLHEP::electron_mass_c2);
278  }
279  // 2./b. loop over the elements at the current kinetic energy point
280  for (G4int ielem=0; ielem<numElems; ++ielem) {
281  const G4Element *elem = (*elemVect)[ielem];
282  G4double zet = elem->GetZ();
283  if (zet>gMaxZet) {
284  zet = (G4double)gMaxZet;
285  }
286  G4int izet = G4lrint(zet);
287  // loaded PWA corrections for the current element
288  DataPerMaterial *perElem = fDataPerElement[izet];
289  //
290  // xi should be one i.e. z(z+1) since total sub-threshold scattering power correction
291  G4double nZZPlus1 = nbAtomsPerVolVect[ielem]*zet*(zet+1.0)/totNbAtomsPerVol;
292  G4double Z23 = std::pow(zet,2./3.);
293  //
294  // 2./b./(i) Add the 3 PWA correction factors
295  G4double mcScrCF = perElem->fCorScreening[iek]; // \kappa_i[1.13+3.76(\alpha Z_i)^2] with \kappa_i=scr_mc/scr_sr
296  // compute the screening parameter correction factor (Z_i contribution to the material)
297  // 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)}
298  // with C = \frac{(mc^2)^\alpha^2} {4(pc)^2 C_{TF}^2} = constFactor/(4*(pc)^2)
299  // 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
300  perMat->fCorScreening[iek] += nZZPlus1*G4Log(Z23*mcScrCF);
301  // compute the corrected screening parameter for the current Z_i and E_{kin}
302  // 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]
303  mcScrCF *= constFactor*Z23/(4.*pt2);
304  // compute first moment correction factor
305  // 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}
306  // where:
307  // 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
308  // 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)
309  // and \gamma_i = \sigma(Z_i)_{el}^(MC-DCS)/\sigma(Z_i,src(Z_i)_{mc})_{el}^(sr)
310  // 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
311  // A_i x B_i is stored in file per e-/e+, E_{kin} and Z_i
312  // here we compute the \sum_i n_i Z_i(Z_i+1) A_i B_i part
313  perMat->fCorFirstMoment[iek] += nZZPlus1*(G4Log(1.+1./mcScrCF)-1./(1.+mcScrCF))*perElem->fCorFirstMoment[iek];
314  // compute the second moment correction factor
315  // [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}
316  // 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}}
317  // here we compute the \sum_i n_i Z_i(Z_i+1) A_i part
318  perMat->fCorSecondMoment[iek] += nZZPlus1*perElem->fCorSecondMoment[iek];
319  //
320  // 2./b./(ii) When the last element has been added:
321  if (ielem==numElems-1) {
322  //
323  // 1. the remaining part of the sreening correction and divide the corrected screening par. with Moliere's one:
324  // (Moliere screening parameter = moliereXc2/(4(pc)^2 moliereBc) )
325  G4double dumScr = G4Exp(perMat->fCorScreening[iek]/zs);
326  perMat->fCorScreening[iek] = constFactor*dumScr*moliereBc/moliereXc2;
327  //
328  // 2. the remaining part of the first moment correction and divide by the one computed by using the corrected
329  // screening parameter (= (mc^2)^\alpha^2/(4(pc)^2C_{TF}^2) dumScr
330  G4double scrCorTed = constFactor*dumScr/(4.*pt2);
331  G4double dum0 = G4Log(1.+1./scrCorTed);
332  perMat->fCorFirstMoment[iek] = perMat->fCorFirstMoment[iek]/(zs*(dum0-1./(1.+scrCorTed)));
333  //
334  // 3. the remaining part of the second moment correction and divide by the one computed by using the corrected
335  // screening parameter
336  G4double G2PerG1 = 3.*(1.+scrCorTed)*((1.+2.*scrCorTed)*dum0-2.)/((1.+scrCorTed)*dum0-1.);
337  perMat->fCorSecondMoment[iek] = perMat->fCorSecondMoment[iek]/(zs*G2PerG1);
338  }
339  }
340  }
341 }
342 
343 
344 
345 void G4GSPWACorrections::ClearDataPerElement() {
346  for (size_t i=0; i<fDataPerElement.size(); ++i) {
347  if (fDataPerElement[i]) {
348  fDataPerElement[i]->fCorScreening.clear();
349  fDataPerElement[i]->fCorFirstMoment.clear();
350  fDataPerElement[i]->fCorSecondMoment.clear();
351  delete fDataPerElement[i];
352  }
353  }
354  fDataPerElement.clear();
355 }
356 
357 
358 void G4GSPWACorrections::ClearDataPerMaterial() {
359  for (size_t i=0; i<fDataPerMaterial.size(); ++i) {
360  if (fDataPerMaterial[i]) {
361  fDataPerMaterial[i]->fCorScreening.clear();
362  fDataPerMaterial[i]->fCorFirstMoment.clear();
363  fDataPerMaterial[i]->fCorSecondMoment.clear();
364  delete fDataPerMaterial[i];
365  }
366  }
367  fDataPerMaterial.clear();
368 }