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 file 29 // GEANT4 Class file 30 // 30 // 31 // 31 // 32 // File name: G4WentzelVIModel 32 // File name: G4WentzelVIModel 33 // 33 // 34 // Author: V.Ivanchenko 34 // Author: V.Ivanchenko 35 // 35 // 36 // Creation date: 09.04.2008 from G4MuMscModel 36 // Creation date: 09.04.2008 from G4MuMscModel 37 // 37 // 38 // Modifications: 38 // Modifications: 39 // 27-05-2010 V.Ivanchenko added G4WentzelOKan 39 // 27-05-2010 V.Ivanchenko added G4WentzelOKandVIxSection class to 40 // compute cross sections and sam 40 // compute cross sections and sample scattering angle 41 // 41 // 42 // 42 // 43 // Class Description: 43 // Class Description: 44 // 44 // 45 // Implementation of the model of multiple sca 45 // Implementation of the model of multiple scattering based on 46 // G.Wentzel, Z. Phys. 40 (1927) 590. 46 // G.Wentzel, Z. Phys. 40 (1927) 590. 47 // H.W.Lewis, Phys Rev 78 (1950) 526. 47 // H.W.Lewis, Phys Rev 78 (1950) 526. 48 // J.M. Fernandez-Varea et al., NIM B73 (1993) 48 // J.M. Fernandez-Varea et al., NIM B73 (1993) 447. 49 // L.Urban, CERN-OPEN-2006-077. 49 // L.Urban, CERN-OPEN-2006-077. 50 50 51 // ------------------------------------------- 51 // ------------------------------------------------------------------- 52 // 52 // 53 53 54 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 55 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 56 56 57 #include "G4WentzelVIModel.hh" 57 #include "G4WentzelVIModel.hh" 58 #include "G4PhysicalConstants.hh" 58 #include "G4PhysicalConstants.hh" 59 #include "G4SystemOfUnits.hh" 59 #include "G4SystemOfUnits.hh" 60 #include "Randomize.hh" 60 #include "Randomize.hh" 61 #include "G4ParticleChangeForMSC.hh" 61 #include "G4ParticleChangeForMSC.hh" 62 #include "G4PhysicsTableHelper.hh" 62 #include "G4PhysicsTableHelper.hh" 63 #include "G4ElementVector.hh" 63 #include "G4ElementVector.hh" 64 #include "G4ProductionCutsTable.hh" 64 #include "G4ProductionCutsTable.hh" 65 #include "G4EmParameters.hh" 65 #include "G4EmParameters.hh" 66 #include "G4Log.hh" 66 #include "G4Log.hh" 67 #include "G4Exp.hh" 67 #include "G4Exp.hh" 68 68 69 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 70 70 71 using namespace std; 71 using namespace std; 72 72 73 const G4double invsqrt12 = 1./std::sqrt(12.); << 74 const G4double numlimit = 0.1; << 75 const G4int minNCollisions = 10; << 76 << 77 G4WentzelVIModel::G4WentzelVIModel(G4bool comb 73 G4WentzelVIModel::G4WentzelVIModel(G4bool comb, const G4String& nam) 78 : G4VMscModel(nam), 74 : G4VMscModel(nam), >> 75 ssFactor(1.05), >> 76 invssFactor(1.0), >> 77 currentCouple(nullptr), >> 78 cosThetaMin(1.0), >> 79 cosThetaMax(-1.0), >> 80 fSecondMoments(nullptr), >> 81 idx2(0), >> 82 numlimit(0.1), 79 singleScatteringMode(false), 83 singleScatteringMode(false), 80 isCombined(comb), 84 isCombined(comb), 81 useSecondMoment(false) 85 useSecondMoment(false) 82 { 86 { 83 tlimitminfix = 1.e-6*CLHEP::mm; << 84 lowEnergyLimit = 1.0*CLHEP::eV; << 85 SetSingleScatteringFactor(1.25); 87 SetSingleScatteringFactor(1.25); 86 wokvi = new G4WentzelOKandVIxSection(isCombi << 88 invsqrt12 = 1./sqrt(12.); >> 89 tlimitminfix = 1.e-6*mm; >> 90 lowEnergyLimit = 1.0*eV; >> 91 particle = nullptr; >> 92 nelments = 5; >> 93 xsecn.resize(nelments); >> 94 prob.resize(nelments); >> 95 wokvi = new G4WentzelOKandVIxSection(isCombined); >> 96 fixedCut = -1.0; >> 97 >> 98 minNCollisions = 10; >> 99 >> 100 preKinEnergy = effKinEnergy = tPathLength = zPathLength = lambdaeff >> 101 = currentRange = xtsec = cosTetMaxNuc = 0.0; >> 102 currentMaterialIndex = 0; >> 103 >> 104 fParticleChange = nullptr; >> 105 currentCuts = nullptr; >> 106 currentMaterial = nullptr; 87 } 107 } 88 108 89 //....oooOO0OOooo........oooOO0OOooo........oo 109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 90 110 91 G4WentzelVIModel::~G4WentzelVIModel() 111 G4WentzelVIModel::~G4WentzelVIModel() 92 { 112 { 93 delete wokvi; 113 delete wokvi; 94 if(IsMaster()) { << 114 if(fSecondMoments && IsMaster()) { 95 delete fSecondMoments; 115 delete fSecondMoments; 96 fSecondMoments = nullptr; 116 fSecondMoments = nullptr; 97 } 117 } 98 } 118 } 99 119 100 //....oooOO0OOooo........oooOO0OOooo........oo 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 101 121 102 void G4WentzelVIModel::Initialise(const G4Part 122 void G4WentzelVIModel::Initialise(const G4ParticleDefinition* p, 103 const G4Data 123 const G4DataVector& cuts) 104 { 124 { 105 // reset parameters 125 // reset parameters 106 SetupParticle(p); 126 SetupParticle(p); 107 InitialiseParameters(p); << 108 currentRange = 0.0; 127 currentRange = 0.0; 109 128 110 if(isCombined) { 129 if(isCombined) { 111 G4double tet = PolarAngleLimit(); 130 G4double tet = PolarAngleLimit(); 112 if(tet <= 0.0) { cosThetaMax = 1 131 if(tet <= 0.0) { cosThetaMax = 1.0; } 113 else if(tet < CLHEP::pi) { cosThetaMax = c 132 else if(tet < CLHEP::pi) { cosThetaMax = cos(tet); } 114 } 133 } 115 //G4cout << "G4WentzelVIModel::Initialise " 134 //G4cout << "G4WentzelVIModel::Initialise " << p->GetParticleName() 116 // << " " << this << " " << wokvi << G4end 135 // << " " << this << " " << wokvi << G4endl; 117 136 118 wokvi->Initialise(p, cosThetaMax); 137 wokvi->Initialise(p, cosThetaMax); 119 /* 138 /* 120 G4cout << "G4WentzelVIModel: " << particle-> 139 G4cout << "G4WentzelVIModel: " << particle->GetParticleName() 121 << " 1-cos(ThetaLimit)= " << 1 - cos 140 << " 1-cos(ThetaLimit)= " << 1 - cosThetaMax 122 << " SingScatFactor= " << ssFactor 141 << " SingScatFactor= " << ssFactor 123 << G4endl; 142 << G4endl; 124 */ 143 */ 125 currentCuts = &cuts; 144 currentCuts = &cuts; 126 145 127 // set values of some data members 146 // set values of some data members 128 fParticleChange = GetParticleChangeForMSC(p) 147 fParticleChange = GetParticleChangeForMSC(p); 129 148 130 // Access to materials << 131 const G4ProductionCutsTable* theCoupleTable << 132 G4ProductionCutsTable::GetProductionCutsTa << 133 G4int numOfCouples = (G4int)theCoupleTable-> << 134 nelments = 0; << 135 for(G4int i=0; i<numOfCouples; ++i) { << 136 G4int nelm = (G4int)theCoupleTable->GetMat << 137 nelments = std::max(nelments, nelm); << 138 } << 139 xsecn.resize(nelments); << 140 prob.resize(nelments); << 141 << 142 // build second moment table only if transpo 149 // build second moment table only if transport table is build 143 G4PhysicsTable* table = GetCrossSectionTable 150 G4PhysicsTable* table = GetCrossSectionTable(); 144 if(useSecondMoment && IsMaster() && nullptr << 151 if(useSecondMoment && IsMaster() && table) { 145 152 146 //G4cout << "### G4WentzelVIModel::Initial 153 //G4cout << "### G4WentzelVIModel::Initialise: build 2nd moment table " 147 // << table << G4endl; 154 // << table << G4endl; 148 fSecondMoments = 155 fSecondMoments = 149 G4PhysicsTableHelper::PreparePhysicsTabl 156 G4PhysicsTableHelper::PreparePhysicsTable(fSecondMoments); >> 157 // Access to materials >> 158 const G4ProductionCutsTable* theCoupleTable = >> 159 G4ProductionCutsTable::GetProductionCutsTable(); >> 160 size_t numOfCouples = theCoupleTable->GetTableSize(); 150 161 151 G4bool splineFlag = true; 162 G4bool splineFlag = true; 152 G4PhysicsVector* aVector = nullptr; 163 G4PhysicsVector* aVector = nullptr; 153 G4PhysicsVector* bVector = nullptr; 164 G4PhysicsVector* bVector = nullptr; 154 G4double emin = std::max(LowEnergyLimit(), 165 G4double emin = std::max(LowEnergyLimit(), LowEnergyActivationLimit()); 155 G4double emax = std::min(HighEnergyLimit() 166 G4double emax = std::min(HighEnergyLimit(), HighEnergyActivationLimit()); 156 if(emin < emax) { 167 if(emin < emax) { 157 std::size_t n = G4EmParameters::Instance << 168 size_t n = G4EmParameters::Instance()->NumberOfBinsPerDecade() 158 *G4lrint(std::log10(emax/emin)); 169 *G4lrint(std::log10(emax/emin)); 159 if(n < 3) { n = 3; } 170 if(n < 3) { n = 3; } 160 171 161 for(G4int i=0; i<numOfCouples; ++i) { << 172 for(size_t i=0; i<numOfCouples; ++i) { 162 173 163 //G4cout<< "i= " << i << " Flag= " << 174 //G4cout<< "i= " << i << " Flag= " << fSecondMoments->GetFlag(i) 164 // << G4endl; 175 // << G4endl; 165 if(fSecondMoments->GetFlag(i)) { 176 if(fSecondMoments->GetFlag(i)) { 166 DefineMaterial(theCoupleTable->GetMa 177 DefineMaterial(theCoupleTable->GetMaterialCutsCouple(i)); 167 178 168 delete (*fSecondMoments)[i]; 179 delete (*fSecondMoments)[i]; 169 if(nullptr == aVector) { << 180 if(!aVector) { 170 aVector = new G4PhysicsLogVector(e << 181 aVector = new G4PhysicsLogVector(emin, emax, n); 171 bVector = aVector; 182 bVector = aVector; 172 } else { 183 } else { 173 bVector = new G4PhysicsVector(*aVe 184 bVector = new G4PhysicsVector(*aVector); 174 } 185 } 175 for(std::size_t j=0; j<n; ++j) { << 186 for(size_t j=0; j<n; ++j) { 176 G4double e = bVector->Energy(j); 187 G4double e = bVector->Energy(j); 177 bVector->PutValue(j, ComputeSecond 188 bVector->PutValue(j, ComputeSecondMoment(p, e)*e*e); 178 } 189 } 179 if(splineFlag) { bVector->FillSecond 190 if(splineFlag) { bVector->FillSecondDerivatives(); } 180 (*fSecondMoments)[i] = bVector; 191 (*fSecondMoments)[i] = bVector; 181 } 192 } 182 } 193 } 183 } 194 } 184 //G4cout << *fSecondMoments << G4endl; 195 //G4cout << *fSecondMoments << G4endl; 185 } 196 } 186 } 197 } 187 198 188 //....oooOO0OOooo........oooOO0OOooo........oo 199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 189 200 190 void G4WentzelVIModel::InitialiseLocal(const G 201 void G4WentzelVIModel::InitialiseLocal(const G4ParticleDefinition*, 191 G4VEmMo 202 G4VEmModel* masterModel) 192 { 203 { 193 fSecondMoments = static_cast<G4WentzelVIMode 204 fSecondMoments = static_cast<G4WentzelVIModel*>(masterModel) 194 ->GetSecondMomentTable(); 205 ->GetSecondMomentTable(); 195 } 206 } 196 207 197 //....oooOO0OOooo........oooOO0OOooo........oo 208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 198 209 199 void G4WentzelVIModel::DefineMaterial(const G4 210 void G4WentzelVIModel::DefineMaterial(const G4MaterialCutsCouple* cup) 200 { 211 { 201 if(cup != currentCouple) { 212 if(cup != currentCouple) { 202 currentCouple = cup; 213 currentCouple = cup; 203 SetCurrentCouple(cup); 214 SetCurrentCouple(cup); 204 currentMaterial = cup->GetMaterial(); 215 currentMaterial = cup->GetMaterial(); 205 currentMaterialIndex = currentCouple->GetI 216 currentMaterialIndex = currentCouple->GetIndex(); 206 } 217 } 207 } 218 } 208 219 209 //....oooOO0OOooo........oooOO0OOooo........oo 220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 210 221 211 G4double G4WentzelVIModel::ComputeCrossSection 222 G4double G4WentzelVIModel::ComputeCrossSectionPerAtom( 212 const G4ParticleD 223 const G4ParticleDefinition* p, 213 G4double kinEnerg 224 G4double kinEnergy, 214 G4double Z, G4dou 225 G4double Z, G4double, 215 G4double cutEnerg 226 G4double cutEnergy, G4double) 216 { 227 { 217 G4double cross = 0.0; 228 G4double cross = 0.0; 218 SetupParticle(p); 229 SetupParticle(p); 219 if(kinEnergy < lowEnergyLimit) { return cros 230 if(kinEnergy < lowEnergyLimit) { return cross; } 220 if(nullptr == CurrentCouple()) { << 231 if(!CurrentCouple()) { 221 G4Exception("G4WentzelVIModel::ComputeCros 232 G4Exception("G4WentzelVIModel::ComputeCrossSectionPerAtom", "em0011", 222 FatalException, " G4MaterialCu 233 FatalException, " G4MaterialCutsCouple is not defined"); 223 return 0.0; 234 return 0.0; 224 } 235 } 225 DefineMaterial(CurrentCouple()); 236 DefineMaterial(CurrentCouple()); 226 cosTetMaxNuc = wokvi->SetupKinematic(kinEner 237 cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial); 227 if(cosTetMaxNuc < 1.0) { 238 if(cosTetMaxNuc < 1.0) { 228 G4double cut = (0.0 < fixedCut) ? fixedCu 239 G4double cut = (0.0 < fixedCut) ? fixedCut : cutEnergy; 229 G4double cost = wokvi->SetupTarget(G4lrint 240 G4double cost = wokvi->SetupTarget(G4lrint(Z), cut); 230 cross = wokvi->ComputeTransportCrossSectio 241 cross = wokvi->ComputeTransportCrossSectionPerAtom(cost); 231 /* 242 /* 232 if(p->GetParticleName() == "e-") 243 if(p->GetParticleName() == "e-") 233 G4cout << "G4WentzelVIModel::CS: Z= " << G 244 G4cout << "G4WentzelVIModel::CS: Z= " << G4int(Z) << " e(MeV)= "<<kinEnergy 234 << " 1-cosN= " << 1 - cosTetMaxNuc 245 << " 1-cosN= " << 1 - cosTetMaxNuc << " cross(bn)= " << cross/barn 235 << " " << particle->GetParticleName 246 << " " << particle->GetParticleName() << G4endl; 236 */ 247 */ 237 } 248 } 238 return cross; 249 return cross; 239 } 250 } 240 251 241 //....oooOO0OOooo........oooOO0OOooo........oo 252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 242 253 243 void G4WentzelVIModel::StartTracking(G4Track* 254 void G4WentzelVIModel::StartTracking(G4Track* track) 244 { 255 { 245 /* 256 /* 246 G4cout << "G4WentzelVIModel::StartTracking " 257 G4cout << "G4WentzelVIModel::StartTracking " << track << " " << this << " " 247 << track->GetParticleDefinition()->GetParti 258 << track->GetParticleDefinition()->GetParticleName() 248 << " workvi: " << wokvi << G4endl; 259 << " workvi: " << wokvi << G4endl; 249 */ 260 */ 250 SetupParticle(track->GetParticleDefinition() 261 SetupParticle(track->GetParticleDefinition()); 251 } 262 } 252 263 253 //....oooOO0OOooo........oooOO0OOooo........oo 264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 254 265 255 G4double G4WentzelVIModel::ComputeTruePathLeng 266 G4double G4WentzelVIModel::ComputeTruePathLengthLimit( 256 const G4Track& tr 267 const G4Track& track, 257 G4double& current 268 G4double& currentMinimalStep) 258 { 269 { 259 G4double tlimit = currentMinimalStep; 270 G4double tlimit = currentMinimalStep; 260 const G4DynamicParticle* dp = track.GetDynam 271 const G4DynamicParticle* dp = track.GetDynamicParticle(); 261 const G4StepPoint* sp = track.GetStep()->Get << 272 G4StepPoint* sp = track.GetStep()->GetPreStepPoint(); 262 G4StepStatus stepStatus = sp->GetStepStatus( 273 G4StepStatus stepStatus = sp->GetStepStatus(); 263 singleScatteringMode = false; 274 singleScatteringMode = false; 264 275 265 //G4cout << "G4WentzelVIModel::ComputeTruePa 276 //G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit stepStatus= " 266 // << stepStatus << " " << track.Ge 277 // << stepStatus << " " << track.GetDefinition()->GetParticleName() 267 // << G4endl; 278 // << G4endl; 268 279 269 // initialisation for each step, lambda may 280 // initialisation for each step, lambda may be computed from scratch 270 preKinEnergy = dp->GetKineticEnergy(); 281 preKinEnergy = dp->GetKineticEnergy(); 271 effKinEnergy = preKinEnergy; 282 effKinEnergy = preKinEnergy; 272 DefineMaterial(track.GetMaterialCutsCouple() 283 DefineMaterial(track.GetMaterialCutsCouple()); 273 const G4double logPreKinEnergy = dp->GetLogK << 284 lambdaeff = GetTransportMeanFreePath(particle,preKinEnergy); 274 lambdaeff = GetTransportMeanFreePath(particl << 285 currentRange = GetRange(particle,preKinEnergy,currentCouple); 275 currentRange = GetRange(particle,preKinEnerg << 276 cosTetMaxNuc = wokvi->SetupKinematic(preKinE 286 cosTetMaxNuc = wokvi->SetupKinematic(preKinEnergy, currentMaterial); 277 287 278 //G4cout << "lambdaeff= " << lambdaeff << " 288 //G4cout << "lambdaeff= " << lambdaeff << " Range= " << currentRange 279 // << " tlimit= " << tlimit << " 1-cost= " < 289 // << " tlimit= " << tlimit << " 1-cost= " << 1 - cosTetMaxNuc << G4endl; 280 290 281 // extra check for abnormal situation 291 // extra check for abnormal situation 282 // this check needed to run MSC with eIoni a 292 // this check needed to run MSC with eIoni and eBrem inactivated 283 if(tlimit > currentRange) { tlimit = current 293 if(tlimit > currentRange) { tlimit = currentRange; } 284 294 285 // stop here if small range particle 295 // stop here if small range particle 286 if(tlimit < tlimitminfix) { 296 if(tlimit < tlimitminfix) { 287 return ConvertTrueToGeom(tlimit, currentMi 297 return ConvertTrueToGeom(tlimit, currentMinimalStep); 288 } 298 } 289 299 290 // pre step 300 // pre step 291 G4double presafety = sp->GetSafety(); 301 G4double presafety = sp->GetSafety(); 292 // far from geometry boundary 302 // far from geometry boundary 293 if(currentRange < presafety) { 303 if(currentRange < presafety) { 294 return ConvertTrueToGeom(tlimit, currentMi 304 return ConvertTrueToGeom(tlimit, currentMinimalStep); 295 } 305 } 296 306 297 // compute presafety again if presafety <= 0 307 // compute presafety again if presafety <= 0 and no boundary 298 // i.e. when it is needed for optimization p 308 // i.e. when it is needed for optimization purposes 299 if(stepStatus != fGeomBoundary && presafety 309 if(stepStatus != fGeomBoundary && presafety < tlimitminfix) { 300 presafety = ComputeSafety(sp->GetPosition( 310 presafety = ComputeSafety(sp->GetPosition(), tlimit); 301 if(currentRange < presafety) { 311 if(currentRange < presafety) { 302 return ConvertTrueToGeom(tlimit, current 312 return ConvertTrueToGeom(tlimit, currentMinimalStep); 303 } 313 } 304 } 314 } 305 /* 315 /* 306 G4cout << "e(MeV)= " << preKinEnergy/MeV 316 G4cout << "e(MeV)= " << preKinEnergy/MeV 307 << " " << particle->GetParticleName( 317 << " " << particle->GetParticleName() 308 << " CurLimit(mm)= " << tlimit/mm <<" 318 << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm 309 << " R(mm)= " <<currentRange/mm 319 << " R(mm)= " <<currentRange/mm 310 << " L0(mm^-1)= " << lambdaeff*mm 320 << " L0(mm^-1)= " << lambdaeff*mm 311 << G4endl; 321 << G4endl; 312 */ 322 */ 313 // natural limit for high energy 323 // natural limit for high energy 314 G4double rlimit = std::max(facrange*currentR 324 G4double rlimit = std::max(facrange*currentRange, 315 (1.0 - cosTetMaxN 325 (1.0 - cosTetMaxNuc)*lambdaeff*invssFactor); 316 326 317 // low-energy e- 327 // low-energy e- 318 if(cosThetaMax > cosTetMaxNuc) { 328 if(cosThetaMax > cosTetMaxNuc) { 319 rlimit = std::min(rlimit, facsafety*presaf 329 rlimit = std::min(rlimit, facsafety*presafety); 320 } 330 } 321 331 322 // cut correction 332 // cut correction 323 G4double rcut = currentCouple->GetProduction 333 G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1); 324 //G4cout << "rcut= " << rcut << " rlimit= " 334 //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= " 325 // << presafety << " 1-cosThetaMax= " <<1-co 335 // << presafety << " 1-cosThetaMax= " <<1-cosThetaMax 326 //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc < 336 //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc << G4endl; 327 if(rcut > rlimit) { rlimit = std::min(rlimit << 337 if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); } 328 338 329 tlimit = std::min(tlimit, rlimit); 339 tlimit = std::min(tlimit, rlimit); 330 tlimit = std::max(tlimit, tlimitminfix); 340 tlimit = std::max(tlimit, tlimitminfix); 331 341 332 // step limit in infinite media 342 // step limit in infinite media 333 tlimit = std::min(tlimit, 50*currentMaterial 343 tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom); 334 344 335 //compute geomlimit and force few steps with 345 //compute geomlimit and force few steps within a volume 336 if (steppingAlgorithm == fUseDistanceToBound 346 if (steppingAlgorithm == fUseDistanceToBoundary 337 && stepStatus == fGeomBoundary) { 347 && stepStatus == fGeomBoundary) { 338 348 339 G4double geomlimit = ComputeGeomLimit(trac 349 G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange); 340 tlimit = std::min(tlimit, geomlimit/facgeo 350 tlimit = std::min(tlimit, geomlimit/facgeom); 341 } 351 } 342 /* 352 /* 343 G4cout << particle->GetParticleName() << " e 353 G4cout << particle->GetParticleName() << " e= " << preKinEnergy 344 << " L0= " << lambdaeff << " R= " << 354 << " L0= " << lambdaeff << " R= " << currentRange 345 << " tlimit= " << tlimit 355 << " tlimit= " << tlimit 346 << " currentMinimalStep= " << curre 356 << " currentMinimalStep= " << currentMinimalStep << G4endl; 347 */ 357 */ 348 return ConvertTrueToGeom(tlimit, currentMini 358 return ConvertTrueToGeom(tlimit, currentMinimalStep); 349 } 359 } 350 360 351 //....oooOO0OOooo........oooOO0OOooo........oo 361 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 352 362 353 G4double G4WentzelVIModel::ComputeGeomPathLeng 363 G4double G4WentzelVIModel::ComputeGeomPathLength(G4double truelength) 354 { 364 { 355 zPathLength = tPathLength = truelength; 365 zPathLength = tPathLength = truelength; 356 366 357 // small step use only single scattering 367 // small step use only single scattering 358 cosThetaMin = 1.0; 368 cosThetaMin = 1.0; 359 ComputeTransportXSectionPerVolume(cosThetaMi 369 ComputeTransportXSectionPerVolume(cosThetaMin); 360 //G4cout << "xtsec= " << xtsec << " Nav= " 370 //G4cout << "xtsec= " << xtsec << " Nav= " 361 // << zPathLength*xtsec << G4endl; 371 // << zPathLength*xtsec << G4endl; 362 if(0.0 >= lambdaeff || G4int(zPathLength*xts 372 if(0.0 >= lambdaeff || G4int(zPathLength*xtsec) < minNCollisions) { 363 singleScatteringMode = true; 373 singleScatteringMode = true; 364 lambdaeff = DBL_MAX; 374 lambdaeff = DBL_MAX; 365 375 366 } else { 376 } else { 367 //G4cout << "ComputeGeomPathLength: tLengt 377 //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength 368 // << " Leff= " << lambdaeff << 378 // << " Leff= " << lambdaeff << G4endl; 369 // small step 379 // small step 370 if(tPathLength < numlimit*lambdaeff) { 380 if(tPathLength < numlimit*lambdaeff) { 371 G4double tau = tPathLength/lambdaeff; 381 G4double tau = tPathLength/lambdaeff; 372 zPathLength *= (1.0 - 0.5*tau + tau*tau/ 382 zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0); 373 383 374 // medium step 384 // medium step 375 } else { 385 } else { 376 G4double e1 = 0.0; 386 G4double e1 = 0.0; 377 if(currentRange > tPathLength) { 387 if(currentRange > tPathLength) { 378 e1 = GetEnergy(particle,currentRange-t 388 e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple); 379 } 389 } 380 effKinEnergy = 0.5*(e1 + preKinEnergy); 390 effKinEnergy = 0.5*(e1 + preKinEnergy); 381 cosTetMaxNuc = wokvi->SetupKinematic(eff 391 cosTetMaxNuc = wokvi->SetupKinematic(effKinEnergy, currentMaterial); 382 lambdaeff = GetTransportMeanFreePath(par 392 lambdaeff = GetTransportMeanFreePath(particle, effKinEnergy); 383 //G4cout << " tLength= "<< tPathLength<< 393 //G4cout << " tLength= "<< tPathLength<< " Leff= " << lambdaeff << G4endl; 384 zPathLength = lambdaeff; 394 zPathLength = lambdaeff; 385 if(tPathLength*numlimit < lambdaeff) { 395 if(tPathLength*numlimit < lambdaeff) { 386 zPathLength *= (1.0 - G4Exp(-tPathLeng 396 zPathLength *= (1.0 - G4Exp(-tPathLength/lambdaeff)); 387 } 397 } 388 } 398 } 389 } 399 } 390 //G4cout << "Comp.geom: zLength= "<<zPathLen 400 //G4cout << "Comp.geom: zLength= "<<zPathLength<<" tLength= " 391 // << tPathLength<< " Leff= " << lam 401 // << tPathLength<< " Leff= " << lambdaeff << G4endl; 392 return zPathLength; 402 return zPathLength; 393 } 403 } 394 404 395 //....oooOO0OOooo........oooOO0OOooo........oo 405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 396 406 397 G4double G4WentzelVIModel::ComputeTrueStepLeng 407 G4double G4WentzelVIModel::ComputeTrueStepLength(G4double geomStepLength) 398 { 408 { 399 // initialisation of single scattering x-sec 409 // initialisation of single scattering x-section 400 /* 410 /* 401 G4cout << "ComputeTrueStepLength: Step= " << 411 G4cout << "ComputeTrueStepLength: Step= " << geomStepLength 402 << " geomL= " << zPathLength 412 << " geomL= " << zPathLength 403 << " Lambda= " << lambdaeff 413 << " Lambda= " << lambdaeff 404 << " 1-cosThetaMaxNuc= " << 1 - cos 414 << " 1-cosThetaMaxNuc= " << 1 - cosTetMaxNuc << G4endl; 405 */ 415 */ 406 if(singleScatteringMode) { 416 if(singleScatteringMode) { 407 zPathLength = tPathLength = geomStepLength 417 zPathLength = tPathLength = geomStepLength; 408 418 409 } else { 419 } else { 410 420 411 // step defined by transportation 421 // step defined by transportation 412 // change both geom and true step lengths 422 // change both geom and true step lengths 413 if(geomStepLength < zPathLength) { 423 if(geomStepLength < zPathLength) { 414 424 415 // single scattering 425 // single scattering 416 if(G4int(geomStepLength*xtsec) < minNCol 426 if(G4int(geomStepLength*xtsec) < minNCollisions) { 417 zPathLength = tPathLength = geomStepLe 427 zPathLength = tPathLength = geomStepLength; 418 lambdaeff = DBL_MAX; 428 lambdaeff = DBL_MAX; 419 singleScatteringMode = true; 429 singleScatteringMode = true; 420 430 421 // multiple scattering 431 // multiple scattering 422 } else { 432 } else { 423 // small step 433 // small step 424 if(geomStepLength < numlimit*lambdaeff 434 if(geomStepLength < numlimit*lambdaeff) { 425 G4double tau = geomStepLength/lambda 435 G4double tau = geomStepLength/lambdaeff; 426 tPathLength = geomStepLength*(1.0 + 436 tPathLength = geomStepLength*(1.0 + 0.5*tau + tau*tau/3.0); 427 437 428 // energy correction for a big step 438 // energy correction for a big step 429 } else { 439 } else { 430 tPathLength *= geomStepLength/zPathL 440 tPathLength *= geomStepLength/zPathLength; 431 G4double e1 = 0.0; 441 G4double e1 = 0.0; 432 if(currentRange > tPathLength) { 442 if(currentRange > tPathLength) { 433 e1 = GetEnergy(particle,currentRan 443 e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple); 434 } 444 } 435 effKinEnergy = 0.5*(e1 + preKinEnerg 445 effKinEnergy = 0.5*(e1 + preKinEnergy); 436 cosTetMaxNuc = wokvi->SetupKinematic 446 cosTetMaxNuc = wokvi->SetupKinematic(effKinEnergy, currentMaterial); 437 lambdaeff = GetTransportMeanFreePath 447 lambdaeff = GetTransportMeanFreePath(particle, effKinEnergy); 438 G4double tau = geomStepLength/lambda 448 G4double tau = geomStepLength/lambdaeff; 439 449 440 if(tau < 0.999999) { tPathLength = - 450 if(tau < 0.999999) { tPathLength = -lambdaeff*G4Log(1.0 - tau); } 441 else { tPathLength = c 451 else { tPathLength = currentRange; } 442 } 452 } 443 zPathLength = geomStepLength; 453 zPathLength = geomStepLength; 444 } 454 } 445 } 455 } 446 } 456 } 447 // check of step length 457 // check of step length 448 // define threshold angle between single and 458 // define threshold angle between single and multiple scattering 449 if(!singleScatteringMode) { 459 if(!singleScatteringMode) { 450 cosThetaMin -= ssFactor*tPathLength/lambda 460 cosThetaMin -= ssFactor*tPathLength/lambdaeff; 451 xtsec = 0.0; 461 xtsec = 0.0; 452 462 453 // recompute transport cross section - do 463 // recompute transport cross section - do not change energy 454 // anymore - cannot be applied for big ste 464 // anymore - cannot be applied for big steps 455 if(cosThetaMin > cosTetMaxNuc) { 465 if(cosThetaMin > cosTetMaxNuc) { 456 // new computation 466 // new computation 457 G4double cross = ComputeTransportXSectio 467 G4double cross = ComputeTransportXSectionPerVolume(cosThetaMin); 458 //G4cout << "%%%% cross= " << cross << " 468 //G4cout << "%%%% cross= " << cross << " xtsec= " << xtsec 459 // << " 1-cosTMin= " << 1.0 - 469 // << " 1-cosTMin= " << 1.0 - cosThetaMin << G4endl; 460 if(cross <= 0.0) { 470 if(cross <= 0.0) { 461 singleScatteringMode = true; 471 singleScatteringMode = true; 462 tPathLength = zPathLength; 472 tPathLength = zPathLength; 463 lambdaeff = DBL_MAX; 473 lambdaeff = DBL_MAX; 464 cosThetaMin = 1.0; 474 cosThetaMin = 1.0; 465 } else if(xtsec > 0.0) { 475 } else if(xtsec > 0.0) { 466 476 467 lambdaeff = 1./cross; 477 lambdaeff = 1./cross; 468 G4double tau = zPathLength*cross; 478 G4double tau = zPathLength*cross; 469 if(tau < numlimit) { 479 if(tau < numlimit) { 470 tPathLength = zPathLength*(1.0 + 0.5 480 tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0); 471 } else if(tau < 0.999999) { 481 } else if(tau < 0.999999) { 472 tPathLength = -lambdaeff*G4Log(1.0 - 482 tPathLength = -lambdaeff*G4Log(1.0 - tau); 473 } else { 483 } else { 474 tPathLength = currentRange; 484 tPathLength = currentRange; 475 } 485 } 476 } 486 } 477 } 487 } 478 } 488 } 479 tPathLength = std::min(tPathLength, currentR 489 tPathLength = std::min(tPathLength, currentRange); 480 /* 490 /* 481 G4cout <<"Comp.true: zLength= "<<zPathLength 491 G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength 482 <<" Leff(mm)= "<<lambdaeff/mm<<" sig0 492 <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl; 483 G4cout << particle->GetParticleName() << " 1 493 G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin 484 << " 1-cosTetMaxNuc= " << 1-cosTetMax 494 << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc 485 << " e(MeV)= " << preKinEnergy/MeV << 495 << " e(MeV)= " << preKinEnergy/MeV << " " 486 << " SSmode= " << singleScatteringMod 496 << " SSmode= " << singleScatteringMode << G4endl; 487 */ 497 */ 488 return tPathLength; 498 return tPathLength; 489 } 499 } 490 500 491 //....oooOO0OOooo........oooOO0OOooo........oo 501 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 492 502 493 G4ThreeVector& 503 G4ThreeVector& 494 G4WentzelVIModel::SampleScattering(const G4Thr 504 G4WentzelVIModel::SampleScattering(const G4ThreeVector& oldDirection, 495 G4double /* 505 G4double /*safety*/) 496 { 506 { 497 fDisplacement.set(0.0,0.0,0.0); 507 fDisplacement.set(0.0,0.0,0.0); 498 //G4cout << "!##! G4WentzelVIModel::SampleSc 508 //G4cout << "!##! G4WentzelVIModel::SampleScattering for " 499 // << particle->GetParticleName() << 509 // << particle->GetParticleName() << G4endl; 500 510 501 // ignore scattering for zero step length an 511 // ignore scattering for zero step length and energy below the limit 502 if(preKinEnergy < lowEnergyLimit || tPathLen 512 if(preKinEnergy < lowEnergyLimit || tPathLength <= 0.0) 503 { return fDisplacement; } 513 { return fDisplacement; } 504 514 505 G4double invlambda = 0.0; 515 G4double invlambda = 0.0; 506 if(lambdaeff < DBL_MAX) { invlambda = 0.5/la 516 if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; } 507 517 508 // use average kinetic energy over the step 518 // use average kinetic energy over the step 509 G4double cut = (*currentCuts)[currentMateria 519 G4double cut = (*currentCuts)[currentMaterialIndex]; 510 if(fixedCut > 0.0) { cut = fixedCut; } 520 if(fixedCut > 0.0) { cut = fixedCut; } 511 /* 521 /* 512 G4cout <<"SampleScat: E0(MeV)= "<< preKinEne 522 G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV 513 << " Leff= " << lambdaeff <<" sig0( 523 << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec 514 << " xmsc= " << tPathLength*invlamb 524 << " xmsc= " << tPathLength*invlambda 515 << " safety= " << safety << G4endl; 525 << " safety= " << safety << G4endl; 516 */ 526 */ 517 // step limit due msc 527 // step limit due msc 518 G4int nMscSteps = 1; 528 G4int nMscSteps = 1; 519 G4double x0 = tPathLength; 529 G4double x0 = tPathLength; 520 G4double z0 = x0*invlambda; 530 G4double z0 = x0*invlambda; >> 531 //G4double zzz = 0.0; 521 G4double prob2 = 0.0; 532 G4double prob2 = 0.0; 522 533 523 CLHEP::HepRandomEngine* rndmEngine = G4Rando 534 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine(); 524 535 525 // large scattering angle case - two step ap 536 // large scattering angle case - two step approach 526 if(!singleScatteringMode) { 537 if(!singleScatteringMode) { >> 538 static const G4double zzmin = 0.05; 527 if(useSecondMoment) { 539 if(useSecondMoment) { 528 G4double z1 = invlambda*invlambda; 540 G4double z1 = invlambda*invlambda; 529 G4double z2 = SecondMoment(particle, cur 541 G4double z2 = SecondMoment(particle, currentCouple, effKinEnergy); 530 prob2 = (z2 - z1)/(1.5*z1 - z2); 542 prob2 = (z2 - z1)/(1.5*z1 - z2); 531 } 543 } 532 // numerical limitation on step length for << 544 // if(z0 > zzmin && safety > tlimitminfix) { 533 if (z0 > 0.05 && z0 < 30.) { << 545 if(z0 > zzmin) { 534 x0 *= 0.5; 546 x0 *= 0.5; 535 z0 *= 0.5; 547 z0 *= 0.5; 536 nMscSteps = 2; 548 nMscSteps = 2; 537 G4double zzz = G4Exp(-1.0/z0); << 549 } >> 550 //if(z0 > zzmin) { zzz = G4Exp(-1.0/z0); } >> 551 G4double zzz = 0.0; >> 552 if(z0 > zzmin) { >> 553 zzz = G4Exp(-1.0/z0); 538 z0 += zzz; 554 z0 += zzz; 539 prob2 *= (1.0 + zzz); << 555 prob2 *= (1 + zzz); 540 } 556 } 541 prob2 /= (1.0 + prob2); << 557 prob2 /= (1 + prob2); 542 } 558 } 543 559 544 // step limit due to single scattering 560 // step limit due to single scattering 545 G4double x1 = 2*tPathLength; 561 G4double x1 = 2*tPathLength; 546 if(0.0 < xtsec) { x1 = -G4Log(rndmEngine->fl 562 if(0.0 < xtsec) { x1 = -G4Log(rndmEngine->flat())/xtsec; } 547 563 548 // no scattering case 564 // no scattering case 549 if(singleScatteringMode && x1 > tPathLength) 565 if(singleScatteringMode && x1 > tPathLength) 550 { return fDisplacement; } 566 { return fDisplacement; } 551 567 552 const G4ElementVector* theElementVector = 568 const G4ElementVector* theElementVector = 553 currentMaterial->GetElementVector(); 569 currentMaterial->GetElementVector(); 554 std::size_t nelm = currentMaterial->GetNumbe << 570 G4int nelm = currentMaterial->GetNumberOfElements(); 555 571 556 // geometry 572 // geometry 557 G4double sint, cost, phi; 573 G4double sint, cost, phi; 558 G4ThreeVector temp(0.0,0.0,1.0); 574 G4ThreeVector temp(0.0,0.0,1.0); 559 575 560 // current position and direction relative t 576 // current position and direction relative to the end point 561 // because of magnetic field geometry is com 577 // because of magnetic field geometry is computed relatively to the 562 // end point of the step 578 // end point of the step 563 G4ThreeVector dir(0.0,0.0,1.0); 579 G4ThreeVector dir(0.0,0.0,1.0); 564 fDisplacement.set(0.0,0.0,-zPathLength); 580 fDisplacement.set(0.0,0.0,-zPathLength); 565 581 566 G4double mscfac = zPathLength/tPathLength; 582 G4double mscfac = zPathLength/tPathLength; 567 583 568 // start a loop 584 // start a loop 569 G4double x2 = x0; 585 G4double x2 = x0; 570 G4double step, z; 586 G4double step, z; 571 G4bool singleScat; 587 G4bool singleScat; 572 /* 588 /* 573 G4cout << "Start of the loop x1(mm)= " << 589 G4cout << "Start of the loop x1(mm)= " << x1 << " x2(mm)= " << x2 574 << " 1-cost1= " << 1 - cosThetaMin << " SS 590 << " 1-cost1= " << 1 - cosThetaMin << " SSmode= " << singleScatteringMode 575 << " xtsec= " << xtsec << " Nst= " 591 << " xtsec= " << xtsec << " Nst= " << nMscSteps << G4endl; 576 */ 592 */ 577 do { 593 do { 578 594 579 //G4cout << "# x1(mm)= "<< x1<< " x2(mm)= 595 //G4cout << "# x1(mm)= "<< x1<< " x2(mm)= "<< x2 << G4endl; 580 // single scattering case 596 // single scattering case 581 if(singleScatteringMode && x1 > x2) { 597 if(singleScatteringMode && x1 > x2) { 582 fDisplacement += x2*mscfac*dir; 598 fDisplacement += x2*mscfac*dir; 583 break; 599 break; 584 } 600 } 585 601 586 // what is next single of multiple? 602 // what is next single of multiple? 587 if(x1 <= x2) { 603 if(x1 <= x2) { 588 step = x1; 604 step = x1; 589 singleScat = true; 605 singleScat = true; 590 } else { 606 } else { 591 step = x2; 607 step = x2; 592 singleScat = false; 608 singleScat = false; 593 } 609 } 594 610 595 //G4cout << "# step(mm)= "<< step<< " sin 611 //G4cout << "# step(mm)= "<< step<< " singlScat= "<< singleScat << G4endl; 596 612 597 // new position 613 // new position 598 fDisplacement += step*mscfac*dir; 614 fDisplacement += step*mscfac*dir; 599 615 600 if(singleScat) { 616 if(singleScat) { 601 617 602 // select element 618 // select element 603 std::size_t i = 0; << 619 G4int i = 0; 604 if(nelm > 1) { 620 if(nelm > 1) { 605 G4double qsec = rndmEngine->flat()*xts 621 G4double qsec = rndmEngine->flat()*xtsec; 606 for (; i<nelm; ++i) { if(xsecn[i] >= q 622 for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } } 607 } 623 } 608 G4double cosTetM = 624 G4double cosTetM = 609 wokvi->SetupTarget((*theElementVector) 625 wokvi->SetupTarget((*theElementVector)[i]->GetZasInt(), cut); 610 //G4cout << "!!! " << cosThetaMin << " 626 //G4cout << "!!! " << cosThetaMin << " " << cosTetM << " " 611 // << prob[i] << G4endl; 627 // << prob[i] << G4endl; 612 temp = wokvi->SampleSingleScattering(cos 628 temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]); 613 629 614 // direction is changed 630 // direction is changed 615 temp.rotateUz(dir); 631 temp.rotateUz(dir); 616 dir = temp; 632 dir = temp; 617 //G4cout << dir << G4endl; 633 //G4cout << dir << G4endl; 618 634 619 // new proposed step length 635 // new proposed step length 620 x2 -= step; 636 x2 -= step; 621 x1 = -G4Log(rndmEngine->flat())/xtsec; 637 x1 = -G4Log(rndmEngine->flat())/xtsec; 622 638 623 // multiple scattering 639 // multiple scattering 624 } else { 640 } else { 625 --nMscSteps; 641 --nMscSteps; 626 x1 -= step; 642 x1 -= step; 627 x2 = x0; 643 x2 = x0; 628 644 629 // sample z in interval 0 - 1 645 // sample z in interval 0 - 1 630 G4bool isFirst = true; 646 G4bool isFirst = true; 631 if(prob2 > 0.0 && rndmEngine->flat() < p 647 if(prob2 > 0.0 && rndmEngine->flat() < prob2) { isFirst = false; } 632 do { 648 do { >> 649 //z = -z0*G4Log(1.0 - (1.0 - zzz)*rndmEngine->flat()); 633 if(isFirst) { z = -G4Log(rndmEngine->f 650 if(isFirst) { z = -G4Log(rndmEngine->flat()); } 634 else { z = G4RandGamma::shoot(r 651 else { z = G4RandGamma::shoot(rndmEngine, 2.0, 2.0); } 635 z *= z0; 652 z *= z0; 636 // Loop checking, 03-Aug-2015, Vladimi 653 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko 637 } while(z > 1.0); 654 } while(z > 1.0); 638 655 639 cost = 1.0 - 2.0*z/*factCM*/; 656 cost = 1.0 - 2.0*z/*factCM*/; 640 if(cost > 1.0) { cost = 1.0; } 657 if(cost > 1.0) { cost = 1.0; } 641 else if(cost < -1.0) { cost =-1.0; } 658 else if(cost < -1.0) { cost =-1.0; } 642 sint = sqrt((1.0 - cost)*(1.0 + cost)); 659 sint = sqrt((1.0 - cost)*(1.0 + cost)); 643 phi = twopi*rndmEngine->flat(); 660 phi = twopi*rndmEngine->flat(); 644 G4double vx1 = sint*std::cos(phi); << 661 G4double vx1 = sint*cos(phi); 645 G4double vy1 = sint*std::sin(phi); << 662 G4double vy1 = sint*sin(phi); 646 663 647 // lateral displacement 664 // lateral displacement 648 if (latDisplasment) { 665 if (latDisplasment) { 649 G4double rms = z0 > 0.0 ? invsqrt12*st << 666 G4double rms = invsqrt12*sqrt(2*z0); 650 G4double r = x0*mscfac; 667 G4double r = x0*mscfac; 651 G4double dx = r*(0.5*vx1 + rms*G4RandG 668 G4double dx = r*(0.5*vx1 + rms*G4RandGauss::shoot(rndmEngine,0.0,1.0)); 652 G4double dy = r*(0.5*vy1 + rms*G4RandG 669 G4double dy = r*(0.5*vy1 + rms*G4RandGauss::shoot(rndmEngine,0.0,1.0)); 653 G4double d = r*r - dx*dx - dy*dy; 670 G4double d = r*r - dx*dx - dy*dy; 654 671 655 // change position 672 // change position 656 if(d >= 0.0) { 673 if(d >= 0.0) { 657 temp.set(dx, dy, std::sqrt(d) - r); << 674 temp.set(dx,dy,sqrt(d) - r); 658 temp.rotateUz(dir); 675 temp.rotateUz(dir); 659 fDisplacement += temp; 676 fDisplacement += temp; 660 } 677 } 661 } 678 } 662 // change direction 679 // change direction 663 temp.set(vx1,vy1,cost); 680 temp.set(vx1,vy1,cost); 664 temp.rotateUz(dir); 681 temp.rotateUz(dir); 665 dir = temp; 682 dir = temp; 666 } 683 } 667 // Loop checking, 03-Aug-2015, Vladimir Iv 684 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko 668 } while (0 < nMscSteps); 685 } while (0 < nMscSteps); 669 686 670 dir.rotateUz(oldDirection); 687 dir.rotateUz(oldDirection); 671 688 672 //G4cout<<"G4WentzelVIModel sampling is done 689 //G4cout<<"G4WentzelVIModel sampling is done 1-cost= "<< 1.-dir.z()<<G4endl; 673 // end of sampling ------------------------- 690 // end of sampling ------------------------------- 674 691 675 fParticleChange->ProposeMomentumDirection(di 692 fParticleChange->ProposeMomentumDirection(dir); 676 693 677 // lateral displacement 694 // lateral displacement 678 fDisplacement.rotateUz(oldDirection); 695 fDisplacement.rotateUz(oldDirection); 679 696 680 /* 697 /* 681 G4cout << " r(mm)= " << fDisplacement 698 G4cout << " r(mm)= " << fDisplacement.mag() 682 << " safety= " << safety 699 << " safety= " << safety 683 << " trueStep(mm)= " << tPathL 700 << " trueStep(mm)= " << tPathLength 684 << " geomStep(mm)= " << zPathL 701 << " geomStep(mm)= " << zPathLength 685 << " x= " << fDisplacement.x() 702 << " x= " << fDisplacement.x() 686 << " y= " << fDisplacement.y() 703 << " y= " << fDisplacement.y() 687 << " z= " << fDisplacement.z() 704 << " z= " << fDisplacement.z() 688 << G4endl; 705 << G4endl; 689 */ 706 */ 690 707 691 //G4cout<< "G4WentzelVIModel::SampleScatteri 708 //G4cout<< "G4WentzelVIModel::SampleScattering end NewDir= " << dir<< G4endl; 692 return fDisplacement; 709 return fDisplacement; 693 } 710 } 694 711 695 //....oooOO0OOooo........oooOO0OOooo........oo 712 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 696 713 697 G4double G4WentzelVIModel::ComputeTransportXSe 714 G4double G4WentzelVIModel::ComputeTransportXSectionPerVolume(G4double cosTheta) 698 { 715 { 699 // prepare recomputation of x-sections 716 // prepare recomputation of x-sections 700 const G4ElementVector* theElementVector = cu 717 const G4ElementVector* theElementVector = currentMaterial->GetElementVector(); 701 const G4double* theAtomNumDensityVector = 718 const G4double* theAtomNumDensityVector = 702 currentMaterial->GetVecNbOfAtomsPerVolume( 719 currentMaterial->GetVecNbOfAtomsPerVolume(); 703 G4int nelm = (G4int)currentMaterial->GetNumb << 720 G4int nelm = currentMaterial->GetNumberOfElements(); 704 if(nelm > nelments) { 721 if(nelm > nelments) { 705 nelments = nelm; 722 nelments = nelm; 706 xsecn.resize(nelm); 723 xsecn.resize(nelm); 707 prob.resize(nelm); 724 prob.resize(nelm); 708 } 725 } 709 726 710 // check consistency 727 // check consistency 711 xtsec = 0.0; 728 xtsec = 0.0; 712 if(cosTetMaxNuc >= cosTheta) { return 0.0; } 729 if(cosTetMaxNuc >= cosTheta) { return 0.0; } 713 730 714 G4double cut = (*currentCuts)[currentMateria 731 G4double cut = (*currentCuts)[currentMaterialIndex]; 715 if(fixedCut > 0.0) { cut = fixedCut; } 732 if(fixedCut > 0.0) { cut = fixedCut; } 716 733 717 // loop over elements 734 // loop over elements 718 G4double xs = 0.0; 735 G4double xs = 0.0; 719 for (G4int i=0; i<nelm; ++i) { 736 for (G4int i=0; i<nelm; ++i) { 720 G4double costm = 737 G4double costm = 721 wokvi->SetupTarget((*theElementVector)[i 738 wokvi->SetupTarget((*theElementVector)[i]->GetZasInt(), cut); 722 G4double density = theAtomNumDensityVector 739 G4double density = theAtomNumDensityVector[i]; 723 740 724 G4double esec = 0.0; 741 G4double esec = 0.0; 725 if(costm < cosTheta) { 742 if(costm < cosTheta) { 726 743 727 // recompute the transport x-section 744 // recompute the transport x-section 728 if(1.0 > cosTheta) { 745 if(1.0 > cosTheta) { 729 xs += density*wokvi->ComputeTransportC 746 xs += density*wokvi->ComputeTransportCrossSectionPerAtom(cosTheta); 730 } 747 } 731 // recompute the total x-section 748 // recompute the total x-section 732 G4double nucsec = wokvi->ComputeNuclearC 749 G4double nucsec = wokvi->ComputeNuclearCrossSection(cosTheta, costm); 733 esec = wokvi->ComputeElectronCrossSectio 750 esec = wokvi->ComputeElectronCrossSection(cosTheta, costm); 734 nucsec += esec; 751 nucsec += esec; 735 if(nucsec > 0.0) { esec /= nucsec; } 752 if(nucsec > 0.0) { esec /= nucsec; } 736 xtsec += nucsec*density; 753 xtsec += nucsec*density; 737 } 754 } 738 xsecn[i] = xtsec; 755 xsecn[i] = xtsec; 739 prob[i] = esec; 756 prob[i] = esec; 740 //G4cout << i << " xs= " << xs << " xtsec 757 //G4cout << i << " xs= " << xs << " xtsec= " << xtsec 741 // << " 1-cosTheta= " << 1-cosTheta 758 // << " 1-cosTheta= " << 1-cosTheta 742 // << " 1-cosTetMaxNuc2= " <<1-c 759 // << " 1-cosTetMaxNuc2= " <<1-cosTetMaxNuc2<< G4endl; 743 } 760 } 744 761 745 //G4cout << "ComputeXS result: xsec(1/mm)= 762 //G4cout << "ComputeXS result: xsec(1/mm)= " << xs 746 // << " txsec(1/mm)= " << xtsec <<G4 763 // << " txsec(1/mm)= " << xtsec <<G4endl; 747 return xs; 764 return xs; 748 } 765 } 749 766 750 //....oooOO0OOooo........oooOO0OOooo........oo 767 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 751 768 752 G4double G4WentzelVIModel::ComputeSecondMoment 769 G4double G4WentzelVIModel::ComputeSecondMoment(const G4ParticleDefinition* p, 753 G4double kinEnergy) 770 G4double kinEnergy) 754 { 771 { 755 G4double xs = 0.0; 772 G4double xs = 0.0; 756 773 757 SetupParticle(p); 774 SetupParticle(p); 758 cosTetMaxNuc = wokvi->SetupKinematic(kinEner 775 cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial); 759 776 760 if(cosTetMaxNuc >= 1.0) { return xs; } 777 if(cosTetMaxNuc >= 1.0) { return xs; } 761 778 762 const G4ElementVector* theElementVector = cu 779 const G4ElementVector* theElementVector = currentMaterial->GetElementVector(); 763 const G4double* theAtomNumDensityVector = 780 const G4double* theAtomNumDensityVector = 764 currentMaterial->GetVecNbOfAtomsPerVolume( 781 currentMaterial->GetVecNbOfAtomsPerVolume(); 765 std::size_t nelm = currentMaterial->GetNumbe << 782 G4int nelm = currentMaterial->GetNumberOfElements(); 766 783 767 G4double cut = (*currentCuts)[currentMateria 784 G4double cut = (*currentCuts)[currentMaterialIndex]; 768 if(fixedCut > 0.0) { cut = fixedCut; } 785 if(fixedCut > 0.0) { cut = fixedCut; } 769 786 770 // loop over elements 787 // loop over elements 771 for (std::size_t i=0; i<nelm; ++i) { << 788 for (G4int i=0; i<nelm; ++i) { 772 G4double costm = 789 G4double costm = 773 wokvi->SetupTarget((*theElementVector)[i 790 wokvi->SetupTarget((*theElementVector)[i]->GetZasInt(), cut); 774 xs += theAtomNumDensityVector[i] 791 xs += theAtomNumDensityVector[i] 775 *wokvi->ComputeSecondTransportMoment(cos 792 *wokvi->ComputeSecondTransportMoment(costm); 776 } 793 } 777 return xs; 794 return xs; 778 } 795 } 779 796 780 //....oooOO0OOooo........oooOO0OOooo........oo 797 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 781 798 782 void G4WentzelVIModel::SetSingleScatteringFact 799 void G4WentzelVIModel::SetSingleScatteringFactor(G4double val) 783 { 800 { 784 if(val > 0.05) { 801 if(val > 0.05) { 785 ssFactor = val; 802 ssFactor = val; 786 invssFactor = 1.0/(val - 0.05); 803 invssFactor = 1.0/(val - 0.05); 787 } 804 } 788 } 805 } 789 806 790 //....oooOO0OOooo........oooOO0OOooo........oo 807 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 791 808