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: G4GSMottCorrection 30 // File name: G4GSMottCorrection 31 // 31 // 32 // Author: Mihaly Novak 32 // Author: Mihaly Novak 33 // 33 // 34 // Creation date: 23.08.2017 34 // Creation date: 23.08.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 "G4GSMottCorrection.hh" 43 #include "G4GSMottCorrection.hh" 44 44 45 #include "G4PhysicalConstants.hh" 45 #include "G4PhysicalConstants.hh" 46 #include "zlib.h" 46 #include "zlib.h" 47 #include "Randomize.hh" 47 #include "Randomize.hh" 48 #include "G4Log.hh" 48 #include "G4Log.hh" 49 #include "G4Exp.hh" 49 #include "G4Exp.hh" 50 50 51 #include "G4ProductionCutsTable.hh" 51 #include "G4ProductionCutsTable.hh" 52 #include "G4MaterialCutsCouple.hh" 52 #include "G4MaterialCutsCouple.hh" 53 #include "G4Material.hh" 53 #include "G4Material.hh" 54 #include "G4ElementVector.hh" 54 #include "G4ElementVector.hh" 55 #include "G4Element.hh" 55 #include "G4Element.hh" 56 #include "G4EmParameters.hh" << 57 56 58 #include <iostream> 57 #include <iostream> 59 #include <fstream> 58 #include <fstream> 60 #include <cmath> 59 #include <cmath> 61 #include <algorithm> 60 #include <algorithm> 62 61 63 62 64 const std::string G4GSMottCorrection::gElemSym 63 const std::string G4GSMottCorrection::gElemSymbols[] = {"H","He","Li","Be","B" , 65 "C" ,"N" ,"O" ,"F" ,"Ne","Na","Mg","Al","Si", 64 "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", 65 "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", 66 "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", 67 "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", 68 "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", 69 "Rn","Fr","Ra","Ac","Th","Pa","U" ,"Np","Pu","Am","Cm","Bk","Cf"}; 71 70 72 G4GSMottCorrection::G4GSMottCorrection(G4bool 71 G4GSMottCorrection::G4GSMottCorrection(G4bool iselectron) : fIsElectron(iselectron) { 73 // init grids related data member values 72 // init grids related data member values 74 fMaxEkin = CLHEP::electron_mass_c2*(1 73 fMaxEkin = CLHEP::electron_mass_c2*(1./std::sqrt(1.-gMaxBeta2)-1.); 75 fLogMinEkin = G4Log(gMinEkin); 74 fLogMinEkin = G4Log(gMinEkin); 76 fInvLogDelEkin = (gNumEkin-gNumBeta2)/G4Log 75 fInvLogDelEkin = (gNumEkin-gNumBeta2)/G4Log(gMidEkin/gMinEkin); 77 G4double pt2 = gMidEkin*(gMidEkin+2.0*CLH 76 G4double pt2 = gMidEkin*(gMidEkin+2.0*CLHEP::electron_mass_c2); 78 fMinBeta2 = pt2/(pt2+CLHEP::electron_m 77 fMinBeta2 = pt2/(pt2+CLHEP::electron_mass_c2*CLHEP::electron_mass_c2); 79 fInvDelBeta2 = (gNumBeta2-1.)/(gMaxBeta2- 78 fInvDelBeta2 = (gNumBeta2-1.)/(gMaxBeta2-fMinBeta2); 80 fInvDelDelta = (gNumDelta-1.)/gMaxDelta; 79 fInvDelDelta = (gNumDelta-1.)/gMaxDelta; 81 fInvDelAngle = gNumAngle-1.; 80 fInvDelAngle = gNumAngle-1.; 82 } 81 } 83 82 84 83 85 G4GSMottCorrection::~G4GSMottCorrection() { 84 G4GSMottCorrection::~G4GSMottCorrection() { 86 ClearMCDataPerElement(); 85 ClearMCDataPerElement(); 87 ClearMCDataPerMaterial(); 86 ClearMCDataPerMaterial(); 88 } 87 } 89 88 90 89 91 void G4GSMottCorrection::GetMottCorrectionFact 90 void G4GSMottCorrection::GetMottCorrectionFactors(G4double logekin, G4double beta2, G4int matindx, G4double &mcToScr, 92 91 G4double &mcToQ1, G4double &mcToG2PerG1) { 93 G4int ekinIndxLow = 0; 92 G4int ekinIndxLow = 0; 94 G4double remRfaction = 0.; 93 G4double remRfaction = 0.; 95 if (beta2>=gMaxBeta2) { 94 if (beta2>=gMaxBeta2) { 96 ekinIndxLow = gNumEkin - 1; 95 ekinIndxLow = gNumEkin - 1; 97 // remRfaction = -1. 96 // remRfaction = -1. 98 } else if (beta2>=fMinBeta2) { // linear in 97 } else if (beta2>=fMinBeta2) { // linear interpolation on \beta^2 99 remRfaction = (beta2 - fMinBeta2) * fInv 98 remRfaction = (beta2 - fMinBeta2) * fInvDelBeta2; 100 ekinIndxLow = (G4int)remRfaction; 99 ekinIndxLow = (G4int)remRfaction; 101 remRfaction -= ekinIndxLow; 100 remRfaction -= ekinIndxLow; 102 ekinIndxLow += (gNumEkin - gNumBeta2); 101 ekinIndxLow += (gNumEkin - gNumBeta2); 103 } else if (logekin>=fLogMinEkin) { 102 } else if (logekin>=fLogMinEkin) { 104 remRfaction = (logekin - fLogMinEkin) * 103 remRfaction = (logekin - fLogMinEkin) * fInvLogDelEkin; 105 ekinIndxLow = (G4int)remRfaction; 104 ekinIndxLow = (G4int)remRfaction; 106 remRfaction -= ekinIndxLow; 105 remRfaction -= ekinIndxLow; 107 } // the defaults otherwise i.e. use the low 106 } // the defaults otherwise i.e. use the lowest energy values when ekin is smaller than the minum ekin 108 // 107 // 109 DataPerEkin *perEkinLow = fMCDataPerMateria 108 DataPerEkin *perEkinLow = fMCDataPerMaterial[matindx]->fDataPerEkin[ekinIndxLow]; 110 mcToScr = perEkinLow->fMCScreening; 109 mcToScr = perEkinLow->fMCScreening; 111 mcToQ1 = perEkinLow->fMCFirstMoment; 110 mcToQ1 = perEkinLow->fMCFirstMoment; 112 mcToG2PerG1 = perEkinLow->fMCSecondMoment; 111 mcToG2PerG1 = perEkinLow->fMCSecondMoment; 113 if (remRfaction>0.) { 112 if (remRfaction>0.) { 114 DataPerEkin *perEkinHigh = fMCDataPerMater 113 DataPerEkin *perEkinHigh = fMCDataPerMaterial[matindx]->fDataPerEkin[ekinIndxLow+1]; 115 mcToScr += remRfaction*(perEkinHigh-> 114 mcToScr += remRfaction*(perEkinHigh->fMCScreening - perEkinLow->fMCScreening); 116 mcToQ1 += remRfaction*(perEkinHigh-> 115 mcToQ1 += remRfaction*(perEkinHigh->fMCFirstMoment - perEkinLow->fMCFirstMoment); 117 mcToG2PerG1 += remRfaction*(perEkinHigh-> 116 mcToG2PerG1 += remRfaction*(perEkinHigh->fMCSecondMoment - perEkinLow->fMCSecondMoment); 118 } 117 } 119 } 118 } 120 119 121 120 122 // accept cost if rndm [0,1] < return value 121 // accept cost if rndm [0,1] < return value 123 double G4GSMottCorrection::GetMottRejectionVal 122 double G4GSMottCorrection::GetMottRejectionValue(G4double logekin, G4double beta2, G4double q1, G4double cost, 124 123 G4int matindx, G4int &ekindx, G4int &deltindx) { 125 G4double val = 1.0; 124 G4double val = 1.0; 126 G4double delta = q1/(0.5+q1); 125 G4double delta = q1/(0.5+q1); 127 // check if converged to 1 for all angles => 126 // check if converged to 1 for all angles => accept cost 128 if (delta>=gMaxDelta) { 127 if (delta>=gMaxDelta) { 129 return val; 128 return val; 130 } 129 } 131 // 130 // 132 // check if kinetic energy index needs to be 131 // check if kinetic energy index needs to be determined 133 if (ekindx<0) { 132 if (ekindx<0) { 134 G4int ekinIndxLow = 0; 133 G4int ekinIndxLow = 0; 135 G4double probIndxHigh = 0.; // will be th 134 G4double probIndxHigh = 0.; // will be the prob. of taking the ekinIndxLow+1 bin 136 if (beta2>gMaxBeta2) { 135 if (beta2>gMaxBeta2) { 137 ekinIndxLow = gNumEkin - 1; 136 ekinIndxLow = gNumEkin - 1; 138 // probIndxHigh = -1. 137 // probIndxHigh = -1. 139 } else if (beta2>=fMinBeta2) { // linea 138 } else if (beta2>=fMinBeta2) { // linear interpolation on \beta^2 140 probIndxHigh = (beta2 - fMinBeta2) * fI 139 probIndxHigh = (beta2 - fMinBeta2) * fInvDelBeta2; 141 ekinIndxLow = (G4int)probIndxHigh; 140 ekinIndxLow = (G4int)probIndxHigh; 142 probIndxHigh -= ekinIndxLow; 141 probIndxHigh -= ekinIndxLow; 143 ekinIndxLow += (gNumEkin - gNumBeta2); 142 ekinIndxLow += (gNumEkin - gNumBeta2); 144 } else if (logekin>fLogMinEkin) { // linea 143 } else if (logekin>fLogMinEkin) { // linear interpolation on \ln(E_{kin}) 145 probIndxHigh = (logekin - fLogMinEkin) 144 probIndxHigh = (logekin - fLogMinEkin) * fInvLogDelEkin; 146 ekinIndxLow = (G4int)probIndxHigh; 145 ekinIndxLow = (G4int)probIndxHigh; 147 probIndxHigh -= ekinIndxLow; 146 probIndxHigh -= ekinIndxLow; 148 } // the defaults otherwise i.e. use the l 147 } // the defaults otherwise i.e. use the lowest energy values when ekin is smaller than the minum ekin 149 // 148 // 150 // check if need to take the higher ekin i 149 // check if need to take the higher ekin index 151 if (G4UniformRand()<probIndxHigh) { 150 if (G4UniformRand()<probIndxHigh) { 152 ++ekinIndxLow; 151 ++ekinIndxLow; 153 } 152 } 154 // set kinetic energy grid index 153 // set kinetic energy grid index 155 ekindx = ekinIndxLow; 154 ekindx = ekinIndxLow; 156 } 155 } 157 // check if delta value index needs to be de 156 // 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: 157 // 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) { 158 if (deltindx<0) { 160 // note: delta is for sure < gMaxDelta at 159 // note: delta is for sure < gMaxDelta at this point ( and minimum delta value is 0) 161 G4double probIndxHigh = delta*fInvDelDelta 160 G4double probIndxHigh = delta*fInvDelDelta; // will be the prob. of taking the deltIndxLow+1 bin 162 G4int deltIndxLow = (G4int)probIndxHig 161 G4int deltIndxLow = (G4int)probIndxHigh; 163 probIndxHigh -= deltIndxLow; 162 probIndxHigh -= deltIndxLow; 164 // check if need to take the higher delta 163 // check if need to take the higher delta index 165 if (G4UniformRand()<probIndxHigh) { 164 if (G4UniformRand()<probIndxHigh) { 166 ++deltIndxLow; 165 ++deltIndxLow; 167 } 166 } 168 // set the delta value grid index 167 // set the delta value grid index 169 deltindx = deltIndxLow; 168 deltindx = deltIndxLow; 170 } 169 } 171 // 170 // 172 // get the corresponding distribution 171 // get the corresponding distribution 173 DataPerDelta *perDelta = fMCDataPerMaterial 172 DataPerDelta *perDelta = fMCDataPerMaterial[matindx]->fDataPerEkin[ekindx]->fDataPerDelta[deltindx]; 174 // 173 // 175 // determine lower index of the angular bin 174 // determine lower index of the angular bin 176 G4double ang = std::sqrt(0.5*(1.-cos 175 G4double ang = std::sqrt(0.5*(1.-cost)); // sin(0.5\theta) in [0,1] 177 G4double remRfaction = ang*fInvDelAngle; 176 G4double remRfaction = ang*fInvDelAngle; 178 G4int angIndx = (G4int)remRfaction; 177 G4int angIndx = (G4int)remRfaction; 179 remRfaction -= angIndx; 178 remRfaction -= angIndx; 180 if (angIndx<gNumAngle-2) { // normal case: l 179 if (angIndx<gNumAngle-2) { // normal case: linear interpolation 181 val = remRfaction*(perDelta->fRej 180 val = remRfaction*(perDelta->fRejFuntion[angIndx+1]-perDelta->fRejFuntion[angIndx]) + perDelta->fRejFuntion[angIndx]; 182 } else { // last bin 181 } else { // last bin 183 G4double dum = ang-1.+1./fInvDelAngle; 182 G4double dum = ang-1.+1./fInvDelAngle; 184 val = perDelta->fSA + dum*(perDel 183 val = perDelta->fSA + dum*(perDelta->fSB + dum*(perDelta->fSC + dum*perDelta->fSD)); 185 } 184 } 186 return val; 185 return val; 187 } 186 } 188 187 189 188 190 void G4GSMottCorrection::Initialise() { 189 void G4GSMottCorrection::Initialise() { 191 // load Mott-correction data for each elemen 190 // load Mott-correction data for each elements that belongs to materials that are used in the detector 192 InitMCDataPerElement(); 191 InitMCDataPerElement(); 193 // clrea Mott-correction data per material 192 // clrea Mott-correction data per material 194 ClearMCDataPerMaterial(); 193 ClearMCDataPerMaterial(); 195 // initialise Mott-correction data for the m 194 // initialise Mott-correction data for the materials that are used in the detector 196 InitMCDataPerMaterials(); 195 InitMCDataPerMaterials(); 197 } 196 } 198 197 199 198 200 void G4GSMottCorrection::InitMCDataPerElement( 199 void G4GSMottCorrection::InitMCDataPerElement() { 201 // do it only once 200 // do it only once 202 if (fMCDataPerElement.size()<gMaxZet+1) { 201 if (fMCDataPerElement.size()<gMaxZet+1) { 203 fMCDataPerElement.resize(gMaxZet+1,nullptr 202 fMCDataPerElement.resize(gMaxZet+1,nullptr); 204 } 203 } 205 // loop over all materials, for those that a 204 // 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 ye 205 // corresponding data has not been loaded yet 207 G4ProductionCutsTable *thePCTable = G4Produc 206 G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable(); 208 G4int numMatCuts = (G4int)thePCTable->GetTab 207 G4int numMatCuts = (G4int)thePCTable->GetTableSize(); 209 for (G4int imc=0; imc<numMatCuts; ++imc) { 208 for (G4int imc=0; imc<numMatCuts; ++imc) { 210 const G4MaterialCutsCouple *matCut = thePC 209 const G4MaterialCutsCouple *matCut = thePCTable->GetMaterialCutsCouple(imc); 211 if (!matCut->IsUsed()) { 210 if (!matCut->IsUsed()) { 212 continue; 211 continue; 213 } 212 } 214 const G4Material *mat = matCut-> 213 const G4Material *mat = matCut->GetMaterial(); 215 const G4ElementVector *elemVect = mat->Get 214 const G4ElementVector *elemVect = mat->GetElementVector(); 216 // 215 // 217 std::size_t numElems = elemVect->size(); 216 std::size_t numElems = elemVect->size(); 218 for (std::size_t ielem=0; ielem<numElems; 217 for (std::size_t ielem=0; ielem<numElems; ++ielem) { 219 const G4Element *elem = (*elemVect)[iele 218 const G4Element *elem = (*elemVect)[ielem]; 220 G4int izet = G4lrint(elem->GetZ()); 219 G4int izet = G4lrint(elem->GetZ()); 221 if (izet>gMaxZet) { 220 if (izet>gMaxZet) { 222 izet = gMaxZet; 221 izet = gMaxZet; 223 } 222 } 224 if (!fMCDataPerElement[izet]) { 223 if (!fMCDataPerElement[izet]) { 225 LoadMCDataElement(elem); 224 LoadMCDataElement(elem); 226 } 225 } 227 } 226 } 228 } 227 } 229 } 228 } 230 229 231 230 232 void G4GSMottCorrection::InitMCDataPerMaterial 231 void G4GSMottCorrection::InitMCDataPerMaterials() { 233 // prepare size of the container 232 // prepare size of the container 234 std::size_t numMaterials = G4Material::GetNu 233 std::size_t numMaterials = G4Material::GetNumberOfMaterials(); 235 if (fMCDataPerMaterial.size()!=numMaterials) 234 if (fMCDataPerMaterial.size()!=numMaterials) { 236 fMCDataPerMaterial.resize(numMaterials); 235 fMCDataPerMaterial.resize(numMaterials); 237 } 236 } 238 // init. Mott-correction data for the Materi 237 // init. Mott-correction data for the Materials that are used in the geometry 239 G4ProductionCutsTable *thePCTable = G4Produc 238 G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable(); 240 G4int numMatCuts = (G4int)thePCTable->GetTab 239 G4int numMatCuts = (G4int)thePCTable->GetTableSize(); 241 for (G4int imc=0; imc<numMatCuts; ++imc) { 240 for (G4int imc=0; imc<numMatCuts; ++imc) { 242 const G4MaterialCutsCouple *matCut = thePC 241 const G4MaterialCutsCouple *matCut = thePCTable->GetMaterialCutsCouple(imc); 243 if (!matCut->IsUsed()) { 242 if (!matCut->IsUsed()) { 244 continue; 243 continue; 245 } 244 } 246 const G4Material *mat = matCut->GetMateria 245 const G4Material *mat = matCut->GetMaterial(); 247 if (!fMCDataPerMaterial[mat->GetIndex()]) 246 if (!fMCDataPerMaterial[mat->GetIndex()]) { 248 InitMCDataMaterial(mat); 247 InitMCDataMaterial(mat); 249 } 248 } 250 } 249 } 251 } 250 } 252 251 253 252 254 // it's called only if data has not been loade 253 // it's called only if data has not been loaded for this element yet 255 void G4GSMottCorrection::LoadMCDataElement(con 254 void G4GSMottCorrection::LoadMCDataElement(const G4Element *elem) { 256 // allocate memory 255 // allocate memory 257 G4int izet = elem->GetZasInt(); 256 G4int izet = elem->GetZasInt(); 258 if (izet>gMaxZet) { 257 if (izet>gMaxZet) { 259 izet = gMaxZet; 258 izet = gMaxZet; 260 } 259 } 261 auto perElem = new DataPerMaterial(); 260 auto perElem = new DataPerMaterial(); 262 AllocateDataPerMaterial(perElem); 261 AllocateDataPerMaterial(perElem); 263 fMCDataPerElement[izet] = perElem; 262 fMCDataPerElement[izet] = perElem; 264 // 263 // 265 // load data from file 264 // load data from file 266 std::string path = G4EmParameters::Instance( << 265 const char* tmppath = G4FindDataDir("G4LEDATA"); >> 266 if (!tmppath) { >> 267 G4Exception("G4GSMottCorrection::LoadMCDataElement()","em0006", >> 268 FatalException, >> 269 "Environment variable G4LEDATA not defined"); >> 270 return; >> 271 } >> 272 std::string path(tmppath); 267 if (fIsElectron) { 273 if (fIsElectron) { 268 path += "/msc_GS/MottCor/el/"; 274 path += "/msc_GS/MottCor/el/"; 269 } else { 275 } else { 270 path += "/msc_GS/MottCor/pos/"; 276 path += "/msc_GS/MottCor/pos/"; 271 } 277 } 272 std::string fname = path+"rej_"+gElemSymbols 278 std::string fname = path+"rej_"+gElemSymbols[izet-1]; 273 std::istringstream infile(std::ios::in); 279 std::istringstream infile(std::ios::in); 274 ReadCompressedFile(fname, infile); 280 ReadCompressedFile(fname, infile); 275 // check if file is open !!! 281 // check if file is open !!! 276 for (G4int iek=0; iek<gNumEkin; ++iek) { 282 for (G4int iek=0; iek<gNumEkin; ++iek) { 277 DataPerEkin *perEkin = perElem->fDataPerEk 283 DataPerEkin *perEkin = perElem->fDataPerEkin[iek]; 278 // 1. get the 3 Mott-correction factors fo 284 // 1. get the 3 Mott-correction factors for the current kinetic energy 279 infile >> perEkin->fMCScreening; 285 infile >> perEkin->fMCScreening; 280 infile >> perEkin->fMCFirstMoment; 286 infile >> perEkin->fMCFirstMoment; 281 infile >> perEkin->fMCSecondMoment; 287 infile >> perEkin->fMCSecondMoment; 282 // 2. load each data per delta: 288 // 2. load each data per delta: 283 for (G4int idel=0; idel<gNumDelta; ++idel) 289 for (G4int idel=0; idel<gNumDelta; ++idel) { 284 DataPerDelta *perDelta = perEkin->fDataP 290 DataPerDelta *perDelta = perEkin->fDataPerDelta[idel]; 285 // 2./a. : first the rejection function 291 // 2./a. : first the rejection function values 286 for (G4int iang=0; iang<gNumAngle; ++ian 292 for (G4int iang=0; iang<gNumAngle; ++iang) { 287 infile >> perDelta->fRejFuntion[iang]; 293 infile >> perDelta->fRejFuntion[iang]; 288 } 294 } 289 // 2./b. : then the 4 spline parameter f 295 // 2./b. : then the 4 spline parameter for the last bin 290 infile >> perDelta->fSA; 296 infile >> perDelta->fSA; 291 infile >> perDelta->fSB; 297 infile >> perDelta->fSB; 292 infile >> perDelta->fSC; 298 infile >> perDelta->fSC; 293 infile >> perDelta->fSD; 299 infile >> perDelta->fSD; 294 } 300 } 295 } 301 } 296 } 302 } 297 303 298 // uncompress one data file into the input str 304 // uncompress one data file into the input string stream 299 void G4GSMottCorrection::ReadCompressedFile(co << 305 void G4GSMottCorrection::ReadCompressedFile(std::string fname, std::istringstream &iss) { 300 std::string *dataString = nullptr; 306 std::string *dataString = nullptr; 301 std::string compfilename(fname+".z"); 307 std::string compfilename(fname+".z"); 302 // create input stream with binary mode oper 308 // create input stream with binary mode operation and positioning at the end of the file 303 std::ifstream in(compfilename, std::ios::bin 309 std::ifstream in(compfilename, std::ios::binary | std::ios::ate); 304 if (in.good()) { 310 if (in.good()) { 305 // get current position in the stream (wa 311 // get current position in the stream (was set to the end) 306 std::streamoff fileSize = in.tellg(); 312 std::streamoff fileSize = in.tellg(); 307 // set current position being the beginni 313 // set current position being the beginning of the stream 308 in.seekg(0,std::ios::beg); 314 in.seekg(0,std::ios::beg); 309 // create (zlib) byte buffer for the data 315 // create (zlib) byte buffer for the data 310 Bytef *compdata = new Bytef[fileSize]; 316 Bytef *compdata = new Bytef[fileSize]; 311 while(in) { 317 while(in) { 312 in.read((char*)compdata, fileSize); 318 in.read((char*)compdata, fileSize); 313 } 319 } 314 // create (zlib) byte buffer for the unco 320 // create (zlib) byte buffer for the uncompressed data 315 uLongf complen = (uLongf)(fileSize*4); 321 uLongf complen = (uLongf)(fileSize*4); 316 Bytef *uncompdata = new Bytef[complen]; 322 Bytef *uncompdata = new Bytef[complen]; 317 while (Z_OK!=uncompress(uncompdata, &comp 323 while (Z_OK!=uncompress(uncompdata, &complen, compdata, fileSize)) { 318 // increase uncompressed byte buffer 324 // increase uncompressed byte buffer 319 delete[] uncompdata; 325 delete[] uncompdata; 320 complen *= 2; 326 complen *= 2; 321 uncompdata = new Bytef[complen]; 327 uncompdata = new Bytef[complen]; 322 } 328 } 323 // delete the compressed data buffer 329 // delete the compressed data buffer 324 delete [] compdata; 330 delete [] compdata; 325 // create a string from the uncompressed 331 // create a string from the uncompressed data (will be deallocated by the caller) 326 dataString = new std::string((char*)uncom 332 dataString = new std::string((char*)uncompdata, (long)complen); 327 // delete the uncompressed data buffer 333 // delete the uncompressed data buffer 328 delete [] uncompdata; 334 delete [] uncompdata; 329 } else { 335 } else { 330 std::string msg = " Problem while trying 336 std::string msg = " Problem while trying to read " + compfilename + " data file.\n"; 331 G4Exception("G4GSMottCorrection::ReadCompr 337 G4Exception("G4GSMottCorrection::ReadCompressedFile","em0006", FatalException,msg.c_str()); 332 return; 338 return; 333 } 339 } 334 // create the input string stream from the d 340 // create the input string stream from the data string 335 if (dataString) { 341 if (dataString) { 336 iss.str(*dataString); 342 iss.str(*dataString); 337 in.close(); 343 in.close(); 338 delete dataString; 344 delete dataString; 339 } 345 } 340 } 346 } 341 347 342 348 343 void G4GSMottCorrection::InitMCDataMaterial(co 349 void G4GSMottCorrection::InitMCDataMaterial(const G4Material *mat) { 344 constexpr G4double const1 = 7821.6; / 350 constexpr G4double const1 = 7821.6; // [cm2/g] 345 constexpr G4double const2 = 0.1569; / 351 constexpr G4double const2 = 0.1569; // [cm2 MeV2 / g] 346 constexpr G4double finstrc2 = 5.325135453E-5 352 constexpr G4double finstrc2 = 5.325135453E-5; // fine-structure const. square 347 353 348 G4double constFactor = CLHEP::electro 354 G4double constFactor = CLHEP::electron_mass_c2*CLHEP::fine_structure_const/0.88534; 349 constFactor *= constFactor; 355 constFactor *= constFactor; // (mc^2)^2\alpha^2/( C_{TF}^2) 350 // allocate memory 356 // allocate memory 351 auto perMat = new DataPerMaterial(); 357 auto perMat = new DataPerMaterial(); 352 AllocateDataPerMaterial(perMat); 358 AllocateDataPerMaterial(perMat); 353 fMCDataPerMaterial[mat->GetIndex()] = perMat 359 fMCDataPerMaterial[mat->GetIndex()] = perMat; 354 // 360 // 355 const G4ElementVector* elemVect = 361 const G4ElementVector* elemVect = mat->GetElementVector(); 356 const G4int numElems = 362 const G4int numElems = (G4int)mat->GetNumberOfElements(); 357 const G4double* nbAtomsPerVolVect = 363 const G4double* nbAtomsPerVolVect = mat->GetVecNbOfAtomsPerVolume(); 358 G4double totNbAtomsPerVol = 364 G4double totNbAtomsPerVol = mat->GetTotNbOfAtomsPerVolume(); 359 // 365 // 360 // 1. Compute material dependent part of Mol 366 // 1. Compute material dependent part of Moliere's b_c \chi_c^2 361 // (with \xi=1 (i.e. total sub-threshold 367 // (with \xi=1 (i.e. total sub-threshold scattering power correction) 362 G4double moliereBc = 0.0; 368 G4double moliereBc = 0.0; 363 G4double moliereXc2 = 0.0; 369 G4double moliereXc2 = 0.0; 364 G4double zs = 0.0; 370 G4double zs = 0.0; 365 G4double ze = 0.0; 371 G4double ze = 0.0; 366 G4double zx = 0.0; 372 G4double zx = 0.0; 367 G4double sa = 0.0; 373 G4double sa = 0.0; 368 G4double xi = 1.0; 374 G4double xi = 1.0; 369 for (G4int ielem=0; ielem<numElems; ++ielem) 375 for (G4int ielem=0; ielem<numElems; ++ielem) { 370 G4double zet = (*elemVect)[ielem]->GetZ(); 376 G4double zet = (*elemVect)[ielem]->GetZ(); 371 if (zet>gMaxZet) { 377 if (zet>gMaxZet) { 372 zet = (G4double)gMaxZet; 378 zet = (G4double)gMaxZet; 373 } 379 } 374 G4double iwa = (*elemVect)[ielem]->GetN() 380 G4double iwa = (*elemVect)[ielem]->GetN(); 375 G4double ipz = nbAtomsPerVolVect[ielem]/t 381 G4double ipz = nbAtomsPerVolVect[ielem]/totNbAtomsPerVol; 376 G4double dum = ipz*zet*(zet+xi); 382 G4double dum = ipz*zet*(zet+xi); 377 zs += dum; 383 zs += dum; 378 ze += dum*(-2.0/3.0)*G4Log(zet); 384 ze += dum*(-2.0/3.0)*G4Log(zet); 379 zx += dum*G4Log(1.0+3.34*finstrc 385 zx += dum*G4Log(1.0+3.34*finstrc2*zet*zet); 380 sa += ipz*iwa; 386 sa += ipz*iwa; 381 } 387 } 382 G4double density = mat->GetDensity()*CLHEP:: 388 G4double density = mat->GetDensity()*CLHEP::cm3/CLHEP::g; // [g/cm3] 383 // 389 // 384 G4double z0 = (0.0 == sa) ? 0.0 : zs/sa; << 390 moliereBc = const1*density*zs/sa*G4Exp(ze/zs)/G4Exp(zx/zs); //[1/cm] 385 G4double z1 = (0.0 == zs) ? 0.0 : (ze - zx)/ << 391 moliereXc2 = const2*density*zs/sa; // [MeV2/cm] 386 << 387 moliereBc = const1*density*z0*G4Exp(z1); / << 388 moliereXc2 = const2*density*z0; // [MeV2/cm << 389 // change to Geant4 internal units of 1/leng 392 // change to Geant4 internal units of 1/length and energ2/length 390 moliereBc *= 1.0/CLHEP::cm; 393 moliereBc *= 1.0/CLHEP::cm; 391 moliereXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::c 394 moliereXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm; 392 // 395 // 393 // 2. loop over the kinetic energy grid 396 // 2. loop over the kinetic energy grid 394 for (G4int iek=0; iek<gNumEkin; ++iek) { 397 for (G4int iek=0; iek<gNumEkin; ++iek) { 395 // 2./a. set current kinetic energy and pt 398 // 2./a. set current kinetic energy and pt2 value 396 G4double ekin = G4Exp(fLogMinEk 399 G4double ekin = G4Exp(fLogMinEkin+iek/fInvLogDelEkin); 397 G4double pt2 = ekin*(ekin+2.0*CLHEP::el 400 G4double pt2 = ekin*(ekin+2.0*CLHEP::electron_mass_c2); 398 if (ekin>gMidEkin) { 401 if (ekin>gMidEkin) { 399 G4double b2 = fMinBeta2+(iek-(gNumEk 402 G4double b2 = fMinBeta2+(iek-(gNumEkin-gNumBeta2))/fInvDelBeta2; 400 ekin = CLHEP::electron_mass_c2*(1./std 403 ekin = CLHEP::electron_mass_c2*(1./std::sqrt(1.-b2)-1.); 401 pt2 = ekin*(ekin+2.0*CLHEP::electron_ 404 pt2 = ekin*(ekin+2.0*CLHEP::electron_mass_c2); 402 } 405 } 403 // 2./b. loop over the elements at the cur 406 // 2./b. loop over the elements at the current kinetic energy point 404 for (G4int ielem=0; ielem<numElems; ++iele 407 for (G4int ielem=0; ielem<numElems; ++ielem) { 405 const G4Element *elem = (*elemVect)[iele 408 const G4Element *elem = (*elemVect)[ielem]; 406 G4double zet = elem->GetZ(); 409 G4double zet = elem->GetZ(); 407 if (zet>gMaxZet) { 410 if (zet>gMaxZet) { 408 zet = (G4double)gMaxZet; 411 zet = (G4double)gMaxZet; 409 } 412 } 410 G4int izet = G4lrint(zet); 413 G4int izet = G4lrint(zet); 411 // xi should be one i.e. z(z+1) since to 414 // xi should be one i.e. z(z+1) since total sub-threshold scattering power correction 412 G4double nZZPlus1 = nbAtomsPerVolVect[i 415 G4double nZZPlus1 = nbAtomsPerVolVect[ielem]*zet*(zet+1.0)/totNbAtomsPerVol; 413 G4double Z23 = std::pow(zet,2./3.) 416 G4double Z23 = std::pow(zet,2./3.); 414 // 417 // 415 DataPerEkin *perElemPerEkin = fMCDataPe 418 DataPerEkin *perElemPerEkin = fMCDataPerElement[izet]->fDataPerEkin[iek]; 416 DataPerEkin *perMatPerEkin = perMat->f 419 DataPerEkin *perMatPerEkin = perMat->fDataPerEkin[iek]; 417 // 420 // 418 // 2./b./(i) Add the 3 Mott-correction f 421 // 2./b./(i) Add the 3 Mott-correction factors 419 G4double mcScrCF = perElemPerEkin->fMCSc 422 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 corre 423 // compute the screening parameter correction factor (Z_i contribution to the material) 421 // src_{mc} = C \exp\left[ \frac{ \sum_i 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)} 422 // with C = \frac{(mc^2)^\alpha^2} {4(pc 425 // 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_ 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 424 perMatPerEkin->fMCScreening += nZZPlus1* 427 perMatPerEkin->fMCScreening += nZZPlus1*G4Log(Z23*mcScrCF); 425 // compute the corrected screening param 428 // compute the corrected screening parameter for the current Z_i and E_{kin} 426 // src(Z_i)_{mc} = \frac{(mc^2)^\alpha^2 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] 427 mcScrCF *= constFactor*Z23/(4.*pt2); 430 mcScrCF *= constFactor*Z23/(4.*pt2); 428 // compute first moment correction facto 431 // compute first moment correction factor 429 // q1_{mc} = \frac{ \sum_i n_i Z_i(Z_i+1 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} 430 // where: 433 // where: 431 // A_i(src(Z_i)_{mc}) = [\ln(1+1/src(Z_i 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 432 // B_i = \beta_i \gamma_i with beta_i(Z_ 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) 433 // and \gamma_i = \sigma(Z_i)_{el}^(MC-D 436 // 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/ 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 435 // A_i x B_i is stored in file per e-/e+ 438 // 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_ 439 // here we compute the \sum_i n_i Z_i(Z_i+1) A_i B_i part 437 perMatPerEkin->fMCFirstMoment += nZZPlus 440 perMatPerEkin->fMCFirstMoment += nZZPlus1*(G4Log(1.+1./mcScrCF)-1./(1.+mcScrCF))*perElemPerEkin->fMCFirstMoment; 438 // compute the second moment correction 441 // compute the second moment correction factor 439 // [G2/G1]_{mc} = \frac{ \sum_i n_i Z_i( 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} 440 // with A_i(Z_i) = G2(Z_i)^{PWA}/G1(Z_i) 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}} 441 // here we compute the \sum_i n_i Z_i(Z_ 444 // here we compute the \sum_i n_i Z_i(Z_i+1) A_i part 442 perMatPerEkin->fMCSecondMoment += nZZPlu 445 perMatPerEkin->fMCSecondMoment += nZZPlus1*perElemPerEkin->fMCSecondMoment; 443 // 446 // 444 // 2./b./(ii) Go for the rejection funti 447 // 2./b./(ii) Go for the rejection funtion part 445 // I. loop over delta values 448 // I. loop over delta values 446 for (G4int idel=0; idel<gNumDelta; ++ide 449 for (G4int idel=0; idel<gNumDelta; ++idel) { 447 DataPerDelta *perMatPerDelta = perMat 450 DataPerDelta *perMatPerDelta = perMatPerEkin->fDataPerDelta[idel]; 448 DataPerDelta *perElemPerDelta = perEle 451 DataPerDelta *perElemPerDelta = perElemPerEkin->fDataPerDelta[idel]; 449 // I./a. loop over angles (i.e. the \s 452 // 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; ++i 453 for (G4int iang=0; iang<gNumAngle; ++iang) { 451 perMatPerDelta->fRejFuntion[iang] += 454 perMatPerDelta->fRejFuntion[iang] += nZZPlus1*perElemPerDelta->fRejFuntion[iang]; 452 } 455 } 453 // I./b. get the last bin spline param 456 // I./b. get the last bin spline parameters and add them (a+bx+cx^2+dx^3) 454 perMatPerDelta->fSA += nZZPlus1*perEle 457 perMatPerDelta->fSA += nZZPlus1*perElemPerDelta->fSA; 455 perMatPerDelta->fSB += nZZPlus1*perEle 458 perMatPerDelta->fSB += nZZPlus1*perElemPerDelta->fSB; 456 perMatPerDelta->fSC += nZZPlus1*perEle 459 perMatPerDelta->fSC += nZZPlus1*perElemPerDelta->fSC; 457 perMatPerDelta->fSD += nZZPlus1*perEle 460 perMatPerDelta->fSD += nZZPlus1*perElemPerDelta->fSD; 458 } 461 } 459 // 462 // 460 // 2./b./(iii) When the last element has 463 // 2./b./(iii) When the last element has been added: 461 if (ielem==numElems-1) { 464 if (ielem==numElems-1) { 462 // 465 // 463 // 1. the remaining part of the sreeni 466 // 1. the remaining part of the sreening correction and divide the corrected screening par. with Moliere's one: 464 // (Moliere screening parameter = m 467 // (Moliere screening parameter = moliereXc2/(4(pc)^2 moliereBc) ) 465 G4double dumScr = G4Exp(perMatPerEki 468 G4double dumScr = G4Exp(perMatPerEkin->fMCScreening/zs); 466 perMatPerEkin->fMCScreening = constFac 469 perMatPerEkin->fMCScreening = constFactor*dumScr*moliereBc/moliereXc2; 467 // 470 // 468 // 2. the remaining part of the first 471 // 2. the remaining part of the first moment correction and divide by the one computed by using the corrected 469 // screening parameter (= (mc^2)^\a 472 // screening parameter (= (mc^2)^\alpha^2/(4(pc)^2C_{TF}^2) dumScr 470 G4double scrCorTed = constFactor*dumSc 473 G4double scrCorTed = constFactor*dumScr/(4.*pt2); 471 G4double dum0 = G4Log(1.+1./scrCo 474 G4double dum0 = G4Log(1.+1./scrCorTed); 472 perMatPerEkin->fMCFirstMoment = perMat 475 perMatPerEkin->fMCFirstMoment = perMatPerEkin->fMCFirstMoment/(zs*(dum0-1./(1.+scrCorTed))); 473 // 476 // 474 // 3. the remaining part of the second 477 // 3. the remaining part of the second moment correction and divide by the one computed by using the corrected 475 // screening parameter 478 // screening parameter 476 G4double G2PerG1 = 3.*(1.+scrCorTed 479 G4double G2PerG1 = 3.*(1.+scrCorTed)*((1.+2.*scrCorTed)*dum0-2.)/((1.+scrCorTed)*dum0-1.); 477 perMatPerEkin->fMCSecondMoment = perMa 480 perMatPerEkin->fMCSecondMoment = perMatPerEkin->fMCSecondMoment/(zs*G2PerG1); 478 // 481 // 479 // 4. scale the maximum of the rejecti 482 // 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 483 // I. loop over delta values 481 for (G4int idel=0; idel<gNumDelta; ++i 484 for (G4int idel=0; idel<gNumDelta; ++idel) { 482 DataPerDelta *perMatPerDelta = perM 485 DataPerDelta *perMatPerDelta = perMatPerEkin->fDataPerDelta[idel]; 483 G4double maxVal = -1.; 486 G4double maxVal = -1.; 484 // II. llop over angles 487 // II. llop over angles 485 for (G4int iang=0; iang<gNumAngle; + 488 for (G4int iang=0; iang<gNumAngle; ++iang) { 486 if (perMatPerDelta->fRejFuntion[ia 489 if (perMatPerDelta->fRejFuntion[iang]>maxVal) 487 maxVal = perMatPerDelta->fRejFun 490 maxVal = perMatPerDelta->fRejFuntion[iang]; 488 } 491 } 489 for (G4int iang=0; iang<gNumAngle; + 492 for (G4int iang=0; iang<gNumAngle; ++iang) { 490 perMatPerDelta->fRejFuntion[iang] 493 perMatPerDelta->fRejFuntion[iang] /=maxVal; 491 } 494 } 492 perMatPerDelta->fSA /= maxVal; 495 perMatPerDelta->fSA /= maxVal; 493 perMatPerDelta->fSB /= maxVal; 496 perMatPerDelta->fSB /= maxVal; 494 perMatPerDelta->fSC /= maxVal; 497 perMatPerDelta->fSC /= maxVal; 495 perMatPerDelta->fSD /= maxVal; 498 perMatPerDelta->fSD /= maxVal; 496 } 499 } 497 } 500 } 498 } 501 } 499 } 502 } 500 } 503 } 501 504 502 505 503 void G4GSMottCorrection::AllocateDataPerMateri 506 void G4GSMottCorrection::AllocateDataPerMaterial(DataPerMaterial *data) { 504 data->fDataPerEkin = new DataPerEkin*[gNumEk 507 data->fDataPerEkin = new DataPerEkin*[gNumEkin](); 505 for (G4int iek=0; iek<gNumEkin; ++iek) { 508 for (G4int iek=0; iek<gNumEkin; ++iek) { 506 auto perEkin = new DataPerEkin(); 509 auto perEkin = new DataPerEkin(); 507 perEkin->fDataPerDelta = new DataPerDelta* 510 perEkin->fDataPerDelta = new DataPerDelta*[gNumDelta](); 508 for (G4int idel=0; idel<gNumDelta; ++idel) 511 for (G4int idel=0; idel<gNumDelta; ++idel) { 509 auto perDelta = new DataP 512 auto perDelta = new DataPerDelta(); 510 perDelta->fRejFuntion = new doubl 513 perDelta->fRejFuntion = new double[gNumAngle](); 511 perEkin->fDataPerDelta[idel] = perDelta; 514 perEkin->fDataPerDelta[idel] = perDelta; 512 } 515 } 513 data->fDataPerEkin[iek] = perEkin; 516 data->fDataPerEkin[iek] = perEkin; 514 } 517 } 515 } 518 } 516 519 517 void G4GSMottCorrection::DeAllocateDataPerMate 520 void G4GSMottCorrection::DeAllocateDataPerMaterial(DataPerMaterial *data) { 518 for (G4int iek=0; iek<gNumEkin; ++iek) { 521 for (G4int iek=0; iek<gNumEkin; ++iek) { 519 DataPerEkin *perEkin = data->fDataPerEkin[ 522 DataPerEkin *perEkin = data->fDataPerEkin[iek]; //new DataPerEkin(); 520 for (G4int idel=0; idel<gNumDelta; ++idel) 523 for (G4int idel=0; idel<gNumDelta; ++idel) { 521 DataPerDelta *perDelta = perEkin->fDataP 524 DataPerDelta *perDelta = perEkin->fDataPerDelta[idel]; 522 delete [] perDelta->fRejFuntion; 525 delete [] perDelta->fRejFuntion; 523 delete perDelta; 526 delete perDelta; 524 } 527 } 525 delete [] perEkin->fDataPerDelta; 528 delete [] perEkin->fDataPerDelta; 526 delete perEkin; 529 delete perEkin; 527 } 530 } 528 delete [] data->fDataPerEkin; 531 delete [] data->fDataPerEkin; 529 } 532 } 530 533 531 534 532 void G4GSMottCorrection::ClearMCDataPerElement 535 void G4GSMottCorrection::ClearMCDataPerElement() { 533 for (std::size_t i=0; i<fMCDataPerElement.si 536 for (std::size_t i=0; i<fMCDataPerElement.size(); ++i) { 534 if (fMCDataPerElement[i]) { 537 if (fMCDataPerElement[i]) { 535 DeAllocateDataPerMaterial(fMCDataPerElem 538 DeAllocateDataPerMaterial(fMCDataPerElement[i]); 536 delete fMCDataPerElement[i]; 539 delete fMCDataPerElement[i]; 537 } 540 } 538 } 541 } 539 fMCDataPerElement.clear(); 542 fMCDataPerElement.clear(); 540 } 543 } 541 544 542 void G4GSMottCorrection::ClearMCDataPerMateria 545 void G4GSMottCorrection::ClearMCDataPerMaterial() { 543 for (std::size_t i=0; i<fMCDataPerMaterial.s 546 for (std::size_t i=0; i<fMCDataPerMaterial.size(); ++i) { 544 if (fMCDataPerMaterial[i]) { 547 if (fMCDataPerMaterial[i]) { 545 DeAllocateDataPerMaterial(fMCDataPerMate 548 DeAllocateDataPerMaterial(fMCDataPerMaterial[i]); 546 delete fMCDataPerMaterial[i]; 549 delete fMCDataPerMaterial[i]; 547 } 550 } 548 } 551 } 549 fMCDataPerMaterial.clear(); 552 fMCDataPerMaterial.clear(); 550 } 553 } 551 554