Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // 27 // ------------------------------------------- 27 // ---------------------------------------------------------------------------- 28 // 28 // 29 // 29 // 30 // File name: G4GSPWACorrections 30 // File name: G4GSPWACorrections 31 // 31 // 32 // Author: Mihaly Novak 32 // Author: Mihaly Novak 33 // 33 // 34 // Creation date: 17.10.2017 34 // Creation date: 17.10.2017 35 // 35 // 36 // Modifications: 36 // Modifications: 37 // 02.02.2018 M.Novak: fixed initialization of 37 // 02.02.2018 M.Novak: fixed initialization of first moment correction. 38 // 38 // 39 // Class description: see the header file. 39 // Class description: see the header file. 40 // 40 // 41 // ------------------------------------------- 41 // ----------------------------------------------------------------------------- 42 42 43 #include "G4GSPWACorrections.hh" 43 #include "G4GSPWACorrections.hh" 44 44 45 #include "G4PhysicalConstants.hh" 45 #include "G4PhysicalConstants.hh" 46 #include "G4Log.hh" 46 #include "G4Log.hh" 47 #include "G4Exp.hh" 47 #include "G4Exp.hh" 48 48 49 #include "G4ProductionCutsTable.hh" 49 #include "G4ProductionCutsTable.hh" 50 #include "G4MaterialCutsCouple.hh" 50 #include "G4MaterialCutsCouple.hh" 51 #include "G4Material.hh" 51 #include "G4Material.hh" 52 #include "G4ElementVector.hh" 52 #include "G4ElementVector.hh" 53 #include "G4Element.hh" 53 #include "G4Element.hh" 54 #include "G4EmParameters.hh" << 55 54 56 55 57 const std::string G4GSPWACorrections::gElemSym 56 const std::string G4GSPWACorrections::gElemSymbols[] = {"H","He","Li","Be","B" , 58 "C" ,"N" ,"O" ,"F" ,"Ne","Na","Mg","Al","Si", 57 "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", 58 "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", 59 "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", 60 "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", 61 "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", 62 "Rn","Fr","Ra","Ac","Th","Pa","U" ,"Np","Pu","Am","Cm","Bk","Cf"}; 64 63 65 G4GSPWACorrections::G4GSPWACorrections(G4bool 64 G4GSPWACorrections::G4GSPWACorrections(G4bool iselectron) : fIsElectron(iselectron) { 66 // init grids related data member values 65 // init grids related data member values 67 fMaxEkin = CLHEP::electron_mass_c2*(1 66 fMaxEkin = CLHEP::electron_mass_c2*(1./std::sqrt(1.-gMaxBeta2)-1.); 68 fLogMinEkin = G4Log(gMinEkin); 67 fLogMinEkin = G4Log(gMinEkin); 69 fInvLogDelEkin = (gNumEkin-gNumBeta2)/G4Log 68 fInvLogDelEkin = (gNumEkin-gNumBeta2)/G4Log(gMidEkin/gMinEkin); 70 G4double pt2 = gMidEkin*(gMidEkin+2.0*CLH 69 G4double pt2 = gMidEkin*(gMidEkin+2.0*CLHEP::electron_mass_c2); 71 fMinBeta2 = pt2/(pt2+CLHEP::electron_m 70 fMinBeta2 = pt2/(pt2+CLHEP::electron_mass_c2*CLHEP::electron_mass_c2); 72 fInvDelBeta2 = (gNumBeta2-1.)/(gMaxBeta2- 71 fInvDelBeta2 = (gNumBeta2-1.)/(gMaxBeta2-fMinBeta2); 73 } 72 } 74 73 75 74 76 G4GSPWACorrections::~G4GSPWACorrections() { 75 G4GSPWACorrections::~G4GSPWACorrections() { 77 ClearDataPerElement(); 76 ClearDataPerElement(); 78 ClearDataPerMaterial(); 77 ClearDataPerMaterial(); 79 } 78 } 80 79 81 80 82 void G4GSPWACorrections::GetPWACorrectionFact 81 void G4GSPWACorrections::GetPWACorrectionFactors(G4double logekin, G4double beta2, G4int matindx, 83 82 G4double &corToScr, G4double &corToQ1, G4double &corToG2PerG1) { 84 G4int ekinIndxLow = 0; 83 G4int ekinIndxLow = 0; 85 G4double remRfaction = 0.; 84 G4double remRfaction = 0.; 86 if (beta2>=gMaxBeta2) { 85 if (beta2>=gMaxBeta2) { 87 ekinIndxLow = gNumEkin - 1; 86 ekinIndxLow = gNumEkin - 1; 88 // remRfaction = -1. 87 // remRfaction = -1. 89 } else if (beta2>=fMinBeta2) { // linear in 88 } else if (beta2>=fMinBeta2) { // linear interpolation on \beta^2 90 remRfaction = (beta2 - fMinBeta2) * fInv 89 remRfaction = (beta2 - fMinBeta2) * fInvDelBeta2; 91 ekinIndxLow = (G4int)remRfaction; 90 ekinIndxLow = (G4int)remRfaction; 92 remRfaction -= ekinIndxLow; 91 remRfaction -= ekinIndxLow; 93 ekinIndxLow += (gNumEkin - gNumBeta2); 92 ekinIndxLow += (gNumEkin - gNumBeta2); 94 } else if (logekin>=fLogMinEkin) { 93 } else if (logekin>=fLogMinEkin) { 95 remRfaction = (logekin - fLogMinEkin) * 94 remRfaction = (logekin - fLogMinEkin) * fInvLogDelEkin; 96 ekinIndxLow = (G4int)remRfaction; 95 ekinIndxLow = (G4int)remRfaction; 97 remRfaction -= ekinIndxLow; 96 remRfaction -= ekinIndxLow; 98 } // the defaults otherwise i.e. use the low 97 } // the defaults otherwise i.e. use the lowest energy values when ekin is smaller than the minum ekin 99 // 98 // 100 DataPerMaterial *data = fDataPerMaterial[mat 99 DataPerMaterial *data = fDataPerMaterial[matindx]; 101 corToScr = data->fCorScreening[ekinIndx 100 corToScr = data->fCorScreening[ekinIndxLow]; 102 corToQ1 = data->fCorFirstMoment[ekinIn 101 corToQ1 = data->fCorFirstMoment[ekinIndxLow]; 103 corToG2PerG1 = data->fCorSecondMoment[ekinI 102 corToG2PerG1 = data->fCorSecondMoment[ekinIndxLow]; 104 if (remRfaction>0.) { 103 if (remRfaction>0.) { 105 corToScr += remRfaction*(data->fCorSc 104 corToScr += remRfaction*(data->fCorScreening[ekinIndxLow+1] - data->fCorScreening[ekinIndxLow]); 106 corToQ1 += remRfaction*(data->fCorFi 105 corToQ1 += remRfaction*(data->fCorFirstMoment[ekinIndxLow+1] - data->fCorFirstMoment[ekinIndxLow]); 107 corToG2PerG1 += remRfaction*(data->fCorSe 106 corToG2PerG1 += remRfaction*(data->fCorSecondMoment[ekinIndxLow+1] - data->fCorSecondMoment[ekinIndxLow]); 108 } 107 } 109 } 108 } 110 109 111 110 112 void G4GSPWACorrections::Initialise() { 111 void G4GSPWACorrections::Initialise() { 113 // load PWA correction data for each element 112 // load PWA correction data for each elements that belongs to materials that are used in the detector 114 InitDataPerElement(); 113 InitDataPerElement(); 115 // clear PWA correction data per material 114 // clear PWA correction data per material 116 ClearDataPerMaterial(); 115 ClearDataPerMaterial(); 117 // initialise PWA correction data for the ma 116 // initialise PWA correction data for the materials that are used in the detector 118 InitDataPerMaterials(); 117 InitDataPerMaterials(); 119 } 118 } 120 119 121 120 122 void G4GSPWACorrections::InitDataPerElement() 121 void G4GSPWACorrections::InitDataPerElement() { 123 // do it only once 122 // do it only once 124 if (fDataPerElement.size()<gMaxZet+1) { 123 if (fDataPerElement.size()<gMaxZet+1) { 125 fDataPerElement.resize(gMaxZet+1,nullptr); 124 fDataPerElement.resize(gMaxZet+1,nullptr); 126 } 125 } 127 // loop over all materials, for those that a 126 // 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 ye 127 // corresponding data has not been loaded yet 129 G4ProductionCutsTable *thePCTable = G4Produc 128 G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable(); 130 G4int numMatCuts = (G4int)thePCTable->GetTab << 129 size_t numMatCuts = thePCTable->GetTableSize(); 131 for (G4int imc=0; imc<numMatCuts; ++imc) { << 130 for (size_t imc=0; imc<numMatCuts; ++imc) { 132 const G4MaterialCutsCouple *matCut = theP 131 const G4MaterialCutsCouple *matCut = thePCTable->GetMaterialCutsCouple(imc); 133 if (!matCut->IsUsed()) { 132 if (!matCut->IsUsed()) { 134 continue; 133 continue; 135 } 134 } 136 const G4Material *mat = matCut-> 135 const G4Material *mat = matCut->GetMaterial(); 137 const G4ElementVector *elemVect = mat->Get 136 const G4ElementVector *elemVect = mat->GetElementVector(); 138 // 137 // 139 std::size_t numElems = elemVect->size(); << 138 size_t numElems = elemVect->size(); 140 for (std::size_t ielem=0; ielem<numElems; << 139 for (size_t ielem=0; ielem<numElems; ++ielem) { 141 const G4Element *elem = (*elemVect)[iele 140 const G4Element *elem = (*elemVect)[ielem]; 142 G4int izet = G4lrint(elem->GetZ()); 141 G4int izet = G4lrint(elem->GetZ()); 143 if (izet>gMaxZet) { 142 if (izet>gMaxZet) { 144 izet = gMaxZet; 143 izet = gMaxZet; 145 } 144 } 146 if (!fDataPerElement[izet]) { 145 if (!fDataPerElement[izet]) { 147 LoadDataElement(elem); 146 LoadDataElement(elem); 148 } 147 } 149 } 148 } 150 } 149 } 151 } 150 } 152 151 153 152 154 void G4GSPWACorrections::InitDataPerMaterials( 153 void G4GSPWACorrections::InitDataPerMaterials() { 155 // prepare size of the container 154 // prepare size of the container 156 std::size_t numMaterials = G4Material::GetNu << 155 size_t numMaterials = G4Material::GetNumberOfMaterials(); 157 if (fDataPerMaterial.size()!=numMaterials) { 156 if (fDataPerMaterial.size()!=numMaterials) { 158 fDataPerMaterial.resize(numMaterials); 157 fDataPerMaterial.resize(numMaterials); 159 } 158 } 160 // init. PWA correction data for the Materia 159 // init. PWA correction data for the Materials that are used in the geometry 161 G4ProductionCutsTable *thePCTable = G4Produc 160 G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable(); 162 G4int numMatCuts = (G4int)thePCTable->GetTab << 161 size_t numMatCuts = thePCTable->GetTableSize(); 163 for (G4int imc=0; imc<numMatCuts; ++imc) { << 162 for (size_t imc=0; imc<numMatCuts; ++imc) { 164 const G4MaterialCutsCouple *matCut = theP 163 const G4MaterialCutsCouple *matCut = thePCTable->GetMaterialCutsCouple(imc); 165 if (!matCut->IsUsed()) { 164 if (!matCut->IsUsed()) { 166 continue; 165 continue; 167 } 166 } 168 const G4Material *mat = matCut->GetMateria 167 const G4Material *mat = matCut->GetMaterial(); 169 if (!fDataPerMaterial[mat->GetIndex()]) { 168 if (!fDataPerMaterial[mat->GetIndex()]) { 170 InitDataMaterial(mat); 169 InitDataMaterial(mat); 171 } 170 } 172 } 171 } 173 } 172 } 174 173 175 174 176 // it's called only if data has not been loade 175 // it's called only if data has not been loaded for this element yet 177 void G4GSPWACorrections::LoadDataElement(const 176 void G4GSPWACorrections::LoadDataElement(const G4Element *elem) { 178 // allocate memory 177 // allocate memory 179 G4int izet = elem->GetZasInt(); 178 G4int izet = elem->GetZasInt(); 180 if (izet>gMaxZet) { 179 if (izet>gMaxZet) { 181 izet = gMaxZet; 180 izet = gMaxZet; 182 } 181 } 183 // load data from file 182 // load data from file 184 std::string path = G4EmParameters::Instance( << 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); 185 if (fIsElectron) { 191 if (fIsElectron) { 186 path += "/msc_GS/PWACor/el/"; 192 path += "/msc_GS/PWACor/el/"; 187 } else { 193 } else { 188 path += "/msc_GS/PWACor/pos/"; 194 path += "/msc_GS/PWACor/pos/"; 189 } 195 } 190 std::string fname = path+"cf_"+gElemSymbols 196 std::string fname = path+"cf_"+gElemSymbols[izet-1]; 191 std::ifstream infile(fname,std::ios::in); 197 std::ifstream infile(fname,std::ios::in); 192 if (!infile.is_open()) { 198 if (!infile.is_open()) { 193 std::string msg = " Problem while trying 199 std::string msg = " Problem while trying to read " + fname + " data file.\n"; 194 G4Exception("G4GSPWACorrection::LoadDataEl 200 G4Exception("G4GSPWACorrection::LoadDataElement","em0006", FatalException,msg.c_str()); 195 return; 201 return; 196 } 202 } 197 // allocate data structure 203 // allocate data structure 198 auto perElem = new DataPerMaterial(); << 204 DataPerMaterial *perElem = new DataPerMaterial(); 199 perElem->fCorScreening.resize(gNumEkin,0.0); 205 perElem->fCorScreening.resize(gNumEkin,0.0); 200 perElem->fCorFirstMoment.resize(gNumEkin,0.0 206 perElem->fCorFirstMoment.resize(gNumEkin,0.0); 201 perElem->fCorSecondMoment.resize(gNumEkin,0. 207 perElem->fCorSecondMoment.resize(gNumEkin,0.0); 202 fDataPerElement[izet] = perElem; 208 fDataPerElement[izet] = perElem; 203 G4double dum0; 209 G4double dum0; 204 for (G4int iek=0; iek<gNumEkin; ++iek) { 210 for (G4int iek=0; iek<gNumEkin; ++iek) { 205 infile >> dum0; 211 infile >> dum0; 206 infile >> perElem->fCorScreening[iek]; 212 infile >> perElem->fCorScreening[iek]; 207 infile >> perElem->fCorFirstMoment[iek]; 213 infile >> perElem->fCorFirstMoment[iek]; 208 infile >> perElem->fCorSecondMoment[iek]; 214 infile >> perElem->fCorSecondMoment[iek]; 209 } 215 } 210 infile.close(); 216 infile.close(); 211 } 217 } 212 218 213 219 214 void G4GSPWACorrections::InitDataMaterial(cons 220 void G4GSPWACorrections::InitDataMaterial(const G4Material *mat) { 215 constexpr G4double const1 = 7821.6; / 221 constexpr G4double const1 = 7821.6; // [cm2/g] 216 constexpr G4double const2 = 0.1569; / 222 constexpr G4double const2 = 0.1569; // [cm2 MeV2 / g] 217 constexpr G4double finstrc2 = 5.325135453E-5 223 constexpr G4double finstrc2 = 5.325135453E-5; // fine-structure const. square 218 224 219 G4double constFactor = CLHEP::electro 225 G4double constFactor = CLHEP::electron_mass_c2*CLHEP::fine_structure_const/0.88534; 220 constFactor *= constFactor; 226 constFactor *= constFactor; // (mc^2)^2\alpha^2/( C_{TF}^2) 221 // allocate memory 227 // allocate memory 222 auto perMat = new DataPerMaterial(); << 228 DataPerMaterial *perMat = new DataPerMaterial(); 223 perMat->fCorScreening.resize(gNumEkin,0.0); 229 perMat->fCorScreening.resize(gNumEkin,0.0); 224 perMat->fCorFirstMoment.resize(gNumEkin,0.0) 230 perMat->fCorFirstMoment.resize(gNumEkin,0.0); 225 perMat->fCorSecondMoment.resize(gNumEkin,0.0 231 perMat->fCorSecondMoment.resize(gNumEkin,0.0); 226 fDataPerMaterial[mat->GetIndex()] = perMat; 232 fDataPerMaterial[mat->GetIndex()] = perMat; 227 // 233 // 228 const G4ElementVector* elemVect = 234 const G4ElementVector* elemVect = mat->GetElementVector(); 229 const G4int numElems = << 235 const G4int numElems = mat->GetNumberOfElements(); 230 const G4double* nbAtomsPerVolVect = 236 const G4double* nbAtomsPerVolVect = mat->GetVecNbOfAtomsPerVolume(); 231 G4double totNbAtomsPerVol = 237 G4double totNbAtomsPerVol = mat->GetTotNbOfAtomsPerVolume(); 232 // 238 // 233 // 1. Compute material dependent part of Mol 239 // 1. Compute material dependent part of Moliere's b_c \chi_c^2 234 // (with \xi=1 (i.e. total sub-threshold 240 // (with \xi=1 (i.e. total sub-threshold scattering power correction) 235 G4double moliereBc = 0.0; 241 G4double moliereBc = 0.0; 236 G4double moliereXc2 = 0.0; 242 G4double moliereXc2 = 0.0; 237 G4double zs = 0.0; 243 G4double zs = 0.0; 238 G4double ze = 0.0; 244 G4double ze = 0.0; 239 G4double zx = 0.0; 245 G4double zx = 0.0; 240 G4double sa = 0.0; 246 G4double sa = 0.0; 241 G4double xi = 1.0; 247 G4double xi = 1.0; 242 for (G4int ielem=0; ielem<numElems; ++ielem) 248 for (G4int ielem=0; ielem<numElems; ++ielem) { 243 G4double zet = (*elemVect)[ielem]->GetZ(); 249 G4double zet = (*elemVect)[ielem]->GetZ(); 244 if (zet>gMaxZet) { 250 if (zet>gMaxZet) { 245 zet = (G4double)gMaxZet; 251 zet = (G4double)gMaxZet; 246 } 252 } 247 G4double iwa = (*elemVect)[ielem]->GetN() 253 G4double iwa = (*elemVect)[ielem]->GetN(); 248 G4double ipz = nbAtomsPerVolVect[ielem]/t 254 G4double ipz = nbAtomsPerVolVect[ielem]/totNbAtomsPerVol; 249 G4double dum = ipz*zet*(zet+xi); 255 G4double dum = ipz*zet*(zet+xi); 250 zs += dum; 256 zs += dum; 251 ze += dum*(-2.0/3.0)*G4Log(zet); 257 ze += dum*(-2.0/3.0)*G4Log(zet); 252 zx += dum*G4Log(1.0+3.34*finstrc 258 zx += dum*G4Log(1.0+3.34*finstrc2*zet*zet); 253 sa += ipz*iwa; 259 sa += ipz*iwa; 254 } 260 } 255 G4double density = mat->GetDensity()*CLHEP:: 261 G4double density = mat->GetDensity()*CLHEP::cm3/CLHEP::g; // [g/cm3] 256 // 262 // 257 G4double z0 = (0.0 == sa) ? 0.0 : zs/sa; << 263 moliereBc = const1*density*zs/sa*G4Exp(ze/zs)/G4Exp(zx/zs); //[1/cm] 258 G4double z1 = (0.0 == zs) ? 0.0 : (ze - zx)/ << 264 moliereXc2 = const2*density*zs/sa; // [MeV2/cm] 259 moliereBc = const1*density*z0*G4Exp(z1); / << 260 moliereXc2 = const2*density*z0; // [MeV2/cm << 261 // change to Geant4 internal units of 1/leng 265 // change to Geant4 internal units of 1/length and energ2/length 262 moliereBc *= 1.0/CLHEP::cm; 266 moliereBc *= 1.0/CLHEP::cm; 263 moliereXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::c 267 moliereXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm; 264 // 268 // 265 // 2. loop over the kinetic energy grid 269 // 2. loop over the kinetic energy grid 266 for (G4int iek=0; iek<gNumEkin; ++iek) { 270 for (G4int iek=0; iek<gNumEkin; ++iek) { 267 // 2./a. set current kinetic energy and pt 271 // 2./a. set current kinetic energy and pt2 value 268 G4double ekin = G4Exp(fLogMinEk 272 G4double ekin = G4Exp(fLogMinEkin+iek/fInvLogDelEkin); 269 G4double pt2 = ekin*(ekin+2.0*CLHEP::el 273 G4double pt2 = ekin*(ekin+2.0*CLHEP::electron_mass_c2); 270 if (ekin>gMidEkin) { 274 if (ekin>gMidEkin) { 271 G4double b2 = fMinBeta2+(iek-(gNumEk 275 G4double b2 = fMinBeta2+(iek-(gNumEkin-gNumBeta2))/fInvDelBeta2; 272 ekin = CLHEP::electron_mass_c2*(1./std 276 ekin = CLHEP::electron_mass_c2*(1./std::sqrt(1.-b2)-1.); 273 pt2 = ekin*(ekin+2.0*CLHEP::electron_ 277 pt2 = ekin*(ekin+2.0*CLHEP::electron_mass_c2); 274 } 278 } 275 // 2./b. loop over the elements at the cur 279 // 2./b. loop over the elements at the current kinetic energy point 276 for (G4int ielem=0; ielem<numElems; ++iele 280 for (G4int ielem=0; ielem<numElems; ++ielem) { 277 const G4Element *elem = (*elemVect)[iele 281 const G4Element *elem = (*elemVect)[ielem]; 278 G4double zet = elem->GetZ(); 282 G4double zet = elem->GetZ(); 279 if (zet>gMaxZet) { 283 if (zet>gMaxZet) { 280 zet = (G4double)gMaxZet; 284 zet = (G4double)gMaxZet; 281 } 285 } 282 G4int izet = G4lrint(zet); 286 G4int izet = G4lrint(zet); 283 // loaded PWA corrections for the curren 287 // loaded PWA corrections for the current element 284 DataPerMaterial *perElem = fDataPerElem 288 DataPerMaterial *perElem = fDataPerElement[izet]; 285 // 289 // 286 // xi should be one i.e. z(z+1) since to 290 // xi should be one i.e. z(z+1) since total sub-threshold scattering power correction 287 G4double nZZPlus1 = nbAtomsPerVolVect[i 291 G4double nZZPlus1 = nbAtomsPerVolVect[ielem]*zet*(zet+1.0)/totNbAtomsPerVol; 288 G4double Z23 = std::pow(zet,2./3.) 292 G4double Z23 = std::pow(zet,2./3.); 289 // 293 // 290 // 2./b./(i) Add the 3 PWA correction fa 294 // 2./b./(i) Add the 3 PWA correction factors 291 G4double mcScrCF = perElem->fCorScreenin 295 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 corre 296 // compute the screening parameter correction factor (Z_i contribution to the material) 293 // src_{mc} = C \exp\left[ \frac{ \sum_i 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)} 294 // with C = \frac{(mc^2)^\alpha^2} {4(pc 298 // 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_ 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 296 perMat->fCorScreening[iek] += nZZPlus1*G 300 perMat->fCorScreening[iek] += nZZPlus1*G4Log(Z23*mcScrCF); 297 // compute the corrected screening param 301 // compute the corrected screening parameter for the current Z_i and E_{kin} 298 // src(Z_i)_{mc} = \frac{(mc^2)^\alpha^2 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] 299 mcScrCF *= constFactor*Z23/(4.*pt2); 303 mcScrCF *= constFactor*Z23/(4.*pt2); 300 // compute first moment correction facto 304 // compute first moment correction factor 301 // q1_{mc} = \frac{ \sum_i n_i Z_i(Z_i+1 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} 302 // where: 306 // where: 303 // A_i(src(Z_i)_{mc}) = [\ln(1+1/src(Z_i 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 304 // B_i = \beta_i \gamma_i with beta_i(Z_ 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) 305 // and \gamma_i = \sigma(Z_i)_{el}^(MC-D 309 // 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/ 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 307 // A_i x B_i is stored in file per e-/e+ 311 // 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_ 312 // here we compute the \sum_i n_i Z_i(Z_i+1) A_i B_i part 309 perMat->fCorFirstMoment[iek] += nZZPlus1 313 perMat->fCorFirstMoment[iek] += nZZPlus1*(G4Log(1.+1./mcScrCF)-1./(1.+mcScrCF))*perElem->fCorFirstMoment[iek]; 310 // compute the second moment correction 314 // compute the second moment correction factor 311 // [G2/G1]_{mc} = \frac{ \sum_i n_i Z_i( 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} 312 // with A_i(Z_i) = G2(Z_i)^{PWA}/G1(Z_i) 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}} 313 // here we compute the \sum_i n_i Z_i(Z_ 317 // here we compute the \sum_i n_i Z_i(Z_i+1) A_i part 314 perMat->fCorSecondMoment[iek] += nZZPlus 318 perMat->fCorSecondMoment[iek] += nZZPlus1*perElem->fCorSecondMoment[iek]; 315 // 319 // 316 // 2./b./(ii) When the last element has 320 // 2./b./(ii) When the last element has been added: 317 if (ielem==numElems-1) { 321 if (ielem==numElems-1) { 318 // 322 // 319 // 1. the remaining part of the sreeni 323 // 1. the remaining part of the sreening correction and divide the corrected screening par. with Moliere's one: 320 // (Moliere screening parameter = m 324 // (Moliere screening parameter = moliereXc2/(4(pc)^2 moliereBc) ) 321 G4double dumScr = G4Exp(perMat->fCor 325 G4double dumScr = G4Exp(perMat->fCorScreening[iek]/zs); 322 perMat->fCorScreening[iek] = constFact 326 perMat->fCorScreening[iek] = constFactor*dumScr*moliereBc/moliereXc2; 323 // 327 // 324 // 2. the remaining part of the first 328 // 2. the remaining part of the first moment correction and divide by the one computed by using the corrected 325 // screening parameter (= (mc^2)^\a 329 // screening parameter (= (mc^2)^\alpha^2/(4(pc)^2C_{TF}^2) dumScr 326 G4double scrCorTed = constFactor*dumSc 330 G4double scrCorTed = constFactor*dumScr/(4.*pt2); 327 G4double dum0 = G4Log(1.+1./scrCo 331 G4double dum0 = G4Log(1.+1./scrCorTed); 328 perMat->fCorFirstMoment[iek] = perMat- 332 perMat->fCorFirstMoment[iek] = perMat->fCorFirstMoment[iek]/(zs*(dum0-1./(1.+scrCorTed))); 329 // 333 // 330 // 3. the remaining part of the second 334 // 3. the remaining part of the second moment correction and divide by the one computed by using the corrected 331 // screening parameter 335 // screening parameter 332 G4double G2PerG1 = 3.*(1.+scrCorTed 336 G4double G2PerG1 = 3.*(1.+scrCorTed)*((1.+2.*scrCorTed)*dum0-2.)/((1.+scrCorTed)*dum0-1.); 333 perMat->fCorSecondMoment[iek] = perMat 337 perMat->fCorSecondMoment[iek] = perMat->fCorSecondMoment[iek]/(zs*G2PerG1); 334 } 338 } 335 } 339 } 336 } 340 } 337 } 341 } 338 342 339 343 340 344 341 void G4GSPWACorrections::ClearDataPerElement() 345 void G4GSPWACorrections::ClearDataPerElement() { 342 for (std::size_t i=0; i<fDataPerElement.size << 346 for (size_t i=0; i<fDataPerElement.size(); ++i) { 343 if (fDataPerElement[i]) { 347 if (fDataPerElement[i]) { 344 fDataPerElement[i]->fCorScreening.clear( 348 fDataPerElement[i]->fCorScreening.clear(); 345 fDataPerElement[i]->fCorFirstMoment.clea 349 fDataPerElement[i]->fCorFirstMoment.clear(); 346 fDataPerElement[i]->fCorSecondMoment.cle 350 fDataPerElement[i]->fCorSecondMoment.clear(); 347 delete fDataPerElement[i]; 351 delete fDataPerElement[i]; 348 } 352 } 349 } 353 } 350 fDataPerElement.clear(); 354 fDataPerElement.clear(); 351 } 355 } 352 356 353 357 354 void G4GSPWACorrections::ClearDataPerMaterial( 358 void G4GSPWACorrections::ClearDataPerMaterial() { 355 for (std::size_t i=0; i<fDataPerMaterial.siz << 359 for (size_t i=0; i<fDataPerMaterial.size(); ++i) { 356 if (fDataPerMaterial[i]) { 360 if (fDataPerMaterial[i]) { 357 fDataPerMaterial[i]->fCorScreening.clear 361 fDataPerMaterial[i]->fCorScreening.clear(); 358 fDataPerMaterial[i]->fCorFirstMoment.cle 362 fDataPerMaterial[i]->fCorFirstMoment.clear(); 359 fDataPerMaterial[i]->fCorSecondMoment.cl 363 fDataPerMaterial[i]->fCorSecondMoment.clear(); 360 delete fDataPerMaterial[i]; 364 delete fDataPerMaterial[i]; 361 } 365 } 362 } 366 } 363 fDataPerMaterial.clear(); 367 fDataPerMaterial.clear(); 364 } 368 } 365 369