Geant4 Cross Reference

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