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 // GEANT4 Class implementation file 29 // GEANT4 Class implementation file 30 // 30 // 31 // File name: G4GoudsmitSaundersonMscModel 31 // File name: G4GoudsmitSaundersonMscModel 32 // 32 // 33 // Author: Mihaly Novak / (Omrane Kadri 33 // Author: Mihaly Novak / (Omrane Kadri) 34 // 34 // 35 // Creation date: 20.02.2009 35 // Creation date: 20.02.2009 36 // 36 // 37 // Modifications: 37 // Modifications: 38 // 04.03.2009 V.Ivanchenko cleanup and format 38 // 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style 39 // 12.05.2010 O.Kadri: adding Qn1 and Qn12 as 39 // 12.05.2010 O.Kadri: adding Qn1 and Qn12 as private doubles 40 // 18.05.2015 M. Novak provide PRELIMINARY ver 40 // 18.05.2015 M. Novak provide PRELIMINARY verison of the revised class 41 // This class has been revised and 41 // This class has been revised and updated, new methods added. 42 // A new version of Kawrakow-Bielaj 42 // A new version of Kawrakow-Bielajew Goudsmit-Saunderson MSC model 43 // based on the screened Rutherford 43 // based on the screened Rutherford DCS for elastic scattering of 44 // electrons/positrons has been int 44 // electrons/positrons has been introduced[1,2]. The corresponding MSC 45 // angular distributions over a 2D 45 // angular distributions over a 2D parameter grid have been recomputed 46 // and the CDFs are now stored in a 46 // and the CDFs are now stored in a variable transformed (smooth) form[2,3] 47 // together with the corresponding 47 // together with the corresponding rational interpolation parameters. 48 // These angular distributions are 48 // These angular distributions are handled by the new 49 // G4GoudsmitSaundersonTable class 49 // G4GoudsmitSaundersonTable class that is responsible to sample if 50 // it was no, single, few or multip 50 // it was no, single, few or multiple scattering case and delivers the 51 // angular deflection (i.e. cos(the 51 // angular deflection (i.e. cos(theta) and sin(theta)). 52 // Two screening options are provid 52 // Two screening options are provided: 53 // - if fIsUsePWATotalXsecData=TRU 53 // - if fIsUsePWATotalXsecData=TRUE i.e. SetOptionPWAScreening(TRUE) 54 // was called before initialisat 54 // was called before initialisation: screening parameter value A is 55 // determined such that the firs 55 // determined such that the first transport coefficient G1(A) 56 // computed according to the scr 56 // computed according to the screened Rutherford DCS for elastic 57 // scattering will reproduce the 57 // scattering will reproduce the one computed from the PWA elastic 58 // and first transport mean free 58 // and first transport mean free paths[4]. 59 // - if fIsUsePWATotalXsecData=FAL 59 // - if fIsUsePWATotalXsecData=FALSE i.e. default value or 60 // SetOptionPWAScreening(FALSE) 60 // SetOptionPWAScreening(FALSE) was called before initialisation: 61 // screening parameter value A i 61 // screening parameter value A is computed according to Moliere's 62 // formula (by using material de 62 // formula (by using material dependent parameters \chi_cc2 and b_c 63 // precomputed for each material 63 // precomputed for each material used at initialization in 64 // G4GoudsmitSaundersonTable) [3 64 // G4GoudsmitSaundersonTable) [3] 65 // Elastic and first trasport mean 65 // Elastic and first trasport mean free paths are used consistently. 66 // The new version is self-consiste 66 // The new version is self-consistent, several times faster, more 67 // robust and accurate compared to 67 // robust and accurate compared to the earlier version. 68 // Spin effects as well as a more a 68 // Spin effects as well as a more accurate energy loss correction and 69 // computations of Lewis moments wi 69 // computations of Lewis moments will be implemented later on. 70 // 02.09.2015 M. Novak: first version of new s 70 // 02.09.2015 M. Novak: first version of new step limit is provided. 71 // fUseSafetyPlus corresponds to Ur 71 // fUseSafetyPlus corresponds to Urban fUseSafety (default) 72 // fUseDistanceToBoundary correspon 72 // fUseDistanceToBoundary corresponds to Urban fUseDistanceToBoundary 73 // fUseSafety corresponds to EGSnr 73 // fUseSafety corresponds to EGSnrc error-free stepping algorithm 74 // Range factor can be significantl 74 // Range factor can be significantly higher at each case than in Urban. 75 // 23.08.2017 M. Novak: added corrections to a 75 // 23.08.2017 M. Novak: added corrections to account spin effects (Mott-correction). 76 // It can be activated by setting t 76 // It can be activated by setting the fIsMottCorrection flag to be true 77 // before initialization using the 77 // before initialization using the SetOptionMottCorrection() public method. 78 // The fMottCorrection member is re 78 // The fMottCorrection member is responsible to handle pre-computed Mott 79 // correction (rejection) functions 79 // correction (rejection) functions obtained by numerically computing 80 // Goudsmit-Saunderson agnular dist 80 // Goudsmit-Saunderson agnular distributions based on a DCS accounting spin 81 // effects and screening correction 81 // effects and screening corrections. The DCS used to compute the accurate 82 // GS angular distributions is: DCS 82 // GS angular distributions is: DCS_{cor} = DCS_{SR}x[ DCS_{R}/DCS_{Mott}] where : 83 // # DCS_{SR} is the relativisti 83 // # DCS_{SR} is the relativistic Screened-Rutherford DCS (first Born approximate 84 // solution of the Klein-Gordo 84 // solution of the Klein-Gordon i.e. relativistic Schrodinger equation => 85 // scattering of spinless e- o 85 // scattering of spinless e- on exponentially screened Coulomb potential) 86 // note: the default (without 86 // note: the default (without using Mott-correction) GS angular distributions 87 // are based on this DCS_{SR} 87 // are based on this DCS_{SR} with Moliere's screening parameter! 88 // # DCS_{R} is the Rutherford D 88 // # DCS_{R} is the Rutherford DCS which is the same as above but without 89 // screening 89 // screening 90 // # DCS_{Mott} is the Mott DCS 90 // # DCS_{Mott} is the Mott DCS i.e. solution of the Dirac equation with a bare 91 // Coulomb potential i.e. scat 91 // Coulomb potential i.e. scattering of particles with spin (e- or e+) on a 92 // point-like unscreened Coulo 92 // point-like unscreened Coulomb potential 93 // # moreover, the screening par 93 // # moreover, the screening parameter of the DCS_{cor} was determined such that 94 // the DCS_{cor} with this cor 94 // the DCS_{cor} with this corrected screening parameter reproduce the first 95 // transport cross sections ob 95 // transport cross sections obtained from the corresponding most accurate DCS 96 // (i.e. from elsepa [4]) 96 // (i.e. from elsepa [4]) 97 // Unlike the default GS, the Mott- 97 // Unlike the default GS, the Mott-corrected angular distributions are particle type 98 // (different for e- and e+ <= the 98 // (different for e- and e+ <= the DCS_{Mott} and the screening correction) and target 99 // (Z and material) dependent. 99 // (Z and material) dependent. 100 // 27.10.2017 M. Novak: 100 // 27.10.2017 M. Novak: 101 // - Mott-correction flag is set no 101 // - Mott-correction flag is set now through the G4EmParameters 102 // - new form of PWA correction to 102 // - new form of PWA correction to integrated quantities and screening (default) 103 // - changed step limit flag conven 103 // - changed step limit flag conventions: 104 // # fUseSafety corresponds to U 104 // # fUseSafety corresponds to Urban's fUseSafety 105 // # fUseDistanceToBoundary corr 105 // # fUseDistanceToBoundary corresponds to Urban's fUseDistanceToBoundary 106 // # fUseSafetyPlus corresponds 106 // # fUseSafetyPlus corresponds to the error-free stepping algorithm 107 // 02.02.2018 M. Novak: implemented CrossSecti 107 // 02.02.2018 M. Novak: implemented CrossSectionPerVolume interface method (used only for testing) 108 // 108 // 109 // Class description: 109 // Class description: 110 // Kawrakow-Bielajew Goudsmit-Saunderson MSC 110 // Kawrakow-Bielajew Goudsmit-Saunderson MSC model based on the screened Rutherford DCS 111 // for elastic scattering of e-/e+. Option, 111 // for elastic scattering of e-/e+. Option, to include (Mott) correction (see above), is 112 // also available now (SetOptionMottCorrecti 112 // also available now (SetOptionMottCorrection(true)). An EGSnrc like error-free stepping 113 // algorithm (UseSafety) is available beyond 113 // algorithm (UseSafety) is available beyond the usual Geant4 step limitation algorithms 114 // and true to geomerty and geometry to true 114 // and true to geomerty and geometry to true step length computations that were adopted 115 // from the Urban model[5]. The most accurat 115 // from the Urban model[5]. The most accurate setting: error-free stepping i.e. the 116 // UseSafetyPlus MSC step limit with Mott-co 116 // UseSafetyPlus MSC step limit with Mott-correction (SetOptionMottCorrection(true)). Both 117 // are expected to be set through the G4EmPa 117 // are expected to be set through the G4EmParameters singleton before initialisation: 118 // # G4EmParameters::Instance()->SetMscStep 118 // # G4EmParameters::Instance()->SetMscStepLimitType(fUseSafetyPlus); 119 // # G4EmParameters::Instance()->SetUseMott 119 // # G4EmParameters::Instance()->SetUseMottCorrection(true); 120 // 120 // 121 // 121 // 122 // References: 122 // References: 123 // [1] A.F.Bielajew, NIMB 111 (1996) 195-208 123 // [1] A.F.Bielajew, NIMB 111 (1996) 195-208 124 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(19 124 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(1998) 325-336 125 // [3] I.Kawrakow, E.Mainegra-Hing, D.W.O.Ro 125 // [3] I.Kawrakow, E.Mainegra-Hing, D.W.O.Rogers, F.Tessier,B.R.B.Walters, NRCC 126 // Report PIRS-701 (2013) 126 // Report PIRS-701 (2013) 127 // [4] F.Salvat, A.Jablonski, C.J. Powell, C 127 // [4] F.Salvat, A.Jablonski, C.J. Powell, CPC 165(2005) 157-190 128 // [5] L.Urban, Preprint CERN-OPEN-2006-077 128 // [5] L.Urban, Preprint CERN-OPEN-2006-077 (2006) 129 // 129 // 130 // ------------------------------------------- 130 // ----------------------------------------------------------------------------- 131 131 132 132 133 #include "G4GoudsmitSaundersonMscModel.hh" 133 #include "G4GoudsmitSaundersonMscModel.hh" 134 134 135 #include "G4GoudsmitSaundersonTable.hh" 135 #include "G4GoudsmitSaundersonTable.hh" 136 #include "G4GSPWACorrections.hh" 136 #include "G4GSPWACorrections.hh" 137 137 138 #include "G4PhysicalConstants.hh" 138 #include "G4PhysicalConstants.hh" 139 #include "G4SystemOfUnits.hh" 139 #include "G4SystemOfUnits.hh" 140 140 141 #include "G4ParticleChangeForMSC.hh" 141 #include "G4ParticleChangeForMSC.hh" 142 #include "G4DynamicParticle.hh" 142 #include "G4DynamicParticle.hh" 143 #include "G4Electron.hh" 143 #include "G4Electron.hh" 144 #include "G4Positron.hh" 144 #include "G4Positron.hh" 145 145 146 #include "G4LossTableManager.hh" 146 #include "G4LossTableManager.hh" 147 #include "G4EmParameters.hh" 147 #include "G4EmParameters.hh" 148 #include "G4Track.hh" 148 #include "G4Track.hh" 149 #include "G4PhysicsTable.hh" 149 #include "G4PhysicsTable.hh" 150 #include "Randomize.hh" 150 #include "Randomize.hh" 151 #include "G4Log.hh" 151 #include "G4Log.hh" 152 #include "G4Exp.hh" 152 #include "G4Exp.hh" 153 #include "G4Pow.hh" 153 #include "G4Pow.hh" 154 #include <fstream> 154 #include <fstream> 155 155 156 156 157 // set accurate energy loss and dispalcement s 157 // set accurate energy loss and dispalcement sampling to be always on now 158 G4bool G4GoudsmitSaundersonMscModel::gIsUseAcc 158 G4bool G4GoudsmitSaundersonMscModel::gIsUseAccurate = true; 159 // set the usual optimization to be always act 159 // set the usual optimization to be always active now 160 G4bool G4GoudsmitSaundersonMscModel::gIsOptimi 160 G4bool G4GoudsmitSaundersonMscModel::gIsOptimizationOn = true; 161 161 162 162 163 G4GoudsmitSaundersonMscModel::G4GoudsmitSaunde 163 G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel(const G4String& nam) 164 : G4VMscModel(nam) { 164 : G4VMscModel(nam) { 165 charge = 0; 165 charge = 0; 166 currentMaterialIndex = -1; 166 currentMaterialIndex = -1; 167 // 167 // >> 168 lambdalimit = 1.*mm; 168 fr = 0.1; 169 fr = 0.1; 169 rangeinit = 1.e+21; 170 rangeinit = 1.e+21; 170 geombig = 1.e+50*mm; 171 geombig = 1.e+50*mm; 171 geomlimit = geombig; 172 geomlimit = geombig; 172 tgeom = geombig; 173 tgeom = geombig; 173 tlimit = 1.e+10*mm; 174 tlimit = 1.e+10*mm; 174 presafety = 0.*mm; 175 presafety = 0.*mm; 175 // 176 // 176 particle = nullptr; << 177 particle = 0; 177 theManager = G4LossTableManager: 178 theManager = G4LossTableManager::Instance(); 178 firstStep = true; 179 firstStep = true; 179 currentKinEnergy = 0.0; 180 currentKinEnergy = 0.0; 180 currentRange = 0.0; 181 currentRange = 0.0; 181 // 182 // 182 tlimitminfix2 = 1.*nm; 183 tlimitminfix2 = 1.*nm; 183 tausmall = 1.e-16; 184 tausmall = 1.e-16; 184 mass = electron_mass_c2; 185 mass = electron_mass_c2; 185 taulim = 1.e-6; 186 taulim = 1.e-6; 186 // 187 // >> 188 facsafety = 0.6; >> 189 187 currentCouple = nullptr; 190 currentCouple = nullptr; 188 fParticleChange = nullptr; 191 fParticleChange = nullptr; 189 // 192 // 190 fZeff = 1.; 193 fZeff = 1.; 191 // 194 // 192 par1 = 0.; 195 par1 = 0.; 193 par2 = 0.; 196 par2 = 0.; 194 par3 = 0.; 197 par3 = 0.; 195 // 198 // 196 // Moliere screeing parameter will be used a 199 // Moliere screeing parameter will be used and (by default) corrections are 197 // appalied to the integrated quantities (sc 200 // appalied to the integrated quantities (screeing parameter, elastic mfp, first 198 // and second moments) derived from the corr 201 // and second moments) derived from the corresponding PWA quantities 199 // this PWA correction is ignored if Mott-co 202 // this PWA correction is ignored if Mott-correction is set to true because 200 // Mott-correction contains all these correc 203 // Mott-correction contains all these corrections as well 201 fIsUsePWACorrection = true; 204 fIsUsePWACorrection = true; 202 // 205 // 203 fIsUseMottCorrection = false; 206 fIsUseMottCorrection = false; 204 // 207 // 205 fLambda0 = 0.0; // elastic mea 208 fLambda0 = 0.0; // elastic mean free path 206 fLambda1 = 0.0; // first trans 209 fLambda1 = 0.0; // first transport mean free path 207 fScrA = 0.0; // screening p 210 fScrA = 0.0; // screening parameter 208 fG1 = 0.0; // first trans 211 fG1 = 0.0; // first transport coef. 209 // 212 // 210 fMCtoScrA = 1.0; 213 fMCtoScrA = 1.0; 211 fMCtoQ1 = 1.0; 214 fMCtoQ1 = 1.0; 212 fMCtoG2PerG1 = 1.0; 215 fMCtoG2PerG1 = 1.0; 213 // 216 // 214 fTheTrueStepLenght = 0.; 217 fTheTrueStepLenght = 0.; 215 fTheTransportDistance = 0.; 218 fTheTransportDistance = 0.; 216 fTheZPathLenght = 0.; 219 fTheZPathLenght = 0.; 217 // 220 // 218 fTheDisplacementVector.set(0.,0.,0.); 221 fTheDisplacementVector.set(0.,0.,0.); 219 fTheNewDirection.set(0.,0.,1.); 222 fTheNewDirection.set(0.,0.,1.); 220 // 223 // 221 fIsEverythingWasDone = false; 224 fIsEverythingWasDone = false; 222 fIsMultipleSacettring = false; 225 fIsMultipleSacettring = false; 223 fIsSingleScattering = false; 226 fIsSingleScattering = false; 224 fIsEndedUpOnBoundary = false; 227 fIsEndedUpOnBoundary = false; 225 fIsNoScatteringInMSC = false; 228 fIsNoScatteringInMSC = false; 226 fIsNoDisplace = false; 229 fIsNoDisplace = false; 227 fIsInsideSkin = false; 230 fIsInsideSkin = false; 228 fIsWasOnBoundary = false; 231 fIsWasOnBoundary = false; 229 fIsFirstRealStep = false; 232 fIsFirstRealStep = false; 230 rndmEngineMod = G4Random::getTheEng 233 rndmEngineMod = G4Random::getTheEngine(); 231 // 234 // 232 fGSTable = nullptr; 235 fGSTable = nullptr; 233 fPWACorrection = nullptr; 236 fPWACorrection = nullptr; 234 } 237 } 235 238 236 239 237 G4GoudsmitSaundersonMscModel::~G4GoudsmitSaund 240 G4GoudsmitSaundersonMscModel::~G4GoudsmitSaundersonMscModel() { 238 if (IsMaster()) { 241 if (IsMaster()) { 239 if (fGSTable) { 242 if (fGSTable) { 240 delete fGSTable; 243 delete fGSTable; 241 fGSTable = nullptr; 244 fGSTable = nullptr; 242 } 245 } 243 if (fPWACorrection) { 246 if (fPWACorrection) { 244 delete fPWACorrection; 247 delete fPWACorrection; 245 fPWACorrection = nullptr; 248 fPWACorrection = nullptr; 246 } 249 } 247 } 250 } 248 } 251 } 249 252 250 253 251 void G4GoudsmitSaundersonMscModel::Initialise( 254 void G4GoudsmitSaundersonMscModel::Initialise(const G4ParticleDefinition* p, const G4DataVector&) { 252 SetParticle(p); 255 SetParticle(p); 253 InitialiseParameters(p); << 254 // -create GoudsmitSaundersonTable and init 256 // -create GoudsmitSaundersonTable and init its Mott-correction member if 255 // Mott-correction was required 257 // Mott-correction was required 256 if (IsMaster()) { 258 if (IsMaster()) { 257 // get the Mott-correction flag from EmPar 259 // get the Mott-correction flag from EmParameters 258 if (G4EmParameters::Instance()->UseMottCor 260 if (G4EmParameters::Instance()->UseMottCorrection()) { 259 fIsUseMottCorrection = true; 261 fIsUseMottCorrection = true; 260 } 262 } 261 // Mott-correction includes other way of P 263 // Mott-correction includes other way of PWA x-section corrections so deactivate it even if it was true 262 // when Mott-correction is activated by th 264 // when Mott-correction is activated by the user 263 if (fIsUseMottCorrection) { 265 if (fIsUseMottCorrection) { 264 fIsUsePWACorrection = false; 266 fIsUsePWACorrection = false; 265 } 267 } 266 // clear GS-table 268 // clear GS-table 267 if (fGSTable) { 269 if (fGSTable) { 268 delete fGSTable; 270 delete fGSTable; 269 fGSTable = nullptr; 271 fGSTable = nullptr; 270 } 272 } 271 // clear PWA corrections table if any 273 // clear PWA corrections table if any 272 if (fPWACorrection) { 274 if (fPWACorrection) { 273 delete fPWACorrection; 275 delete fPWACorrection; 274 fPWACorrection = nullptr; 276 fPWACorrection = nullptr; 275 } 277 } 276 // create GS-table 278 // create GS-table 277 G4bool isElectron = true; 279 G4bool isElectron = true; 278 if (p->GetPDGCharge()>0.) { 280 if (p->GetPDGCharge()>0.) { 279 isElectron = false; 281 isElectron = false; 280 } 282 } 281 fGSTable = new G4GoudsmitSaundersonTable(i 283 fGSTable = new G4GoudsmitSaundersonTable(isElectron); 282 // G4GSTable will be initialised: 284 // G4GSTable will be initialised: 283 // - Screened-Rutherford DCS based GS angu 285 // - Screened-Rutherford DCS based GS angular distributions will be loaded only if they are not there yet 284 // - Mott-correction will be initialised i 286 // - Mott-correction will be initialised if Mott-correction was requested to be used 285 fGSTable->SetOptionMottCorrection(fIsUseMo 287 fGSTable->SetOptionMottCorrection(fIsUseMottCorrection); 286 // - set PWA correction (correction to int 288 // - set PWA correction (correction to integrated quantites from Dirac-PWA) 287 fGSTable->SetOptionPWACorrection(fIsUsePWA 289 fGSTable->SetOptionPWACorrection(fIsUsePWACorrection); 288 // init 290 // init 289 fGSTable->Initialise(LowEnergyLimit(),High 291 fGSTable->Initialise(LowEnergyLimit(),HighEnergyLimit()); 290 // create PWA corrections table if it was 292 // create PWA corrections table if it was requested (and not disactivated because active Mott-correction) 291 if (fIsUsePWACorrection) { 293 if (fIsUsePWACorrection) { 292 fPWACorrection = new G4GSPWACorrections( 294 fPWACorrection = new G4GSPWACorrections(isElectron); 293 fPWACorrection->Initialise(); 295 fPWACorrection->Initialise(); 294 } 296 } 295 } 297 } 296 fParticleChange = GetParticleChangeForMSC(p) 298 fParticleChange = GetParticleChangeForMSC(p); 297 } 299 } 298 300 299 301 300 void G4GoudsmitSaundersonMscModel::InitialiseL 302 void G4GoudsmitSaundersonMscModel::InitialiseLocal(const G4ParticleDefinition*, G4VEmModel* masterModel) { 301 fGSTable = static_cast<G4Goud 303 fGSTable = static_cast<G4GoudsmitSaundersonMscModel*>(masterModel)->GetGSTable(); 302 fIsUseMottCorrection = static_cast<G4Goud 304 fIsUseMottCorrection = static_cast<G4GoudsmitSaundersonMscModel*>(masterModel)->GetOptionMottCorrection(); 303 fIsUsePWACorrection = static_cast<G4Goud 305 fIsUsePWACorrection = static_cast<G4GoudsmitSaundersonMscModel*>(masterModel)->GetOptionPWACorrection(); 304 fPWACorrection = static_cast<G4Goud 306 fPWACorrection = static_cast<G4GoudsmitSaundersonMscModel*>(masterModel)->GetPWACorrection(); 305 } 307 } 306 308 307 309 308 // computes macroscopic first transport cross 310 // computes macroscopic first transport cross section: used only in testing not during mc transport 309 G4double G4GoudsmitSaundersonMscModel::CrossSe 311 G4double G4GoudsmitSaundersonMscModel::CrossSectionPerVolume(const G4Material* mat, 310 const 312 const G4ParticleDefinition*, 311 G4dou 313 G4double kineticEnergy, 312 G4dou 314 G4double, 313 G4dou 315 G4double) { 314 G4double xsecTr1 = 0.; // cross section per 316 G4double xsecTr1 = 0.; // cross section per volume i.e. macroscopic 1st transport cross section >> 317 G4double efEnergy = kineticEnergy; 315 // 318 // 316 fLambda0 = 0.0; // elastic mean free path 319 fLambda0 = 0.0; // elastic mean free path 317 fLambda1 = 0.0; // first transport mean free 320 fLambda1 = 0.0; // first transport mean free path 318 fScrA = 0.0; // screening parameter 321 fScrA = 0.0; // screening parameter 319 fG1 = 0.0; // first transport coef. 322 fG1 = 0.0; // first transport coef. 320 // use Moliere's screening (with Mott-corret 323 // use Moliere's screening (with Mott-corretion if it was requested) 321 G4double efEnergy = std::max(kineticEnergy, << 324 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.*CLHEP::eV; 322 // total mometum square 325 // total mometum square 323 G4double pt2 = efEnergy*(efEnergy+2.0*el 326 G4double pt2 = efEnergy*(efEnergy+2.0*electron_mass_c2); 324 // beta square 327 // beta square 325 G4double beta2 = pt2/(pt2+electron_mass_c2 328 G4double beta2 = pt2/(pt2+electron_mass_c2*electron_mass_c2); 326 // current material index 329 // current material index 327 G4int matindx = (G4int)mat->GetIndex(); << 330 G4int matindx = mat->GetIndex(); 328 // Moliere's b_c 331 // Moliere's b_c 329 G4double bc = fGSTable->GetMoliereBc(ma 332 G4double bc = fGSTable->GetMoliereBc(matindx); 330 // get the Mott-correcton factors if Mott-co 333 // get the Mott-correcton factors if Mott-correcton was requested by the user 331 fMCtoScrA = 1.0; 334 fMCtoScrA = 1.0; 332 fMCtoQ1 = 1.0; 335 fMCtoQ1 = 1.0; 333 fMCtoG2PerG1 = 1.0; 336 fMCtoG2PerG1 = 1.0; 334 G4double scpCor = 1.0; 337 G4double scpCor = 1.0; 335 if (fIsUseMottCorrection) { 338 if (fIsUseMottCorrection) { 336 fGSTable->GetMottCorrectionFactors(G4Log(e 339 fGSTable->GetMottCorrectionFactors(G4Log(efEnergy), beta2, matindx, fMCtoScrA, fMCtoQ1, fMCtoG2PerG1); 337 // ! no scattering power correction since 340 // ! no scattering power correction since the current couple is not set before this interface method is called 338 // scpCor = fGSTable->ComputeScatteringPow 341 // scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy); 339 } else if (fIsUsePWACorrection) { 342 } else if (fIsUsePWACorrection) { 340 fPWACorrection->GetPWACorrectionFactors(G4 343 fPWACorrection->GetPWACorrectionFactors(G4Log(efEnergy), beta2, matindx, fMCtoScrA, fMCtoQ1, fMCtoG2PerG1); 341 // scpCor = fGSTable->ComputeScatteringPow 344 // scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy); 342 } 345 } 343 // screening parameter: 346 // screening parameter: 344 // - if Mott-corretioncorrection: the Screen 347 // - if Mott-corretioncorrection: the Screened-Rutherford times Mott-corretion DCS with this 345 // screening parameter gives back the (els 348 // screening parameter gives back the (elsepa) PWA first transport cross section 346 // - if PWA correction: he Screened-Rutherfo 349 // - if PWA correction: he Screened-Rutherford DCS with this screening parameter 347 // gives back the (elsepa) PWA first trans 350 // gives back the (elsepa) PWA first transport cross section 348 fScrA = fGSTable->GetMoliereXc2(matindx)/ 351 fScrA = fGSTable->GetMoliereXc2(matindx)/(4.0*pt2*bc)*fMCtoScrA; 349 // elastic mean free path in Geant4 internal 352 // elastic mean free path in Geant4 internal lenght units: the neglected (1+screening parameter) term is corrected 350 // (if Mott-corretion: the corrected screeni 353 // (if Mott-corretion: the corrected screening parameter is used for this (1+A) correction + Moliere b_c is also 351 // corrected with the screening parameter co 354 // corrected with the screening parameter correction) 352 fLambda0 = beta2*(1.+fScrA)*fMCtoScrA/bc/scp 355 fLambda0 = beta2*(1.+fScrA)*fMCtoScrA/bc/scpCor; 353 // first transport coefficient (if Mott-corr 356 // first transport coefficient (if Mott-corretion: the corrected screening parameter is used (it will be fully 354 // consistent with the one used during the p 357 // consistent with the one used during the pre-computation of the Mott-correted GS angular distributions)) 355 fG1 = 2.0*fScrA*((1.0+fScrA)*G4Log(1.0/ 358 fG1 = 2.0*fScrA*((1.0+fScrA)*G4Log(1.0/fScrA+1.0)-1.0); 356 // first transport mean free path 359 // first transport mean free path 357 fLambda1 = fLambda0/fG1; 360 fLambda1 = fLambda0/fG1; 358 xsecTr1 = 1./fLambda1; 361 xsecTr1 = 1./fLambda1; 359 return xsecTr1; 362 return xsecTr1; 360 } 363 } 361 364 362 365 363 // gives back the first transport mean free pa 366 // gives back the first transport mean free path in internal G4 units 364 G4double 367 G4double 365 G4GoudsmitSaundersonMscModel::GetTransportMean 368 G4GoudsmitSaundersonMscModel::GetTransportMeanFreePath(const G4ParticleDefinition* /*partdef*/, 366 369 G4double kineticEnergy) { 367 // kinetic energy is assumed to be in Geant4 370 // kinetic energy is assumed to be in Geant4 internal energy unit which is MeV 368 G4double efEnergy = kineticEnergy; 371 G4double efEnergy = kineticEnergy; 369 // 372 // 370 const G4Material* mat = currentCouple->GetM 373 const G4Material* mat = currentCouple->GetMaterial(); 371 // 374 // 372 fLambda0 = 0.0; // elastic mean free path 375 fLambda0 = 0.0; // elastic mean free path 373 fLambda1 = 0.0; // first transport mean free 376 fLambda1 = 0.0; // first transport mean free path 374 fScrA = 0.0; // screening parameter 377 fScrA = 0.0; // screening parameter 375 fG1 = 0.0; // first transport coef. 378 fG1 = 0.0; // first transport coef. 376 379 377 // use Moliere's screening (with Mott-corret 380 // use Moliere's screening (with Mott-corretion if it was requested) 378 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.* 381 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.*CLHEP::eV; 379 // total mometum square 382 // total mometum square 380 G4double pt2 = efEnergy*(efEnergy+2.0*el 383 G4double pt2 = efEnergy*(efEnergy+2.0*electron_mass_c2); 381 // beta square 384 // beta square 382 G4double beta2 = pt2/(pt2+electron_mass_c2 385 G4double beta2 = pt2/(pt2+electron_mass_c2*electron_mass_c2); 383 // current material index 386 // current material index 384 G4int matindx = (G4int)mat->GetIndex(); << 387 G4int matindx = mat->GetIndex(); 385 // Moliere's b_c 388 // Moliere's b_c 386 G4double bc = fGSTable->GetMoliereBc(ma 389 G4double bc = fGSTable->GetMoliereBc(matindx); 387 // get the Mott-correcton factors if Mott-co 390 // get the Mott-correcton factors if Mott-correcton was requested by the user 388 fMCtoScrA = 1.0; 391 fMCtoScrA = 1.0; 389 fMCtoQ1 = 1.0; 392 fMCtoQ1 = 1.0; 390 fMCtoG2PerG1 = 1.0; 393 fMCtoG2PerG1 = 1.0; 391 G4double scpCor = 1.0; 394 G4double scpCor = 1.0; 392 if (fIsUseMottCorrection) { 395 if (fIsUseMottCorrection) { 393 fGSTable->GetMottCorrectionFactors(G4Log(e 396 fGSTable->GetMottCorrectionFactors(G4Log(efEnergy), beta2, matindx, fMCtoScrA, fMCtoQ1, fMCtoG2PerG1); 394 scpCor = fGSTable->ComputeScatteringPowerC 397 scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy); 395 } else if (fIsUsePWACorrection) { 398 } else if (fIsUsePWACorrection) { 396 fPWACorrection->GetPWACorrectionFactors(G4 399 fPWACorrection->GetPWACorrectionFactors(G4Log(efEnergy), beta2, matindx, fMCtoScrA, fMCtoQ1, fMCtoG2PerG1); 397 // scpCor = fGSTable->ComputeScatteringPow 400 // scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy); 398 } 401 } 399 // screening parameter: 402 // screening parameter: 400 // - if Mott-corretioncorrection: the Screen 403 // - if Mott-corretioncorrection: the Screened-Rutherford times Mott-corretion DCS with this 401 // screening parameter gives back the (els 404 // screening parameter gives back the (elsepa) PWA first transport cross section 402 // - if PWA correction: he Screened-Rutherfo 405 // - if PWA correction: he Screened-Rutherford DCS with this screening parameter 403 // gives back the (elsepa) PWA first trans 406 // gives back the (elsepa) PWA first transport cross section 404 fScrA = fGSTable->GetMoliereXc2(matindx)/ 407 fScrA = fGSTable->GetMoliereXc2(matindx)/(4.0*pt2*bc)*fMCtoScrA; 405 // elastic mean free path in Geant4 internal 408 // elastic mean free path in Geant4 internal lenght units: the neglected (1+screening parameter) term is corrected 406 // (if Mott-corretion: the corrected screeni 409 // (if Mott-corretion: the corrected screening parameter is used for this (1+A) correction + Moliere b_c is also 407 // corrected with the screening parameter co 410 // corrected with the screening parameter correction) 408 fLambda0 = beta2*(1.+fScrA)*fMCtoScrA/bc/scp 411 fLambda0 = beta2*(1.+fScrA)*fMCtoScrA/bc/scpCor; 409 // first transport coefficient (if Mott-corr 412 // first transport coefficient (if Mott-corretion: the corrected screening parameter is used (it will be fully 410 // consistent with the one used during the p 413 // consistent with the one used during the pre-computation of the Mott-correted GS angular distributions)) 411 fG1 = 2.0*fScrA*((1.0+fScrA)*G4Log(1.0/ 414 fG1 = 2.0*fScrA*((1.0+fScrA)*G4Log(1.0/fScrA+1.0)-1.0); 412 // first transport mean free path 415 // first transport mean free path 413 fLambda1 = fLambda0/fG1; 416 fLambda1 = fLambda0/fG1; 414 417 415 return fLambda1; 418 return fLambda1; 416 } 419 } 417 420 418 421 419 G4double 422 G4double 420 G4GoudsmitSaundersonMscModel::GetTransportMean 423 G4GoudsmitSaundersonMscModel::GetTransportMeanFreePathOnly(const G4ParticleDefinition* /*partdef*/, 421 424 G4double kineticEnergy) { 422 // kinetic energy is assumed to be in Geant4 425 // kinetic energy is assumed to be in Geant4 internal energy unit which is MeV 423 G4double efEnergy = kineticEnergy; 426 G4double efEnergy = kineticEnergy; 424 // 427 // 425 const G4Material* mat = currentCouple->GetM 428 const G4Material* mat = currentCouple->GetMaterial(); 426 // 429 // 427 G4double lambda0 = 0.0; // elastc mean free 430 G4double lambda0 = 0.0; // elastc mean free path 428 G4double lambda1 = 0.0; // first transport m 431 G4double lambda1 = 0.0; // first transport mean free path 429 G4double scrA = 0.0; // screening paramet 432 G4double scrA = 0.0; // screening parametr 430 G4double g1 = 0.0; // first transport m 433 G4double g1 = 0.0; // first transport mean free path 431 434 432 // use Moliere's screening (with Mott-corret 435 // use Moliere's screening (with Mott-corretion if it was requested) 433 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.* 436 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.*CLHEP::eV; 434 // total mometum square in Geant4 internal e 437 // total mometum square in Geant4 internal energy2 units which is MeV2 435 G4double pt2 = efEnergy*(efEnergy+2.0*el 438 G4double pt2 = efEnergy*(efEnergy+2.0*electron_mass_c2); 436 G4double beta2 = pt2/(pt2+electron_mass_c2 439 G4double beta2 = pt2/(pt2+electron_mass_c2*electron_mass_c2); 437 G4int matindx = (G4int)mat->GetIndex(); << 440 G4int matindx = mat->GetIndex(); 438 G4double bc = fGSTable->GetMoliereBc(ma 441 G4double bc = fGSTable->GetMoliereBc(matindx); 439 // get the Mott-correcton factors if Mott-co 442 // get the Mott-correcton factors if Mott-correcton was requested by the user 440 G4double mctoScrA = 1.0; 443 G4double mctoScrA = 1.0; 441 G4double mctoQ1 = 1.0; 444 G4double mctoQ1 = 1.0; 442 G4double mctoG2PerG1 = 1.0; 445 G4double mctoG2PerG1 = 1.0; 443 G4double scpCor = 1.0; 446 G4double scpCor = 1.0; 444 if (fIsUseMottCorrection) { 447 if (fIsUseMottCorrection) { 445 fGSTable->GetMottCorrectionFactors(G4Log(e 448 fGSTable->GetMottCorrectionFactors(G4Log(efEnergy), beta2, matindx, mctoScrA, mctoQ1, mctoG2PerG1); 446 scpCor = fGSTable->ComputeScatteringPowerC 449 scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy); 447 } else if (fIsUsePWACorrection) { 450 } else if (fIsUsePWACorrection) { 448 fPWACorrection->GetPWACorrectionFactors(G4 451 fPWACorrection->GetPWACorrectionFactors(G4Log(efEnergy), beta2, matindx, mctoScrA, mctoQ1, mctoG2PerG1); 449 // scpCor = fGSTable->ComputeScatteringPow 452 // scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy); 450 } 453 } 451 scrA = fGSTable->GetMoliereXc2(matindx)/( 454 scrA = fGSTable->GetMoliereXc2(matindx)/(4.0*pt2*bc)*mctoScrA; 452 // total elastic mean free path in Geant4 in 455 // total elastic mean free path in Geant4 internal lenght units 453 lambda0 = beta2*(1.+scrA)*mctoScrA/bc/scpCor 456 lambda0 = beta2*(1.+scrA)*mctoScrA/bc/scpCor; 454 g1 = 2.0*scrA*((1.0+scrA)*G4Log(1.0/scr 457 g1 = 2.0*scrA*((1.0+scrA)*G4Log(1.0/scrA+1.0)-1.0); 455 lambda1 = lambda0/g1; 458 lambda1 = lambda0/g1; 456 459 457 return lambda1; 460 return lambda1; 458 } 461 } 459 462 460 463 461 void G4GoudsmitSaundersonMscModel::StartTracki 464 void G4GoudsmitSaundersonMscModel::StartTracking(G4Track* track) { 462 SetParticle(track->GetDynamicParticle()->Get 465 SetParticle(track->GetDynamicParticle()->GetDefinition()); 463 firstStep = true; 466 firstStep = true; 464 tlimit = tgeom = rangeinit = geombig; 467 tlimit = tgeom = rangeinit = geombig; 465 rangeinit = 1.e+21; 468 rangeinit = 1.e+21; 466 } 469 } 467 470 468 471 469 G4double G4GoudsmitSaundersonMscModel::Compute 472 G4double G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit(const G4Track& track, 470 473 G4double& currentMinimalStep) { 471 G4double skindepth = 0.; 474 G4double skindepth = 0.; 472 // 475 // 473 const G4DynamicParticle* dp = track.GetDynam 476 const G4DynamicParticle* dp = track.GetDynamicParticle(); 474 G4StepPoint* sp = track.GetStep( 477 G4StepPoint* sp = track.GetStep()->GetPreStepPoint(); 475 G4StepStatus stepStatus = sp->GetStepSta 478 G4StepStatus stepStatus = sp->GetStepStatus(); 476 currentCouple = track.GetMater 479 currentCouple = track.GetMaterialCutsCouple(); 477 SetCurrentCouple(currentCouple); 480 SetCurrentCouple(currentCouple); 478 currentMaterialIndex = (G4int)current << 481 currentMaterialIndex = currentCouple->GetMaterial()->GetIndex(); 479 currentKinEnergy = dp->GetKinetic 482 currentKinEnergy = dp->GetKineticEnergy(); 480 currentRange = GetRange(parti << 483 currentRange = GetRange(particle,currentKinEnergy,currentCouple); 481 dp->G << 482 // elastic and first transport mfp, screenin 484 // elastic and first transport mfp, screening parameter and G1 are also set 483 // (Mott-correction will be used if it was r 485 // (Mott-correction will be used if it was requested by the user) 484 fLambda1 = GetTransportMeanFreePath(particle 486 fLambda1 = GetTransportMeanFreePath(particle,currentKinEnergy); 485 // Set initial values: 487 // Set initial values: 486 // : lengths are initialised to currentMini 488 // : lengths are initialised to currentMinimalStep which is the true, minimum 487 // step length from all other physics 489 // step length from all other physics 488 fTheTrueStepLenght = currentMinimalStep; 490 fTheTrueStepLenght = currentMinimalStep; 489 fTheTransportDistance = currentMinimalStep; 491 fTheTransportDistance = currentMinimalStep; 490 fTheZPathLenght = currentMinimalStep; 492 fTheZPathLenght = currentMinimalStep; // will need to be converted 491 fTheDisplacementVector.set(0.,0.,0.); 493 fTheDisplacementVector.set(0.,0.,0.); 492 fTheNewDirection.set(0.,0.,1.); 494 fTheNewDirection.set(0.,0.,1.); 493 495 494 // Can everything be done in the step limit 496 // Can everything be done in the step limit phase ? 495 fIsEverythingWasDone = false; 497 fIsEverythingWasDone = false; 496 // Multiple scattering needs to be sample ? 498 // Multiple scattering needs to be sample ? 497 fIsMultipleSacettring = false; 499 fIsMultipleSacettring = false; 498 // Single scattering needs to be sample ? 500 // Single scattering needs to be sample ? 499 fIsSingleScattering = false; 501 fIsSingleScattering = false; 500 // Was zero deflection in multiple scatterin 502 // Was zero deflection in multiple scattering sampling ? 501 fIsNoScatteringInMSC = false; 503 fIsNoScatteringInMSC = false; 502 // Do not care about displacement in MSC sam 504 // Do not care about displacement in MSC sampling 503 // ( used only in the case of gIsOptimizatio 505 // ( used only in the case of gIsOptimizationOn = true) 504 fIsNoDisplace = false; 506 fIsNoDisplace = false; 505 // get pre-step point safety 507 // get pre-step point safety 506 presafety = sp->GetSafety(); 508 presafety = sp->GetSafety(); 507 // 509 // 508 fZeff = currentCouple->GetMaterial()->GetIon 510 fZeff = currentCouple->GetMaterial()->GetIonisation()->GetZeffective(); 509 // distance will take into account max-fluct 511 // distance will take into account max-fluct. 510 G4double distance = currentRange; 512 G4double distance = currentRange; 511 distance *= (1.20-fZeff*(1.62e-2-9.22e-5*fZe 513 distance *= (1.20-fZeff*(1.62e-2-9.22e-5*fZeff)); 512 // 514 // 513 // Possible optimization : if the distance i 515 // Possible optimization : if the distance is samller than the safety -> the 514 // particle will never leave this volume -> 516 // particle will never leave this volume -> dispalcement 515 // as the effect of multiple elastic scatter 517 // as the effect of multiple elastic scattering can be skipped 516 // Important : this optimization can cause p 518 // Important : this optimization can cause problems if one does scoring 517 // in a bigger volume since MSC won't be don 519 // in a bigger volume since MSC won't be done deep inside the volume when 518 // distance < safety so don't use optimized- 520 // distance < safety so don't use optimized-mode in such case. 519 if (gIsOptimizationOn && (distance<presafety 521 if (gIsOptimizationOn && (distance<presafety)) { 520 // Indicate that we need to do MSC after 522 // Indicate that we need to do MSC after transportation and no dispalcement. 521 fIsMultipleSacettring = true; 523 fIsMultipleSacettring = true; 522 fIsNoDisplace = true; 524 fIsNoDisplace = true; 523 } else if (steppingAlgorithm==fUseDistanceTo 525 } else if (steppingAlgorithm==fUseDistanceToBoundary) { 524 //Compute geomlimit (and presafety) : 526 //Compute geomlimit (and presafety) : 525 // - geomlimit will be: 527 // - geomlimit will be: 526 // == the straight line distance to the 528 // == the straight line distance to the boundary if currentRange is 527 // longer than that 529 // longer than that 528 // == a big value [geombig = 1.e50*mm] 530 // == a big value [geombig = 1.e50*mm] if currentRange is shorter than 529 // the straight line distance to the 531 // the straight line distance to the boundary 530 // - presafety will be updated as well 532 // - presafety will be updated as well 531 // So the particle can travell 'gemlimit' 533 // So the particle can travell 'gemlimit' distance (along a straight 532 // line!) in its current direction: 534 // line!) in its current direction: 533 // (1) before reaching a boundary (geomli 535 // (1) before reaching a boundary (geomlimit < geombig) OR 534 // (2) before reaching its current range 536 // (2) before reaching its current range (geomlimit == geombig) 535 geomlimit = ComputeGeomLimit(track, presaf 537 geomlimit = ComputeGeomLimit(track, presafety, currentRange); 536 // Record that the particle is on a bounda 538 // Record that the particle is on a boundary 537 if ( (stepStatus==fGeomBoundary) || (stepS 539 if ( (stepStatus==fGeomBoundary) || (stepStatus==fUndefined && presafety==0.0)) { 538 fIsWasOnBoundary = true; 540 fIsWasOnBoundary = true; 539 } 541 } 540 // Set skin depth = skin x elastic_mean_fr 542 // Set skin depth = skin x elastic_mean_free_path 541 skindepth = skin*fLambda0; 543 skindepth = skin*fLambda0; 542 // Init the flag that indicates that the p 544 // Init the flag that indicates that the particle are within a skindepth 543 // distance from a boundary 545 // distance from a boundary 544 fIsInsideSkin = false; 546 fIsInsideSkin = false; 545 // Check if we can try Single Scattering b 547 // Check if we can try Single Scattering because we are within skindepth 546 // distance from/to a boundary OR the curr 548 // distance from/to a boundary OR the current minimum true-step-length is 547 // shorter than skindepth. NOTICE: the lat 549 // shorter than skindepth. NOTICE: the latest has only efficieny reasons 548 // because the MSC angular sampling is fin 550 // because the MSC angular sampling is fine for any short steps but much 549 // faster to try single scattering in case 551 // faster to try single scattering in case of short steps. 550 if ((stepStatus==fGeomBoundary) || (presaf 552 if ((stepStatus==fGeomBoundary) || (presafety<skindepth) || (fTheTrueStepLenght<skindepth)) { 551 // check if we are within skindepth dist 553 // check if we are within skindepth distance from a boundary 552 if ((stepStatus == fGeomBoundary) || (pr 554 if ((stepStatus == fGeomBoundary) || (presafety < skindepth)) { 553 fIsInsideSkin = true; 555 fIsInsideSkin = true; 554 fIsWasOnBoundary = true; 556 fIsWasOnBoundary = true; 555 } 557 } 556 //Try single scattering: 558 //Try single scattering: 557 // - sample distance to next single scat 559 // - sample distance to next single scattering interaction (sslimit) 558 // - compare to current minimum length 560 // - compare to current minimum length 559 // == if sslimit is the shorter: 561 // == if sslimit is the shorter: 560 // - set the step length to ssl 562 // - set the step length to sslimit 561 // - indicate that single scatt 563 // - indicate that single scattering needs to be done 562 // == else : nothing to do 564 // == else : nothing to do 563 //- in both cases, the step length was v 565 //- in both cases, the step length was very short so geometrical and 564 // true path length are the same 566 // true path length are the same 565 G4double sslimit = -1.*fLambda0*G4Log(G4 567 G4double sslimit = -1.*fLambda0*G4Log(G4UniformRand()); 566 // compare to current minimum step lengt 568 // compare to current minimum step length 567 if (sslimit<fTheTrueStepLenght) { 569 if (sslimit<fTheTrueStepLenght) { 568 fTheTrueStepLenght = sslimit; 570 fTheTrueStepLenght = sslimit; 569 fIsSingleScattering = true; 571 fIsSingleScattering = true; 570 } 572 } 571 // short step -> true step length equal 573 // short step -> true step length equal to geometrical path length 572 fTheZPathLenght = fTheTrueStepLeng 574 fTheZPathLenght = fTheTrueStepLenght; 573 // Set taht everything is done in step-l 575 // Set taht everything is done in step-limit phase so no MSC call 574 // We will check if we need to perform t 576 // We will check if we need to perform the single-scattering angular 575 // sampling i.e. if single elastic scatt 577 // sampling i.e. if single elastic scattering was the winer! 576 fIsEverythingWasDone = true; 578 fIsEverythingWasDone = true; 577 } else { 579 } else { 578 // After checking we know that we cannot 580 // After checking we know that we cannot try single scattering so we will 579 // need to make an MSC step 581 // need to make an MSC step 580 // Indicate that we need to make and MSC 582 // Indicate that we need to make and MSC step. We do not check if we can 581 // do it now i.e. if presafety>final_tru 583 // do it now i.e. if presafety>final_true_step_length so we let the 582 // fIsEverythingWasDone = false which in 584 // fIsEverythingWasDone = false which indicates that we will perform 583 // MSC after transportation. 585 // MSC after transportation. 584 fIsMultipleSacettring = true; 586 fIsMultipleSacettring = true; 585 // Init the first-real-step falg: it wil 587 // Init the first-real-step falg: it will indicate if we do the first 586 // non-single scattering step in this vo 588 // non-single scattering step in this volume with this particle 587 fIsFirstRealStep = false; 589 fIsFirstRealStep = false; 588 // If previously the partcile was on bou 590 // If previously the partcile was on boundary it was within skin as 589 // well. When it is not within skin anym 591 // well. When it is not within skin anymore it has just left the skin 590 // so we make the first real MSC step wi 592 // so we make the first real MSC step with the particle. 591 if (fIsWasOnBoundary && !fIsInsideSkin) 593 if (fIsWasOnBoundary && !fIsInsideSkin) { 592 // reset the 'was on boundary' indicat 594 // reset the 'was on boundary' indicator flag 593 fIsWasOnBoundary = false; 595 fIsWasOnBoundary = false; 594 fIsFirstRealStep = true; 596 fIsFirstRealStep = true; 595 } 597 } 596 // If this is the first-real msc step (t 598 // If this is the first-real msc step (the partcile has just left the 597 // skin) or this is the first step with 599 // skin) or this is the first step with the particle (was born or 598 // primary): 600 // primary): 599 // - set the initial range that will b 601 // - set the initial range that will be used later to limit its step 600 // (only in this volume, because aft 602 // (only in this volume, because after boundary crossing at the 601 // first-real MSC step we will reset 603 // first-real MSC step we will reset) 602 // - don't let the partcile to cross th 604 // - don't let the partcile to cross the volume just in one step 603 if (firstStep || fIsFirstRealStep || ran 605 if (firstStep || fIsFirstRealStep || rangeinit>1.e+20) { 604 rangeinit = currentRange; 606 rangeinit = currentRange; 605 // If geomlimit < geombig than the par 607 // If geomlimit < geombig than the particle might reach the boundary 606 // along its initial direction before 608 // along its initial direction before losing its energy (in this step) 607 // Otherwise we can be sure that the p 609 // Otherwise we can be sure that the particle will lose it energy 608 // before reaching the boundary along 610 // before reaching the boundary along a starigth line so there is no 609 // geometrical limit appalied. [Howeve 611 // geometrical limit appalied. [However, tgeom is set only in the 610 // first or the first-real MSC step. A 612 // first or the first-real MSC step. After the first or first real 611 // MSC step the direction will change 613 // MSC step the direction will change tgeom won't guaranty anything! 612 // But we will try to end up within sk 614 // But we will try to end up within skindepth from the boundary using 613 // the actual value of geomlimit(See l 615 // the actual value of geomlimit(See later at step reduction close to 614 // boundary).] 616 // boundary).] 615 if (geomlimit<geombig) { 617 if (geomlimit<geombig) { 616 // transfrom straight line distance 618 // transfrom straight line distance to the boundary to real step 617 // length based on the mean values ( 619 // length based on the mean values (using the prestep point 618 // first-transport mean free path i. 620 // first-transport mean free path i.e. no energy loss correction) 619 if ((1.-geomlimit/fLambda1)> 0.) { 621 if ((1.-geomlimit/fLambda1)> 0.) { 620 geomlimit = -fLambda1*G4Log(1.-geo 622 geomlimit = -fLambda1*G4Log(1.-geomlimit/fLambda1); 621 } 623 } 622 // the 2-different case that could l 624 // the 2-different case that could lead us here 623 if (firstStep) { 625 if (firstStep) { 624 tgeom = 2.*geomlimit/facgeom; 626 tgeom = 2.*geomlimit/facgeom; 625 } else { 627 } else { 626 tgeom = geomlimit/facgeom; 628 tgeom = geomlimit/facgeom; 627 } 629 } 628 } else { 630 } else { 629 tgeom = geombig; 631 tgeom = geombig; 630 } 632 } 631 } 633 } 632 // True step length limit from range fac 634 // True step length limit from range factor. Noteice, that the initial 633 // range is used that was set at the fir 635 // range is used that was set at the first step or first-real MSC step 634 // in this volume with this particle. 636 // in this volume with this particle. 635 tlimit = facrange*rangeinit; 637 tlimit = facrange*rangeinit; 636 // Take the minimum of the true step len 638 // Take the minimum of the true step length limits coming from 637 // geometrical constraint or range-facto 639 // geometrical constraint or range-factor limitation 638 tlimit = std::min(tlimit,tgeom); 640 tlimit = std::min(tlimit,tgeom); 639 // Step reduction close to boundary: we 641 // Step reduction close to boundary: we try to end up within skindepth 640 // from the boundary ( Notice: in case o 642 // from the boundary ( Notice: in case of mag. field it might not work 641 // because geomlimit is the straigth lin 643 // because geomlimit is the straigth line distance to the boundary in 642 // the currect direction (if geomlimit<g 644 // the currect direction (if geomlimit<geombig) and mag. field can 643 // change the initial direction. So te p 645 // change the initial direction. So te particle might hit some boundary 644 // before in a different direction. Howe 646 // before in a different direction. However, here we restrict the true 645 // path length to this (straight line) l 647 // path length to this (straight line) lenght so the corresponding 646 // transport distance (straight line) wi 648 // transport distance (straight line) will be even shorter than 647 // geomlimit-0.999*skindepth after the c 649 // geomlimit-0.999*skindepth after the change of true->geom. 648 if (geomlimit<geombig) { 650 if (geomlimit<geombig) { 649 tlimit = std::min(tlimit, geomlimit-0. 651 tlimit = std::min(tlimit, geomlimit-0.999*skindepth); 650 } 652 } 651 // randomize 1st step or 1st 'normal' st 653 // randomize 1st step or 1st 'normal' step in volume 652 if (firstStep || fIsFirstRealStep) { 654 if (firstStep || fIsFirstRealStep) { 653 fTheTrueStepLenght = std::min(fTheTrue 655 fTheTrueStepLenght = std::min(fTheTrueStepLenght, Randomizetlimit()); 654 } else { 656 } else { 655 fTheTrueStepLenght = std::min(fTheTrue 657 fTheTrueStepLenght = std::min(fTheTrueStepLenght, tlimit); 656 } 658 } 657 } 659 } 658 } else if (steppingAlgorithm==fUseSafetyPlus 660 } else if (steppingAlgorithm==fUseSafetyPlus) { // THE ERROR_FREE stepping alg. 659 presafety = ComputeSafety(sp->GetPosition 661 presafety = ComputeSafety(sp->GetPosition(),fTheTrueStepLenght); 660 geomlimit = presafety; 662 geomlimit = presafety; 661 // Set skin depth = skin x elastic_mean_fr 663 // Set skin depth = skin x elastic_mean_free_path 662 skindepth = skin*fLambda0; 664 skindepth = skin*fLambda0; 663 // Check if we can try Single Scattering b 665 // Check if we can try Single Scattering because we are within skindepth 664 // distance from/to a boundary OR the curr 666 // distance from/to a boundary OR the current minimum true-step-length is 665 // shorter than skindepth. NOTICE: the lat 667 // shorter than skindepth. NOTICE: the latest has only efficieny reasons 666 // because the MSC angular sampling is fin 668 // because the MSC angular sampling is fine for any short steps but much 667 // faster to try single scattering in case 669 // faster to try single scattering in case of short steps. 668 if ((stepStatus==fGeomBoundary) || (presaf 670 if ((stepStatus==fGeomBoundary) || (presafety<skindepth) || (fTheTrueStepLenght<skindepth)) { 669 //Try single scattering: 671 //Try single scattering: 670 // - sample distance to next single scat 672 // - sample distance to next single scattering interaction (sslimit) 671 // - compare to current minimum length 673 // - compare to current minimum length 672 // == if sslimit is the shorter: 674 // == if sslimit is the shorter: 673 // - set the step length to ssl 675 // - set the step length to sslimit 674 // - indicate that single scatt 676 // - indicate that single scattering needs to be done 675 // == else : nothing to do 677 // == else : nothing to do 676 //- in both cases, the step length was v 678 //- in both cases, the step length was very short so geometrical and 677 // true path length are the same 679 // true path length are the same 678 G4double sslimit = -1.*fLambda0*G4Log(G4 680 G4double sslimit = -1.*fLambda0*G4Log(G4UniformRand()); 679 // compare to current minimum step lengt 681 // compare to current minimum step length 680 if (sslimit<fTheTrueStepLenght) { 682 if (sslimit<fTheTrueStepLenght) { 681 fTheTrueStepLenght = sslimit; 683 fTheTrueStepLenght = sslimit; 682 fIsSingleScattering = true; 684 fIsSingleScattering = true; 683 } 685 } 684 // short step -> true step length equal 686 // short step -> true step length equal to geometrical path length 685 fTheZPathLenght = fTheTrueStepLeng 687 fTheZPathLenght = fTheTrueStepLenght; 686 // Set taht everything is done in step-l 688 // Set taht everything is done in step-limit phase so no MSC call 687 // We will check if we need to perform t 689 // We will check if we need to perform the single-scattering angular 688 // sampling i.e. if single elastic scatt 690 // sampling i.e. if single elastic scattering was the winer! 689 fIsEverythingWasDone = true; 691 fIsEverythingWasDone = true; 690 } else { 692 } else { 691 // After checking we know that we cannot 693 // After checking we know that we cannot try single scattering so we will 692 // need to make an MSC step 694 // need to make an MSC step 693 // Indicate that we need to make and MSC 695 // Indicate that we need to make and MSC step. 694 fIsMultipleSacettring = true; 696 fIsMultipleSacettring = true; 695 fIsEverythingWasDone = true; 697 fIsEverythingWasDone = true; 696 // limit from range factor 698 // limit from range factor 697 fTheTrueStepLenght = std::min(fTheTru 699 fTheTrueStepLenght = std::min(fTheTrueStepLenght, facrange*currentRange); 698 // never let the particle go further tha 700 // never let the particle go further than the safety if we are out of the skin 699 // if we are here we are out of the skin 701 // if we are here we are out of the skin, presafety > 0. 700 if (fTheTrueStepLenght>presafety) { 702 if (fTheTrueStepLenght>presafety) { 701 fTheTrueStepLenght = std::min(fTheTrue 703 fTheTrueStepLenght = std::min(fTheTrueStepLenght, presafety); 702 } 704 } 703 // make sure that we are still within th 705 // make sure that we are still within the aplicability of condensed histry model 704 // i.e. true step length is not longer t 706 // i.e. true step length is not longer than first transport mean free path. 705 // We schould take into account energy l 707 // We schould take into account energy loss along 0.5x lambda_transport1 706 // step length as well. So let it 0.5 x 708 // step length as well. So let it 0.5 x lambda_transport1 707 fTheTrueStepLenght = std::min(fTheTrueSt 709 fTheTrueStepLenght = std::min(fTheTrueStepLenght, fLambda1*0.5); 708 } 710 } 709 } else { 711 } else { 710 // This is the default stepping algorithm: 712 // This is the default stepping algorithm: the fastest but the least 711 // accurate that corresponds to fUseSafety 713 // accurate that corresponds to fUseSafety in Urban model. Note, that GS 712 // model can handle any short steps so we 714 // model can handle any short steps so we do not need the minimum limits 713 // 715 // 714 // NO single scattering in case of skin or 716 // NO single scattering in case of skin or short steps (by defult the MSC 715 // model will be single or even no scatter 717 // model will be single or even no scattering in case of short steps 716 // compared to the elastic mean free path. 718 // compared to the elastic mean free path.) 717 // 719 // 718 // indicate that MSC needs to be done (alw 720 // indicate that MSC needs to be done (always and always after transportation) 719 fIsMultipleSacettring = true; 721 fIsMultipleSacettring = true; 720 if (stepStatus!=fGeomBoundary) { 722 if (stepStatus!=fGeomBoundary) { 721 presafety = ComputeSafety(sp->GetPositio 723 presafety = ComputeSafety(sp->GetPosition(),fTheTrueStepLenght); 722 } 724 } 723 // Far from boundary-> in optimized mode d 725 // Far from boundary-> in optimized mode do not sample dispalcement. 724 if ((distance<presafety) && (gIsOptimizati 726 if ((distance<presafety) && (gIsOptimizationOn)) { 725 fIsNoDisplace = true; 727 fIsNoDisplace = true; 726 } else { 728 } else { 727 // Urban like 729 // Urban like 728 if (firstStep || (stepStatus==fGeomBound 730 if (firstStep || (stepStatus==fGeomBoundary) || rangeinit>1.e+20) { 729 rangeinit = currentRange; 731 rangeinit = currentRange; 730 fr = facrange; 732 fr = facrange; 731 // We don't use this: we won't converge to the 733 // We don't use this: we won't converge to the single scattering results with 732 // decreasing range-factor. 734 // decreasing range-factor. 733 // rangeinit = std::max(rangeinit 735 // rangeinit = std::max(rangeinit, fLambda1); 734 // if(fLambda1 > lambdalimit) { 736 // if(fLambda1 > lambdalimit) { 735 // fr *= (0.75+0.25*fLambda1/la 737 // fr *= (0.75+0.25*fLambda1/lambdalimit); 736 // } 738 // } 737 739 738 } 740 } 739 //step limit 741 //step limit 740 tlimit = std::max(fr*rangeinit, facsafet 742 tlimit = std::max(fr*rangeinit, facsafety*presafety); 741 // first step randomization 743 // first step randomization 742 if (firstStep || stepStatus==fGeomBounda 744 if (firstStep || stepStatus==fGeomBoundary) { 743 fTheTrueStepLenght = std::min(fTheTrue 745 fTheTrueStepLenght = std::min(fTheTrueStepLenght, Randomizetlimit()); 744 } else { 746 } else { 745 fTheTrueStepLenght = std::min(fTheTrue 747 fTheTrueStepLenght = std::min(fTheTrueStepLenght, tlimit); 746 } 748 } 747 } 749 } 748 } 750 } 749 // 751 // 750 // unset first-step 752 // unset first-step 751 firstStep =false; 753 firstStep =false; 752 // performe single scattering, multiple scat 754 // performe single scattering, multiple scattering if this later can be done safely here 753 if (fIsEverythingWasDone) { 755 if (fIsEverythingWasDone) { 754 if (fIsSingleScattering) { 756 if (fIsSingleScattering) { 755 // sample single scattering 757 // sample single scattering 756 //G4double ekin = 0.5*(currentKinEnerg 758 //G4double ekin = 0.5*(currentKinEnergy + GetEnergy(particle,currentRange-fTheTrueStepLenght,currentCouple)); 757 G4double lekin = G4Log(currentKinEnergy 759 G4double lekin = G4Log(currentKinEnergy); 758 G4double pt2 = currentKinEnergy*(curr 760 G4double pt2 = currentKinEnergy*(currentKinEnergy+2.0*CLHEP::electron_mass_c2); 759 G4double beta2 = pt2/(pt2+CLHEP::electr 761 G4double beta2 = pt2/(pt2+CLHEP::electron_mass_c2*CLHEP::electron_mass_c2); 760 G4double cost = fGSTable->SingleScatte 762 G4double cost = fGSTable->SingleScattering(1., fScrA, lekin, beta2, currentMaterialIndex); 761 // protection 763 // protection 762 if (cost<-1.) cost = -1.; 764 if (cost<-1.) cost = -1.; 763 if (cost> 1.) cost = 1.; 765 if (cost> 1.) cost = 1.; 764 // compute sint 766 // compute sint 765 G4double dum = 1.-cost; 767 G4double dum = 1.-cost; 766 G4double sint = std::sqrt(dum*(2.-dum) 768 G4double sint = std::sqrt(dum*(2.-dum)); 767 G4double phi = CLHEP::twopi*G4Uniform 769 G4double phi = CLHEP::twopi*G4UniformRand(); 768 G4double sinPhi = std::sin(phi); 770 G4double sinPhi = std::sin(phi); 769 G4double cosPhi = std::cos(phi); 771 G4double cosPhi = std::cos(phi); 770 fTheNewDirection.set(sint*cosPhi,sint*si 772 fTheNewDirection.set(sint*cosPhi,sint*sinPhi,cost); 771 } else if (fIsMultipleSacettring) { 773 } else if (fIsMultipleSacettring) { 772 // sample multiple scattering 774 // sample multiple scattering 773 SampleMSC(); // fTheZPathLenght, fTheDis 775 SampleMSC(); // fTheZPathLenght, fTheDisplacementVector and fTheNewDirection will be set 774 } // and if single scattering but it was l 776 } // and if single scattering but it was longer => nothing to do 775 } //else { do nothing here but after transpo 777 } //else { do nothing here but after transportation 776 // 778 // 777 return ConvertTrueToGeom(fTheTrueStepLenght, 779 return ConvertTrueToGeom(fTheTrueStepLenght,currentMinimalStep); 778 } 780 } 779 781 780 782 781 G4double G4GoudsmitSaundersonMscModel::Compute 783 G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength(G4double) { 782 // convert true ->geom 784 // convert true ->geom 783 // It is called from the step limitation Com 785 // It is called from the step limitation ComputeTruePathLengthLimit if 784 // !fIsEverythingWasDone but protect: 786 // !fIsEverythingWasDone but protect: 785 par1 = -1.; 787 par1 = -1.; 786 par2 = par3 = 0.; 788 par2 = par3 = 0.; 787 // if fIsEverythingWasDone = TRUE => fTheZP 789 // if fIsEverythingWasDone = TRUE => fTheZPathLenght is already set 788 // so return with the already known value 790 // so return with the already known value 789 // Otherwise: 791 // Otherwise: 790 if (!fIsEverythingWasDone) { 792 if (!fIsEverythingWasDone) { 791 // this correction needed to run MSC with 793 // this correction needed to run MSC with eIoni and eBrem inactivated 792 // and makes no harm for a normal run 794 // and makes no harm for a normal run 793 fTheTrueStepLenght = std::min(fTheTrueStep 795 fTheTrueStepLenght = std::min(fTheTrueStepLenght,currentRange); 794 // do the true -> geom transformation 796 // do the true -> geom transformation 795 fTheZPathLenght = fTheTrueStepLenght; 797 fTheZPathLenght = fTheTrueStepLenght; 796 // z = t for very small true-path-length 798 // z = t for very small true-path-length 797 if (fTheTrueStepLenght<tlimitminfix2) { 799 if (fTheTrueStepLenght<tlimitminfix2) { 798 return fTheZPathLenght; 800 return fTheZPathLenght; 799 } 801 } 800 G4double tau = fTheTrueStepLenght/fLambda1 802 G4double tau = fTheTrueStepLenght/fLambda1; 801 if (tau<=tausmall) { 803 if (tau<=tausmall) { 802 fTheZPathLenght = std::min(fTheTrueStepL 804 fTheZPathLenght = std::min(fTheTrueStepLenght, fLambda1); 803 } else if (fTheTrueStepLenght<currentRang 805 } else if (fTheTrueStepLenght<currentRange*dtrl) { 804 if (tau<taulim) fTheZPathLenght = fTheTr 806 if (tau<taulim) fTheZPathLenght = fTheTrueStepLenght*(1.-0.5*tau) ; 805 else fTheZPathLenght = fLambd 807 else fTheZPathLenght = fLambda1*(1.-G4Exp(-tau)); 806 } else if (currentKinEnergy<mass || fTheTr 808 } else if (currentKinEnergy<mass || fTheTrueStepLenght==currentRange) { 807 par1 = 1./currentRange ; // alpha =1 809 par1 = 1./currentRange ; // alpha =1/range_init for Ekin<mass 808 par2 = 1./(par1*fLambda1) ; // 1/(alpha 810 par2 = 1./(par1*fLambda1) ; // 1/(alphaxlambda01) 809 par3 = 1.+par2 ; // 1+1/ 811 par3 = 1.+par2 ; // 1+1/ 810 if (fTheTrueStepLenght<currentRange) { 812 if (fTheTrueStepLenght<currentRange) { 811 fTheZPathLenght = 1./(par1*par3) * (1. 813 fTheZPathLenght = 1./(par1*par3) * (1.-std::pow(1.-par1*fTheTrueStepLenght,par3)); 812 } else { 814 } else { 813 fTheZPathLenght = 1./(par1*par3); 815 fTheZPathLenght = 1./(par1*par3); 814 } 816 } 815 } else { 817 } else { 816 G4double rfin = std::max(currentRange 818 G4double rfin = std::max(currentRange-fTheTrueStepLenght, 0.01*currentRange); 817 G4double T1 = GetEnergy(particle,rf 819 G4double T1 = GetEnergy(particle,rfin,currentCouple); 818 G4double lambda1 = GetTransportMeanFreeP 820 G4double lambda1 = GetTransportMeanFreePathOnly(particle,T1); 819 // 821 // 820 par1 = (fLambda1-lambda1)/(fLambda1*fThe 822 par1 = (fLambda1-lambda1)/(fLambda1*fTheTrueStepLenght); // alpha 821 par2 = 1./(par1*fLambda1); 823 par2 = 1./(par1*fLambda1); 822 par3 = 1.+par2 ; 824 par3 = 1.+par2 ; 823 G4Pow *g4calc = G4Pow::GetInstance(); 825 G4Pow *g4calc = G4Pow::GetInstance(); 824 fTheZPathLenght = 1./(par1*par3) * (1.-g 826 fTheZPathLenght = 1./(par1*par3) * (1.-g4calc->powA(1.-par1*fTheTrueStepLenght,par3)); 825 } 827 } 826 } 828 } 827 fTheZPathLenght = std::min(fTheZPathLenght, 829 fTheZPathLenght = std::min(fTheZPathLenght, fLambda1); 828 // 830 // 829 return fTheZPathLenght; 831 return fTheZPathLenght; 830 } 832 } 831 833 832 834 833 G4double G4GoudsmitSaundersonMscModel::Compute 835 G4double G4GoudsmitSaundersonMscModel::ComputeTrueStepLength(G4double geomStepLength) { 834 // init 836 // init 835 fIsEndedUpOnBoundary = false; 837 fIsEndedUpOnBoundary = false; 836 // step defined other than transportation 838 // step defined other than transportation 837 if (geomStepLength==fTheZPathLenght) { 839 if (geomStepLength==fTheZPathLenght) { 838 return fTheTrueStepLenght; 840 return fTheTrueStepLenght; 839 } 841 } 840 // else :: 842 // else :: 841 // - set the flag that transportation was th 843 // - set the flag that transportation was the winer so DoNothin in DOIT !! 842 // - convert geom -> true by using the mean 844 // - convert geom -> true by using the mean value 843 fIsEndedUpOnBoundary = true; // OR LAST STEP 845 fIsEndedUpOnBoundary = true; // OR LAST STEP 844 fTheZPathLenght = geomStepLength; 846 fTheZPathLenght = geomStepLength; 845 // was a short single scattering step 847 // was a short single scattering step 846 if (fIsEverythingWasDone && !fIsMultipleSace 848 if (fIsEverythingWasDone && !fIsMultipleSacettring) { 847 fTheTrueStepLenght = geomStepLength; 849 fTheTrueStepLenght = geomStepLength; 848 return fTheTrueStepLenght; 850 return fTheTrueStepLenght; 849 } 851 } 850 // t = z for very small step 852 // t = z for very small step 851 if (geomStepLength<tlimitminfix2) { 853 if (geomStepLength<tlimitminfix2) { 852 fTheTrueStepLenght = geomStepLength; 854 fTheTrueStepLenght = geomStepLength; 853 // recalculation 855 // recalculation 854 } else { 856 } else { 855 G4double tlength = geomStepLength; 857 G4double tlength = geomStepLength; 856 if (geomStepLength>fLambda1*tausmall) { 858 if (geomStepLength>fLambda1*tausmall) { 857 if (par1< 0.) { 859 if (par1< 0.) { 858 tlength = -fLambda1*G4Log(1.-geomStepL 860 tlength = -fLambda1*G4Log(1.-geomStepLength/fLambda1) ; 859 } else { 861 } else { 860 if (par1*par3*geomStepLength<1.) { 862 if (par1*par3*geomStepLength<1.) { 861 G4Pow *g4calc = G4Pow::GetInstance() 863 G4Pow *g4calc = G4Pow::GetInstance(); 862 tlength = (1.-g4calc->powA( 1.-par1* 864 tlength = (1.-g4calc->powA( 1.-par1*par3*geomStepLength,1./par3))/par1; 863 } else { 865 } else { 864 tlength = currentRange; 866 tlength = currentRange; 865 } 867 } 866 } 868 } 867 if (tlength<geomStepLength || tlength>fT 869 if (tlength<geomStepLength || tlength>fTheTrueStepLenght) { 868 tlength = geomStepLength; 870 tlength = geomStepLength; 869 } 871 } 870 } 872 } 871 fTheTrueStepLenght = tlength; 873 fTheTrueStepLenght = tlength; 872 } 874 } 873 // 875 // 874 return fTheTrueStepLenght; 876 return fTheTrueStepLenght; 875 } 877 } 876 878 877 G4ThreeVector& 879 G4ThreeVector& 878 G4GoudsmitSaundersonMscModel::SampleScattering 880 G4GoudsmitSaundersonMscModel::SampleScattering(const G4ThreeVector& oldDirection, G4double) { 879 if (steppingAlgorithm==fUseDistanceToBoundar 881 if (steppingAlgorithm==fUseDistanceToBoundary && fIsEverythingWasDone && fIsSingleScattering) { 880 // single scattering was and scattering ha 882 // single scattering was and scattering happend 881 fTheNewDirection.rotateUz(oldDirection); 883 fTheNewDirection.rotateUz(oldDirection); 882 fParticleChange->ProposeMomentumDirection( 884 fParticleChange->ProposeMomentumDirection(fTheNewDirection); 883 return fTheDisplacementVector; 885 return fTheDisplacementVector; 884 } else if (steppingAlgorithm==fUseSafetyPlus 886 } else if (steppingAlgorithm==fUseSafetyPlus) { // error-free stepping 885 if (fIsEndedUpOnBoundary) { // do nothing 887 if (fIsEndedUpOnBoundary) { // do nothing on the boundary 886 return fTheDisplacementVector; 888 return fTheDisplacementVector; 887 } else if (fIsEverythingWasDone) { // evry 889 } else if (fIsEverythingWasDone) { // evrything is done if not optimizations case !!! 888 // check single scattering and see if it 890 // check single scattering and see if it happened 889 if (fIsSingleScattering) { 891 if (fIsSingleScattering) { 890 fTheNewDirection.rotateUz(oldDirection 892 fTheNewDirection.rotateUz(oldDirection); 891 fParticleChange->ProposeMomentumDirect 893 fParticleChange->ProposeMomentumDirection(fTheNewDirection); 892 return fTheDisplacementVector; 894 return fTheDisplacementVector; 893 } 895 } 894 // check if multiple scattering happened 896 // check if multiple scattering happened and do things only if scattering was really happening 895 if (fIsMultipleSacettring && !fIsNoScatt 897 if (fIsMultipleSacettring && !fIsNoScatteringInMSC) { 896 fTheNewDirection.rotateUz(oldDirect 898 fTheNewDirection.rotateUz(oldDirection); 897 fTheDisplacementVector.rotateUz(old 899 fTheDisplacementVector.rotateUz(oldDirection); 898 fParticleChange->ProposeMomentumDir 900 fParticleChange->ProposeMomentumDirection(fTheNewDirection); 899 } 901 } 900 // The only thing that could happen if w 902 // The only thing that could happen if we are here (fUseSafety and fIsEverythingWasDone) 901 // is that single scattering was tried 903 // is that single scattering was tried but did not win so scattering did not happen. 902 // So no displacement and no scattering 904 // So no displacement and no scattering 903 return fTheDisplacementVector; 905 return fTheDisplacementVector; 904 } 906 } 905 // 907 // 906 // The only thing that could still happen 908 // The only thing that could still happen with fUseSafetyPlus is that we are in the 907 // optimization branch: so sample MSC angl 909 // optimization branch: so sample MSC angle here (no displacement) 908 } 910 } 909 //else MSC needs to be done here 911 //else MSC needs to be done here 910 SampleMSC(); 912 SampleMSC(); 911 if (!fIsNoScatteringInMSC) { 913 if (!fIsNoScatteringInMSC) { 912 fTheNewDirection.rotateUz(oldDirection); 914 fTheNewDirection.rotateUz(oldDirection); 913 fParticleChange->ProposeMomentumDirection( 915 fParticleChange->ProposeMomentumDirection(fTheNewDirection); 914 if (!fIsNoDisplace) { 916 if (!fIsNoDisplace) { 915 fTheDisplacementVector.rotateUz(oldDirec 917 fTheDisplacementVector.rotateUz(oldDirection); 916 } 918 } 917 } 919 } 918 // 920 // 919 return fTheDisplacementVector; 921 return fTheDisplacementVector; 920 } 922 } 921 923 922 924 923 void G4GoudsmitSaundersonMscModel::SampleMSC() 925 void G4GoudsmitSaundersonMscModel::SampleMSC() { 924 fIsNoScatteringInMSC = false; 926 fIsNoScatteringInMSC = false; 925 // kinetic energy is assumed to be in Geant4 927 // kinetic energy is assumed to be in Geant4 internal energy unit which is MeV 926 G4double kineticEnergy = currentKinEnergy; 928 G4double kineticEnergy = currentKinEnergy; 927 // 929 // 928 // Energy loss correction: 2 version 930 // Energy loss correction: 2 version 929 G4double eloss = 0.0; 931 G4double eloss = 0.0; 930 // if (fTheTrueStepLenght > currentRange*dtrl 932 // if (fTheTrueStepLenght > currentRange*dtrl) { 931 eloss = kineticEnergy - GetEnergy(particle,c 933 eloss = kineticEnergy - GetEnergy(particle,currentRange-fTheTrueStepLenght,currentCouple); 932 // } else { 934 // } else { 933 // eloss = fTheTrueStepLenght*GetDEDX(parti 935 // eloss = fTheTrueStepLenght*GetDEDX(particle,kineticEnergy,currentCouple); 934 // } 936 // } 935 937 936 G4double tau = 0.;// = kineticEnergy/ele 938 G4double tau = 0.;// = kineticEnergy/electron_mass_c2; // where kinEnergy is the mean kinetic energy 937 G4double tau2 = 0.;// = tau*tau; 939 G4double tau2 = 0.;// = tau*tau; 938 G4double eps0 = 0.;// = eloss/kineticEnerg 940 G4double eps0 = 0.;// = eloss/kineticEnergy0; // energy loss fraction to the begin step energy 939 G4double epsm = 0.;// = eloss/kineticEnerg 941 G4double epsm = 0.;// = eloss/kineticEnergy; // energy loss fraction to the mean step energy 940 942 941 // - init. 943 // - init. 942 G4double efEnergy = kineticEnergy; 944 G4double efEnergy = kineticEnergy; 943 G4double efStep = fTheTrueStepLenght; 945 G4double efStep = fTheTrueStepLenght; 944 946 945 G4double kineticEnergy0 = kineticEnergy; 947 G4double kineticEnergy0 = kineticEnergy; 946 if (gIsUseAccurate) { // - use accurate e 948 if (gIsUseAccurate) { // - use accurate energy loss correction 947 kineticEnergy -= 0.5*eloss; // mean ener 949 kineticEnergy -= 0.5*eloss; // mean energy along the full step 948 // other parameters for energy loss correc 950 // other parameters for energy loss corrections 949 tau = kineticEnergy/electron_m 951 tau = kineticEnergy/electron_mass_c2; // where kinEnergy is the mean kinetic energy 950 tau2 = tau*tau; 952 tau2 = tau*tau; 951 eps0 = eloss/kineticEnergy0; // 953 eps0 = eloss/kineticEnergy0; // energy loss fraction to the begin step energy 952 epsm = eloss/kineticEnergy; // 954 epsm = eloss/kineticEnergy; // energy loss fraction to the mean step energy 953 955 954 efEnergy = kineticEnergy * (1.-epsm 956 efEnergy = kineticEnergy * (1.-epsm*epsm*(6.+10.*tau+5.*tau2)/(24.*tau2+48.*tau+72.)); 955 G4double dum = 0.166666*(4.+tau*(6.+tau 957 G4double dum = 0.166666*(4.+tau*(6.+tau*(7.+tau*(4.+tau))))*(epsm/((tau+1.)*(tau+2.)))*(epsm/((tau+1.)*(tau+2.))); 956 efStep = fTheTrueStepLenght*(1.-d 958 efStep = fTheTrueStepLenght*(1.-dum); 957 } else { // - t 959 } else { // - take only mean energy 958 kineticEnergy -= 0.5*eloss; // mean ener 960 kineticEnergy -= 0.5*eloss; // mean energy along the full step 959 efEnergy = kineticEnergy; 961 efEnergy = kineticEnergy; 960 G4double factor = 1./(1.+0.9784671*kinetic 962 G4double factor = 1./(1.+0.9784671*kineticEnergy); //0.9784671 = 1/(2*m_e) 961 eps0 = eloss/kineticEnergy0; 963 eps0 = eloss/kineticEnergy0; 962 epsm = eps0/(1.-0.5*eps0); 964 epsm = eps0/(1.-0.5*eps0); 963 G4double temp = 0.3*(1 -factor*(1.-0.333 965 G4double temp = 0.3*(1 -factor*(1.-0.333333*factor))*eps0*eps0; 964 efStep = fTheTrueStepLenght*(1.+t 966 efStep = fTheTrueStepLenght*(1.+temp); 965 } 967 } 966 // 968 // 967 // compute elastic mfp, first transport mfp, 969 // compute elastic mfp, first transport mfp, screening parameter, and G1 (with Mott-correction 968 // if it was requested by the user) 970 // if it was requested by the user) 969 fLambda1 = GetTransportMeanFreePath(particle 971 fLambda1 = GetTransportMeanFreePath(particle, efEnergy); 970 // s/lambda_el 972 // s/lambda_el 971 G4double lambdan=0.; 973 G4double lambdan=0.; 972 if (fLambda0>0.0) { 974 if (fLambda0>0.0) { 973 lambdan=efStep/fLambda0; 975 lambdan=efStep/fLambda0; 974 } 976 } 975 if (lambdan<=1.0e-12) { 977 if (lambdan<=1.0e-12) { 976 if (fIsEverythingWasDone) { 978 if (fIsEverythingWasDone) { 977 fTheZPathLenght = fTheTrueStepLenght; 979 fTheZPathLenght = fTheTrueStepLenght; 978 } 980 } 979 fIsNoScatteringInMSC = true; 981 fIsNoScatteringInMSC = true; 980 return; 982 return; 981 } 983 } 982 // first moment: 2.* lambdan *scrA*((1.+scrA 984 // first moment: 2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.); 983 G4double Qn1 = lambdan *fG1; 985 G4double Qn1 = lambdan *fG1; 984 // sample scattering angles 986 // sample scattering angles 985 // new direction, relative to the orriginal 987 // new direction, relative to the orriginal one is in {uss,vss,wss} 986 G4double cosTheta1 = 1.0, sinTheta1 = 0.0, c 988 G4double cosTheta1 = 1.0, sinTheta1 = 0.0, cosTheta2 = 1.0, sinTheta2 = 0.0; 987 G4double cosPhi1 = 1.0, sinPhi1 = 0.0, c 989 G4double cosPhi1 = 1.0, sinPhi1 = 0.0, cosPhi2 = 1.0, sinPhi2 = 0.0; 988 G4double uss = 0.0, vss = 0.0, w 990 G4double uss = 0.0, vss = 0.0, wss = 1.0; 989 G4double x_coord = 0.0, y_coord = 0.0, z 991 G4double x_coord = 0.0, y_coord = 0.0, z_coord = 1.0; 990 G4double u2 = 0.0, v2 = 0.0; 992 G4double u2 = 0.0, v2 = 0.0; 991 // if we are above the upper grid limit with 993 // if we are above the upper grid limit with lambdaxG1=true-length/first-trans-mfp 992 // => izotropic distribution: lambG1_max =7. 994 // => izotropic distribution: lambG1_max =7.992 but set it to 7 993 if (0.5*Qn1 > 7.0){ 995 if (0.5*Qn1 > 7.0){ 994 cosTheta1 = 1.-2.*G4UniformRand(); 996 cosTheta1 = 1.-2.*G4UniformRand(); 995 sinTheta1 = std::sqrt((1.-cosTheta1)*(1.+c 997 sinTheta1 = std::sqrt((1.-cosTheta1)*(1.+cosTheta1)); 996 cosTheta2 = 1.-2.*G4UniformRand(); 998 cosTheta2 = 1.-2.*G4UniformRand(); 997 sinTheta2 = std::sqrt((1.-cosTheta2)*(1.+c 999 sinTheta2 = std::sqrt((1.-cosTheta2)*(1.+cosTheta2)); 998 } else { 1000 } else { 999 // sample 2 scattering cost1, sint1, cost 1001 // sample 2 scattering cost1, sint1, cost2 and sint2 for half path 1000 G4double lekin = G4Log(efEnergy); 1002 G4double lekin = G4Log(efEnergy); 1001 G4double pt2 = efEnergy*(efEnergy+2.0 1003 G4double pt2 = efEnergy*(efEnergy+2.0*CLHEP::electron_mass_c2); 1002 G4double beta2 = pt2/(pt2+CLHEP::electr 1004 G4double beta2 = pt2/(pt2+CLHEP::electron_mass_c2*CLHEP::electron_mass_c2); 1003 // backup GS angular dtr pointer (kineti 1005 // backup GS angular dtr pointer (kinetic energy and delta index in case of Mott-correction) 1004 // if the first was an msc sampling (the 1006 // if the first was an msc sampling (the same will be used if the second is also an msc step) 1005 G4GoudsmitSaundersonTable::GSMSCAngularD 1007 G4GoudsmitSaundersonTable::GSMSCAngularDtr *gsDtr = nullptr; 1006 G4int mcEkinIdx = -1; 1008 G4int mcEkinIdx = -1; 1007 G4int mcDeltIdx = -1; 1009 G4int mcDeltIdx = -1; 1008 G4double transfPar = 0.; 1010 G4double transfPar = 0.; 1009 G4bool isMsc = fGSTable->Sampling(0.5*la 1011 G4bool isMsc = fGSTable->Sampling(0.5*lambdan, 0.5*Qn1, fScrA, cosTheta1, sinTheta1, lekin, beta2, 1010 curren 1012 currentMaterialIndex, &gsDtr, mcEkinIdx, mcDeltIdx, transfPar, 1011 true); 1013 true); 1012 fGSTable->Sampling(0.5*lambdan, 0.5*Qn1, 1014 fGSTable->Sampling(0.5*lambdan, 0.5*Qn1, fScrA, cosTheta2, sinTheta2, lekin, beta2, 1013 currentMaterialIndex, 1015 currentMaterialIndex, &gsDtr, mcEkinIdx, mcDeltIdx, transfPar, !isMsc); 1014 if (cosTheta1+cosTheta2==2.) { // no sca 1016 if (cosTheta1+cosTheta2==2.) { // no scattering happened 1015 if (fIsEverythingWasDone) 1017 if (fIsEverythingWasDone) 1016 fTheZPathLenght = fTheTrueStepLeng 1018 fTheZPathLenght = fTheTrueStepLenght; 1017 fIsNoScatteringInMSC = true; 1019 fIsNoScatteringInMSC = true; 1018 return; 1020 return; 1019 } 1021 } 1020 } 1022 } 1021 // sample 2 azimuthal angles 1023 // sample 2 azimuthal angles 1022 G4double phi1 = CLHEP::twopi*G4UniformRand( 1024 G4double phi1 = CLHEP::twopi*G4UniformRand(); 1023 sinPhi1 = std::sin(phi1); 1025 sinPhi1 = std::sin(phi1); 1024 cosPhi1 = std::cos(phi1); 1026 cosPhi1 = std::cos(phi1); 1025 G4double phi2 = CLHEP::twopi*G4UniformRand( 1027 G4double phi2 = CLHEP::twopi*G4UniformRand(); 1026 sinPhi2 = std::sin(phi2); 1028 sinPhi2 = std::sin(phi2); 1027 cosPhi2 = std::cos(phi2); 1029 cosPhi2 = std::cos(phi2); 1028 1030 1029 // compute final direction realtive to z-di 1031 // compute final direction realtive to z-dir 1030 u2 = sinTheta2*cosPhi2; 1032 u2 = sinTheta2*cosPhi2; 1031 v2 = sinTheta2*sinPhi2; 1033 v2 = sinTheta2*sinPhi2; 1032 G4double u2p = cosTheta1*u2 + sinTheta1*cos 1034 G4double u2p = cosTheta1*u2 + sinTheta1*cosTheta2; 1033 uss = u2p*cosPhi1 - v2*sinPhi1; 1035 uss = u2p*cosPhi1 - v2*sinPhi1; 1034 vss = u2p*sinPhi1 + v2*cosPhi1; 1036 vss = u2p*sinPhi1 + v2*cosPhi1; 1035 wss = cosTheta1*cosTheta2 - sinTheta1*u2; 1037 wss = cosTheta1*cosTheta2 - sinTheta1*u2; 1036 1038 1037 // set new direction (is scattering frame) 1039 // set new direction (is scattering frame) 1038 fTheNewDirection.set(uss,vss,wss); 1040 fTheNewDirection.set(uss,vss,wss); 1039 1041 1040 // set the fTheZPathLenght if we don't samp 1042 // set the fTheZPathLenght if we don't sample displacement and 1041 // we should do everything at the step-limi 1043 // we should do everything at the step-limit-phase before we return 1042 if(fIsNoDisplace && fIsEverythingWasDone) 1044 if(fIsNoDisplace && fIsEverythingWasDone) 1043 fTheZPathLenght = fTheTrueStepLenght; 1045 fTheZPathLenght = fTheTrueStepLenght; 1044 1046 1045 // in optimized-mode if the current-safety 1047 // in optimized-mode if the current-safety > current-range we do not use dispalcement 1046 if(fIsNoDisplace) 1048 if(fIsNoDisplace) 1047 return; 1049 return; 1048 1050 1049 /////////////////////////////////////////// 1051 ////////////////////////////////////////////////////////////////////// 1050 // Compute final position 1052 // Compute final position 1051 Qn1 *= fMCtoQ1; 1053 Qn1 *= fMCtoQ1; 1052 if (gIsUseAccurate) { 1054 if (gIsUseAccurate) { 1053 // correction parameter 1055 // correction parameter 1054 G4double par =1.; 1056 G4double par =1.; 1055 if(Qn1<0.7) par = 1.; 1057 if(Qn1<0.7) par = 1.; 1056 else if (Qn1<7.0) par = -0.031376*Qn1+1. 1058 else if (Qn1<7.0) par = -0.031376*Qn1+1.01356; 1057 else par = 0.79; 1059 else par = 0.79; 1058 1060 1059 // Moments with energy loss correction 1061 // Moments with energy loss correction 1060 // --first the uncorrected (for energy l 1062 // --first the uncorrected (for energy loss) values of gamma, eta, a1=a2=0.5*(1-eta), delta 1061 // gamma = G_2/G_1 based on G2 computed 1063 // gamma = G_2/G_1 based on G2 computed from A by using the Wentzel DCS form of G2 1062 G4double loga = G4Log(1.0+1.0/fScrA); 1064 G4double loga = G4Log(1.0+1.0/fScrA); 1063 G4double gamma = 6.0*fScrA*(1.0 + fScrA 1065 G4double gamma = 6.0*fScrA*(1.0 + fScrA)*(loga*(1.0 + 2.0*fScrA) - 2.0)/fG1; 1064 gamma *= fMCtoG2PerG1; 1066 gamma *= fMCtoG2PerG1; 1065 // sample eta from p(eta)=2*eta i.e. P(e 1067 // sample eta from p(eta)=2*eta i.e. P(eta) = eta_square ;-> P(eta) = rand --> eta = sqrt(rand) 1066 G4double eta = std::sqrt(G4UniformRan 1068 G4double eta = std::sqrt(G4UniformRand()); 1067 G4double eta1 = 0.5*(1 - eta); // use 1069 G4double eta1 = 0.5*(1 - eta); // used more than once 1068 // 0.5 +sqrt(6)/6 = 0.9082483; 1070 // 0.5 +sqrt(6)/6 = 0.9082483; 1069 // 1/(4*sqrt(6)) = 0.1020621; 1071 // 1/(4*sqrt(6)) = 0.1020621; 1070 // (4-sqrt(6)/(24*sqrt(6))) = 0.02637471 1072 // (4-sqrt(6)/(24*sqrt(6))) = 0.026374715 1071 // delta = 0.9082483-(0.1020621-0.026374 1073 // delta = 0.9082483-(0.1020621-0.0263747*gamma)*Qn1 without energy loss cor. 1072 G4double delta = 0.9082483-(0.1020621-0 1074 G4double delta = 0.9082483-(0.1020621-0.0263747*gamma)*Qn1; 1073 1075 1074 // compute alpha1 and alpha2 for energy 1076 // compute alpha1 and alpha2 for energy loss correction 1075 G4double temp1 = 2.0 + tau; 1077 G4double temp1 = 2.0 + tau; 1076 G4double temp = (2.0+tau*temp1)/((tau+1 1078 G4double temp = (2.0+tau*temp1)/((tau+1.0)*temp1); 1077 //Take logarithmic dependence 1079 //Take logarithmic dependence 1078 temp = temp - (tau+1.0)/((tau+2.0)*(loga 1080 temp = temp - (tau+1.0)/((tau+2.0)*(loga*(1.0+fScrA)-1.0)); 1079 temp = temp * epsm; 1081 temp = temp * epsm; 1080 temp1 = 1.0 - temp; 1082 temp1 = 1.0 - temp; 1081 delta = delta + 0.40824829*(eps0*(tau+1. 1083 delta = delta + 0.40824829*(eps0*(tau+1.0)/((tau+2.0)* 1082 (loga*(1.0+fScrA)-1.0)*(loga*(1. 1084 (loga*(1.0+fScrA)-1.0)*(loga*(1.0+2.0*fScrA)-2.0)) - 0.25*temp*temp); 1083 G4double b = eta*delta; 1085 G4double b = eta*delta; 1084 G4double c = eta*(1.0-delta); 1086 G4double c = eta*(1.0-delta); 1085 1087 1086 //Calculate transport direction cosines: 1088 //Calculate transport direction cosines: 1087 // ut,vt,wt is the final position divide 1089 // ut,vt,wt is the final position divided by the true step length 1088 G4double w1v2 = cosTheta1*v2; 1090 G4double w1v2 = cosTheta1*v2; 1089 G4double ut = b*sinTheta1*cosPhi1 + c* 1091 G4double ut = b*sinTheta1*cosPhi1 + c*(cosPhi1*u2 - sinPhi1*w1v2) + eta1*uss*temp1; 1090 G4double vt = b*sinTheta1*sinPhi1 + c* 1092 G4double vt = b*sinTheta1*sinPhi1 + c*(sinPhi1*u2 + cosPhi1*w1v2) + eta1*vss*temp1; 1091 G4double wt = eta1*(1+temp) + b* 1093 G4double wt = eta1*(1+temp) + b*cosTheta1 + c*cosTheta2 + eta1*wss*temp1; 1092 1094 1093 // long step correction 1095 // long step correction 1094 ut *=par; 1096 ut *=par; 1095 vt *=par; 1097 vt *=par; 1096 wt *=par; 1098 wt *=par; 1097 1099 1098 // final position relative to the pre-st 1100 // final position relative to the pre-step point in the scattering frame 1099 // ut = x_f/s so needs to multiply by s 1101 // ut = x_f/s so needs to multiply by s 1100 x_coord = ut*fTheTrueStepLenght; 1102 x_coord = ut*fTheTrueStepLenght; 1101 y_coord = vt*fTheTrueStepLenght; 1103 y_coord = vt*fTheTrueStepLenght; 1102 z_coord = wt*fTheTrueStepLenght; 1104 z_coord = wt*fTheTrueStepLenght; 1103 1105 1104 if(fIsEverythingWasDone){ 1106 if(fIsEverythingWasDone){ 1105 // We sample in the step limit so set 1107 // We sample in the step limit so set fTheZPathLenght = transportDistance 1106 // and lateral displacement (x_coord,y 1108 // and lateral displacement (x_coord,y_coord,z_coord-transportDistance) 1107 //Calculate transport distance 1109 //Calculate transport distance 1108 G4double transportDistance = std::sqr 1110 G4double transportDistance = std::sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord); 1109 // protection 1111 // protection 1110 if(transportDistance>fTheTrueStepLengh 1112 if(transportDistance>fTheTrueStepLenght) 1111 transportDistance = fTheTrueStepLen 1113 transportDistance = fTheTrueStepLenght; 1112 fTheZPathLenght = transportDistance; 1114 fTheZPathLenght = transportDistance; 1113 } 1115 } 1114 // else:: we sample in the DoIt so 1116 // else:: we sample in the DoIt so 1115 // the fTheZPathLenght was already 1117 // the fTheZPathLenght was already set and was taken as transport along zet 1116 fTheDisplacementVector.set(x_coord,y_coo 1118 fTheDisplacementVector.set(x_coord,y_coord,z_coord-fTheZPathLenght); 1117 } else { 1119 } else { 1118 // compute zz = <z>/tPathLength 1120 // compute zz = <z>/tPathLength 1119 // s -> true-path-length 1121 // s -> true-path-length 1120 // z -> geom-path-length:: when PRESTA i 1122 // z -> geom-path-length:: when PRESTA is used z =(def.) <z> 1121 // r -> lateral displacement = s/2 sin(t 1123 // r -> lateral displacement = s/2 sin(theta) => x_f = r cos(phi); y_f = r sin(phi) 1122 G4double zz = 0.0; 1124 G4double zz = 0.0; 1123 if(fIsEverythingWasDone){ 1125 if(fIsEverythingWasDone){ 1124 // We sample in the step limit so set 1126 // We sample in the step limit so set fTheZPathLenght = transportDistance 1125 // and lateral displacement (x_coord, 1127 // and lateral displacement (x_coord,y_coord,z_coord-transportDistance) 1126 if(Qn1<0.1) { // use 3-order Taylor a 1128 if(Qn1<0.1) { // use 3-order Taylor approximation of (1-exp(-x))/x around x=0 1127 zz = 1.0 - Qn1*(0.5 - Qn1*(0.166666 1129 zz = 1.0 - Qn1*(0.5 - Qn1*(0.166666667 - 0.041666667*Qn1)); // 1/6 =0.166..7 ; 1/24=0.041.. 1128 } else { 1130 } else { 1129 zz = (1.-G4Exp(-Qn1))/Qn1; 1131 zz = (1.-G4Exp(-Qn1))/Qn1; 1130 } 1132 } 1131 } else { 1133 } else { 1132 // we sample in the DoIt so 1134 // we sample in the DoIt so 1133 // the fTheZPathLenght was already se 1135 // the fTheZPathLenght was already set and was taken as transport along zet 1134 zz = fTheZPathLenght/fTheTrueStepLeng 1136 zz = fTheZPathLenght/fTheTrueStepLenght; 1135 } 1137 } 1136 1138 1137 G4double rr = (1.-zz*zz)/(1.-wss*wss); / 1139 G4double rr = (1.-zz*zz)/(1.-wss*wss); // s^2 >= <z>^2+r^2 :: where r^2 = s^2/4 sin^2(theta) 1138 if(rr >= 0.25) rr = 0.25; // 1140 if(rr >= 0.25) rr = 0.25; // (1-<z>^2/s^2)/sin^2(theta) >= r^2/(s^2 sin^2(theta)) = 1/4 must hold 1139 G4double rperp = fTheTrueStepLenght*std: 1141 G4double rperp = fTheTrueStepLenght*std::sqrt(rr); // this is r/sint 1140 x_coord = rperp*uss; 1142 x_coord = rperp*uss; 1141 y_coord = rperp*vss; 1143 y_coord = rperp*vss; 1142 z_coord = zz*fTheTrueStepLenght; 1144 z_coord = zz*fTheTrueStepLenght; 1143 1145 1144 if(fIsEverythingWasDone){ 1146 if(fIsEverythingWasDone){ 1145 G4double transportDistance = std::sqrt 1147 G4double transportDistance = std::sqrt(x_coord*x_coord + y_coord*y_coord + z_coord*z_coord); 1146 fTheZPathLenght = transportDistance; 1148 fTheZPathLenght = transportDistance; 1147 } 1149 } 1148 1150 1149 fTheDisplacementVector.set(x_coord,y_coo 1151 fTheDisplacementVector.set(x_coord,y_coord,z_coord- fTheZPathLenght); 1150 } 1152 } 1151 } 1153 } 1152 1154