Changeset 1196 for trunk/source/processes/electromagnetic/standard/src
- Timestamp:
- Nov 25, 2009, 5:13:58 PM (15 years ago)
- Location:
- trunk/source/processes/electromagnetic/standard/src
- Files:
-
- 48 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/electromagnetic/standard/src/G4ASTARStopping.cc
r1007 r1196 25 25 // 26 26 // $Id: G4ASTARStopping.cc,v 1.8 2008/11/24 18:28:09 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-0 2$27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 29 29 //--------------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4BetheBlochModel.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4BetheBlochModel.cc,v 1.3 1 2009/04/24 14:24:00vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4BetheBlochModel.cc,v 1.34 2009/11/11 23:22:27 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 234 234 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy(); 235 235 G4double eexc2 = eexc*eexc; 236 G4double cden = material->GetIonisation()->GetCdensity();237 G4double mden = material->GetIonisation()->GetMdensity();238 G4double aden = material->GetIonisation()->GetAdensity();239 G4double x0den = material->GetIonisation()->GetX0density();240 G4double x1den = material->GetIonisation()->GetX1density();236 //G4double cden = material->GetIonisation()->GetCdensity(); 237 //G4double mden = material->GetIonisation()->GetMdensity(); 238 //G4double aden = material->GetIonisation()->GetAdensity(); 239 //G4double x0den = material->GetIonisation()->GetX0density(); 240 //G4double x1den = material->GetIonisation()->GetX1density(); 241 241 242 242 G4double eDensity = material->GetElectronDensity(); … … 252 252 // density correction 253 253 G4double x = log(bg2)/twoln10; 254 if ( x >= x0den ) { 255 dedx -= twoln10*x - cden ; 256 if ( x < x1den ) dedx -= aden*pow((x1den-x),mden) ; 257 } 254 //if ( x >= x0den ) { 255 // dedx -= twoln10*x - cden ; 256 // if ( x < x1den ) dedx -= aden*pow((x1den-x),mden) ; 257 //} 258 dedx -= material->GetIonisation()->DensityCorrection(x); 258 259 259 260 // shell correction … … 283 284 G4double length) 284 285 { 285 const G4ParticleDefinition* p = dp->GetDefinition();286 const G4Material* mat = couple->GetMaterial();287 G4double preKinEnergy = dp->GetKineticEnergy();288 G4double e = preKinEnergy - eloss*0.5;289 if(e < 0.0) e = preKinEnergy*0.5;290 291 286 if(isIon) { 287 const G4ParticleDefinition* p = dp->GetDefinition(); 288 const G4Material* mat = couple->GetMaterial(); 289 G4double preKinEnergy = dp->GetKineticEnergy(); 290 G4double e = preKinEnergy - eloss*0.5; 291 if(e < 0.0) e = preKinEnergy*0.5; 292 292 293 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e); 293 294 GetModelOfFluctuations()->SetParticleAndCharge(p, q2); … … 295 296 eloss += length*corr->IonHighOrderCorrections(p,couple,e); 296 297 } 297 298 /* 298 299 if(nuclearStopping && preKinEnergy*proton_mass_c2/mass < chargeSquare*100.*MeV) { 299 300 … … 307 308 eloss += nloss; 308 309 } 309 /*310 310 311 G4cout << "G4ionIonisation::CorrectionsAlongStep: e= " << preKinEnergy 311 312 << " de= " << eloss << " NIEL= " << nloss 312 313 << " dynQ= " << dp->GetCharge()/eplus << G4endl; 313 */314 314 315 fParticleChange->ProposeNonIonizingEnergyDeposit(nloss); 315 316 } 316 317 */ 317 318 } 318 319 -
trunk/source/processes/electromagnetic/standard/src/G4BetheHeitlerModel.cc
r1055 r1196 25 25 // 26 26 // $Id: G4BetheHeitlerModel.cc,v 1.13 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4BohrFluctuations.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4BohrFluctuations.cc,v 1. 7 2009/02/19 19:17:50vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4BohrFluctuations.cc,v 1.8 2009/09/29 11:33:22 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 63 63 :G4VEmFluctuationModel(nam), 64 64 particle(0), 65 minNumberInteractionsBohr( 10.0),65 minNumberInteractionsBohr(2.0), 66 66 minFraction(0.2), 67 67 xmin(0.2), … … 110 110 } 111 111 siga = sqrt(siga); 112 //G4cout << "siga= " << siga << G4endl;113 112 G4double twomeanLoss = meanLoss + meanLoss; 113 //G4cout << "siga= " << siga << " 2edp= " << twomeanLoss <<G4endl; 114 114 115 115 if(twomeanLoss < siga) { … … 130 130 loss = meanLoss*n/navr; 131 131 } 132 // 132 //G4cout << "loss= " << loss << G4endl; 133 133 134 134 return loss; -
trunk/source/processes/electromagnetic/standard/src/G4BraggIonModel.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4BraggIonModel.cc,v 1.2 4 2009/04/09 18:41:18vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4BraggIonModel.cc,v 1.26 2009/11/10 19:25:47 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 136 136 G4double kineticEnergy) 137 137 { 138 //G4cout << "G4BraggIonModel::GetChargeSquareRatio e= " << kineticEnergy << G4endl; 138 139 // this method is called only for ions 139 140 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy); … … 148 149 G4double kineticEnergy) 149 150 { 151 //G4cout << "G4BraggIonModel::GetParticleCharge e= " << kineticEnergy << G4endl; 150 152 // this method is called only for ions 151 153 return corr->GetParticleCharge(p,mat,kineticEnergy); … … 252 254 G4double& eloss, 253 255 G4double&, 254 G4double length)256 G4double /*length*/) 255 257 { 256 258 // this method is called only for ions … … 261 263 if(e < 0.0) e = preKinEnergy*0.5; 262 264 265 //G4cout << "G4BraggIonModel::CorrectionsAlongStep e= " << e << G4endl; 266 263 267 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e); 264 268 GetModelOfFluctuations()->SetParticleAndCharge(p, q2); 265 269 eloss *= q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor; 266 270 /* 267 271 if(nuclearStopping) { 268 272 … … 276 280 eloss += nloss; 277 281 } 278 /*282 279 283 G4cout << "G4ionIonisation::CorrectionsAlongStep: e= " << preKinEnergy 280 284 << " de= " << eloss << " NIEL= " << nloss 281 285 << " dynQ= " << dp->GetCharge()/eplus << G4endl; 282 */286 283 287 fParticleChange->ProposeNonIonizingEnergyDeposit(nloss); 284 288 } 289 */ 290 //G4cout << "G4BraggIonModel::CorrectionsAlongStep end" << G4endl; 285 291 } 286 292 -
trunk/source/processes/electromagnetic/standard/src/G4BraggModel.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4BraggModel.cc,v 1.2 2 2009/04/09 18:41:18vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4BraggModel.cc,v 1.23 2009/11/10 19:25:47 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 240 240 241 241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 242 242 /* 243 243 void G4BraggModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple, 244 244 const G4DynamicParticle* dp, … … 263 263 eloss += nloss; 264 264 } 265 /*265 266 266 G4cout << "G4ionIonisation::CorrectionsAlongStep: e= " << preKinEnergy 267 267 << " de= " << eloss << " NIEL= " << nloss 268 268 << " dynQ= " << dp->GetCharge()/eplus << G4endl; 269 */269 270 270 fParticleChange->ProposeNonIonizingEnergyDeposit(nloss); 271 271 } 272 272 } 273 */ 273 274 274 275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -
trunk/source/processes/electromagnetic/standard/src/G4ComptonScattering.cc
r1055 r1196 25 25 // 26 26 // $Id: G4ComptonScattering.cc,v 1.31 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // -
trunk/source/processes/electromagnetic/standard/src/G4ComptonScattering52.cc
r1007 r1196 25 25 // 26 26 // $Id: G4ComptonScattering52.cc,v 1.7 2008/10/15 17:53:44 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-0 2$27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // -
trunk/source/processes/electromagnetic/standard/src/G4CoulombScattering.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4CoulombScattering.cc,v 1.2 0 2009/02/20 12:06:37vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4CoulombScattering.cc,v 1.25 2009/10/28 10:14:13 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 52 52 #include "G4CoulombScatteringModel.hh" 53 53 #include "G4eCoulombScatteringModel.hh" 54 //#include "G4hCoulombScatteringModel.hh" 54 55 #include "G4Electron.hh" 56 #include "G4Proton.hh" 55 57 56 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... … … 103 105 aParticle = p; 104 106 G4double mass = p->GetPDGMass(); 107 G4String name = p->GetParticleName(); 105 108 if (mass > GeV || p->GetParticleType() == "nucleus") { 106 109 SetBuildTableFlag(false); 107 110 verboseLevel = 0; 108 111 } else { 109 G4String name = p->GetParticleName();110 112 if(name != "e-" && name != "e+" && 111 113 name != "mu+" && name != "mu-" && name != "pi+" && … … 115 117 G4double emin = MinKinEnergy(); 116 118 G4double emax = MaxKinEnergy(); 119 // G4gCoulombScatteringModel* model = new G4hCoulombScatteringModel(); 120 G4eCoulombScatteringModel* model = new G4eCoulombScatteringModel(); 121 model->SetPolarAngleLimit(PolarAngleLimit()); 122 model->SetLowEnergyLimit(emin); 123 model->SetHighEnergyLimit(emax); 124 AddEmModel(1, model); 125 /* 126 117 127 G4double eth = thEnergy; 118 128 if(mass < MeV) eth = thEnergyElec; … … 131 141 AddEmModel(2, model); 132 142 } 143 */ 133 144 } 134 145 } … … 139 150 { 140 151 G4cout << " " << PolarAngleLimit()/degree 141 << " < Theta(degree) < 180" 152 << " < Theta(degree) < 180" 142 153 << ", Eth(MeV)= "; 154 143 155 if(aParticle->GetPDGMass() < MeV) G4cout << thEnergyElec; 144 156 else G4cout << thEnergy; -
trunk/source/processes/electromagnetic/standard/src/G4CoulombScatteringModel.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4CoulombScatteringModel.cc,v 1. 39 2009/05/10 16:09:29vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4CoulombScatteringModel.cc,v 1.43 2009/10/28 10:14:13 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 39 39 // 40 40 // Modifications: 41 // 41 42 // 01.08.06 V.Ivanchenko extend upper limit of table to TeV and review the 42 43 // logic of building - only elements from G4ElementTable … … 45 46 // 09.10.07 V.Ivanchenko reorganized methods, add cut dependence in scattering off e- 46 47 // 09.06.08 V.Ivanchenko SelectIsotope is moved to the base class 48 // 16.06.09 Consolandi rows 109, 111-112, 183, 185-186 49 // 47 50 // 48 51 // Class Description: … … 84 87 { 85 88 SetupParticle(p); 86 G4double ekin = std::max(lowEnergyLimit, kinEnergy);87 SetupKinematic( ekin, cutEnergy);89 if(kinEnergy < lowEnergyLimit) return 0.0; 90 SetupKinematic(kinEnergy, cutEnergy); 88 91 89 92 // save lab system kinematics … … 103 106 mom2 = momCM*momCM; 104 107 tkin = sqrt(mom2 + m12) - mass; 105 invbeta2 = 1.0 + m12/mom2; 108 109 //invbeta2 = 1.0 + m12/mom2; 110 111 G4double fm = m2/(mass + m2); 112 invbeta2 = 1.0 + m12*fm*fm/mom2; 106 113 107 114 SetupTarget(Z, tkin); … … 127 134 { 128 135 G4double kinEnergy = dp->GetKineticEnergy(); 129 if(kinEnergy < = DBL_MIN) return;136 if(kinEnergy < lowEnergyLimit) return; 130 137 DefineMaterial(couple); 131 138 SetupParticle(dp->GetDefinition()); 132 G4double ekin = std::max(lowEnergyLimit, kinEnergy); 133 SetupKinematic(ekin, cutEnergy); 139 SetupKinematic(kinEnergy, cutEnergy); 134 140 135 141 // Choose nucleus 136 currentElement = SelectRandomAtom(couple,particle,ekin,ecut,tkin); 142 currentElement = SelectRandomAtom(couple,particle, 143 kinEnergy,ecut,kinEnergy); 137 144 138 145 G4double Z = currentElement->GetZ(); … … 160 167 G4double z1 = 1.0 - cost; 161 168 if(z1 < 0.0) return; 162 163 169 G4double sint = sqrt(z1*(1.0 + cost)); 164 170 G4double phi = twopi * G4UniformRand(); … … 175 181 fParticleChange->ProposeMomentumDirection(newDirection); 176 182 177 G4double elab = gam*(eCM + bet*pzCM); 178 ekin = elab - mass; 179 if(ekin < 0.0) ekin = 0.0; 180 fParticleChange->SetProposedKineticEnergy(ekin); 183 // G4double elab = gam*(eCM + bet*pzCM); 184 185 G4double Ecm = sqrt(mass*mass + m2*m2 + 2.0*etot*m2); 186 G4double elab = etot - m2*(ptot/Ecm)*(ptot/Ecm)*(1.-cost) ; 187 188 189 G4double finalT = elab - mass; 190 if(finalT < 0.0) finalT = 0.0; 191 fParticleChange->SetProposedKineticEnergy(finalT); 181 192 182 193 // recoil 183 G4double erec = kinEnergy - ekin; 184 185 if(erec > recoilThreshold) { 194 G4double erec = kinEnergy - finalT; 195 196 G4double tcut = recoilThreshold; 197 if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); } 198 199 if(erec > tcut) { 186 200 G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz); 187 G4double plab = sqrt( ekin*(ekin+ 2.0*mass));201 G4double plab = sqrt(finalT*(finalT + 2.0*mass)); 188 202 G4ThreeVector p2 = (ptot*dir - plab*newDirection).unit(); 189 203 G4DynamicParticle* newdp = new G4DynamicParticle(ion, p2, erec); -
trunk/source/processes/electromagnetic/standard/src/G4GammaConversion.cc
r1055 r1196 25 25 // 26 26 // $Id: G4GammaConversion.cc,v 1.31 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // -
trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonMscModel.cc
r1058 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4GoudsmitSaundersonMscModel.cc,v 1. 7 2009/04/20 19:22:29vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4GoudsmitSaundersonMscModel.cc,v 1.19 2009/11/09 18:40:25 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 44 44 // step into two sub-steps occur only for msc regime 45 45 // 46 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 46 // 12.06.2009 O.Kadri: linear log-log extrapolation of lambda0 & lambda1 between 1 GeV - 100 TeV 47 // adding a theta min limit due to screening effect of the atomic nucleus 48 // 26.08.2009 O.Kadri: Cubic Spline interpolation was replaced with polynomial method 49 // within CalculateIntegrals method 50 // 05.10.2009 O.Kadri: tuning small angle theta distributions 51 // assuming the case of lambdan<1 as single scattering regime 52 // tuning theta sampling for theta below the screening angle 53 // 54 // 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 48 56 //REFERENCES: 49 57 //Ref.1:E. Benedito et al.,"Mixed simulation ... cross-sections", NIMB 174 (2001) pp 91-110; … … 62 70 #include "G4MaterialCutsCouple.hh" 63 71 #include "G4DynamicParticle.hh" 64 #include "G4DataInterpolation.hh"65 72 #include "G4Electron.hh" 66 73 #include "G4Positron.hh" … … 82 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 83 90 G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel(const G4String& nam) 84 : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy( GeV),isInitialized(false)91 : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy(100.*TeV),isInitialized(false) 85 92 { 86 fr=0.02,rangeinit=0.,masslimite=0.6*MeV ;93 fr=0.02,rangeinit=0.,masslimite=0.6*MeV, 87 94 particle=0;tausmall=1.e-16;taulim=1.e-6;tlimit=1.e10*mm; 88 95 tlimitmin=10.e-6*mm;geombig=1.e50*mm;geommin=1.e-3*mm,tgeom=geombig; … … 117 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 118 125 119 G4double G4GoudsmitSaundersonMscModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p, 120 G4double kineticEnergy,G4double Z, G4double, G4double, G4double) 126 G4double 127 G4GoudsmitSaundersonMscModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p, 128 G4double kineticEnergy,G4double Z, G4double, G4double, G4double) 121 129 { 122 130 //Build cross section table : Taken from Ref.7 … … 126 134 if(kinEnergy>highKEnergy)kinEnergy=highKEnergy; 127 135 128 //value0=Lambda0;value1=Lambda1129 136 G4double value0,value1; 130 137 CalculateIntegrals(p,Z,kinEnergy,value0,value1); … … 136 143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 137 144 138 void G4GoudsmitSaundersonMscModel::SampleScattering(const G4DynamicParticle* dynParticle, 139 G4double safety) 145 void 146 G4GoudsmitSaundersonMscModel::SampleScattering(const G4DynamicParticle* dynParticle, 147 G4double safety) 140 148 { 141 149 G4double kineticEnergy = dynParticle->GetKineticEnergy(); … … 144 152 G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2; 145 153 G4double phi1,phi2,cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0; 146 G4double q1,Gamma,Eta,delta,nu,nu0,nu1,nu2 ,nu_interm;154 G4double q1,Gamma,Eta,delta,nu,nu0,nu1,nu2; 147 155 148 156 /////////////////////////////////////////// 149 157 // Effective energy and path-length from Eq. 4.7.15+16 of Ref.4 150 158 G4double eloss = theManager->GetEnergy(particle,tPathLength,currentCouple); 151 if(eloss>0.5*kineticEnergy) return;159 if(eloss>0.5*kineticEnergy)eloss=kineticEnergy-eloss;//exchange possibility between target atomic e- and incident particle 152 160 G4double ee = kineticEnergy - 0.5*eloss; 153 161 G4double ttau = ee/electron_mass_c2; 154 162 G4double ttau2 = ttau*ttau; 155 163 G4double epsilonpp= eloss/ee; 156 G4double temp2 = 0.166666*(4+ttau*(6+ttau*(7+ttau*(4+ttau))))*(epsilonpp/(ttau+1)/(ttau+2))*(epsilonpp/(ttau+1)/(ttau+2));157 164 G4double cst1=epsilonpp*epsilonpp*(6+10*ttau+5*ttau2)/(24*ttau2+48*ttau+72); 158 165 159 166 kineticEnergy *= (1 - cst1); 160 tPathLength *= (1 - temp2);161 167 /////////////////////////////////////////// 162 // additivity rule for mixture xsection calculation168 // additivity rule for mixture and compound xsection calculation 163 169 const G4Material* mat = currentCouple->GetMaterial(); 164 170 G4int nelm = mat->GetNumberOfElements(); … … 166 172 const G4double* theFraction = mat->GetFractionVector(); 167 173 G4double atomPerVolume = mat->GetTotNbOfAtomsPerVolume(); 168 G4double scrA,llambda0,llambda1; 169 scrA=0.0;llambda0 =0.;llambda1=0.; 174 G4double llambda0=0.,llambda1=0.; 170 175 for(G4int i=0;i<nelm;i++) 171 176 { … … 175 180 llambda1 += (theFraction[i]/l1); 176 181 } 177 if(llambda0>DBL_MIN)llambda0 =1./llambda0; 178 if(llambda1>DBL_MIN)llambda1 =1./llambda1; 179 G4double g1=llambda0/llambda1; 182 if(llambda0>DBL_MIN) llambda0 =1./llambda0; 183 if(llambda1>DBL_MIN) llambda1 =1./llambda1; 184 G4double g1=0.0; 185 if(llambda1>DBL_MIN) g1 = llambda0/llambda1; 186 180 187 G4double x1,x0; 181 182 188 x0=g1/2.; 183 189 do 184 190 { 185 x1 = x0-(x0*((1.+x0)* std::log(1.+1./x0)-1.0)-g1/2.)/( (1.+2.*x0)*std::log(1.+1./x0)-2.0);// x1=x0-f(x0)/f'(x0)191 x1 = x0-(x0*((1.+x0)*log(1.+1./x0)-1.0)-g1/2.)/( (1.+2.*x0)*log(1.+1./x0)-2.0);// x1=x0-f(x0)/f'(x0) 186 192 delta = std::abs( x1 - x0 ); 187 193 x0 = x1; // new approximation becomes the old approximation for the next iteration 188 194 } while (delta > 1e-10); 189 scrA = x1;195 G4double scrA = x1; 190 196 191 197 G4double us=0.0,vs=0.0,ws=1.0,x_coord=0.0,y_coord=0.0,z_coord=1.0; 192 198 G4double lambdan=0.; 199 G4bool mscatt=false,noscatt=false; 200 193 201 if(llambda0>0.)lambdan=atomPerVolume*tPathLength/llambda0; 194 if((lambdan<=1.0e-12)||(lambdan>1.0e+5))return; 195 G4bool noscatt=false; 196 G4bool singlescatt=false; 197 G4bool mscatt=false; 202 if((lambdan<=1.0e-12))return; 198 203 199 204 G4double epsilon1=G4UniformRand(); 200 if(epsilon1<(exp(-lambdan)))noscatt=true;// no scattering 201 else if(epsilon1<((1.+lambdan)*exp(-lambdan)))//single scattering 202 {singlescatt=true; 203 ws=G4UniformRand(); 204 ws= 1.-2.*scrA*ws/(1.-ws + scrA); 205 G4double phi0=twopi*G4UniformRand(); 206 us=sqrt(1.-ws*ws)*cos(phi0); 207 vs=sqrt(1.-ws*ws)*sin(phi0); 208 G4double rr=G4UniformRand(); 209 x_coord=(rr*us); 210 y_coord=(rr*vs); 211 z_coord=((1.-rr)+rr*ws); 205 G4double expn = exp(-lambdan); 206 if((epsilon1<expn)||insideskin)// no scattering 207 {noscatt=true;} 208 else if((epsilon1<((1.+lambdan)*expn)||(lambdan<1.))) 209 { 210 mscatt=false; 211 ws=G4UniformRand(); 212 ws= 1.-2.*scrA*ws/(1.-ws + scrA); 213 if(acos(ws)<sqrt(scrA))//small angle approximation for theta less than screening angle 214 {G4int i=0; 215 do{i++; 216 ws=1.+0.5*atomPerVolume*tPathLength*log(G4UniformRand())/llambda1; 217 }while((fabs(ws)>1.)&&(i<20));//i<20 to avoid time consuming during the run 218 if(i==19)ws=cos(sqrt(scrA)); 219 } 220 G4double phi0=twopi*G4UniformRand(); 221 us=sqrt(1.-ws*ws)*cos(phi0); 222 vs=sqrt(1.-ws*ws)*sin(phi0); 223 G4double rr=G4UniformRand(); 224 x_coord=(rr*us); 225 y_coord=(rr*vs); 226 z_coord=((1.-rr)+rr*ws); 212 227 } 213 228 else 214 {mscatt=true; 215 // Ref.2 subsection 4.4 "The best solution found" 216 // Sample first substep scattering angle 217 SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1); 218 phi1 = twopi*G4UniformRand(); 219 cosPhi1 = cos(phi1); 220 sinPhi1 = sin(phi1); 221 222 // Sample second substep scattering angle 223 SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2); 224 phi2 = twopi*G4UniformRand(); 225 cosPhi2 = cos(phi2); 226 sinPhi2 = sin(phi2); 227 228 // Scattering direction 229 us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1; 230 vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1; 231 ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2; 232 } 229 { 230 mscatt=true; 231 // Ref.2 subsection 4.4 "The best solution found" 232 // Sample first substep scattering angle 233 SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1); 234 phi1 = twopi*G4UniformRand(); 235 cosPhi1 = cos(phi1); 236 sinPhi1 = sin(phi1); 237 238 // Sample second substep scattering angle 239 SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2); 240 phi2 = twopi*G4UniformRand(); 241 cosPhi2 = cos(phi2); 242 sinPhi2 = sin(phi2); 243 244 // Scattering direction 245 us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1; 246 vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1; 247 ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2; 248 } 249 233 250 G4ThreeVector oldDirection = dynParticle->GetMomentumDirection(); 234 251 G4ThreeVector newDirection(us,vs,ws); … … 238 255 if((safety > tlimitminfix)&&(latDisplasment)) 239 256 { 240 // Scattering coordinates 257 241 258 if(mscatt) 242 { 243 if(scrA<DBL_MIN)scrA=DBL_MIN; 244 q1 = 2.*scrA*((1. + scrA)*log(1. + 1./scrA) - 1.); 245 if(q1<DBL_MIN)q1=DBL_MIN; 246 Gamma = 6.*scrA*(1. + scrA)*((1. + 2.*scrA)*log(1. + 1./scrA) - 2.)/q1; 247 Eta = atomPerVolume*tPathLength/llambda1; 248 delta = 0.90824829 - Eta*(0.102062073-Gamma*0.026374715); 249 250 nu = G4UniformRand(); 251 nu = std::sqrt(nu); 252 nu0 = (1.0 - nu)/2.; 253 nu1 = nu*delta; 254 nu2 = nu*(1.0-delta); 255 nu_interm = 1.0 - nu0 - nu1 - nu2; 256 x_coord=(nu1*sinTheta1*cosPhi1+nu2*sinTheta2*(cosPhi1*cosPhi2-cosTheta1*sinPhi1*sinPhi2)+nu_interm*us); 257 y_coord=(nu1*sinTheta1*sinPhi1+nu2*sinTheta2*(sinPhi1*cosPhi2+cosTheta1*cosPhi1*sinPhi2)+nu_interm*vs); 258 z_coord=(nu0+nu1*cosTheta1+nu2*cosTheta2+ nu_interm*ws) ; 259 } 260 G4double r=sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord); 261 262 G4double check= 1.- tPathLength/zPathLength; 263 if(check<=0.) return; 264 else if(r>check) {r=check;x_coord /=r;y_coord /=r;z_coord /=r;} 265 266 r *=tPathLength; 267 259 { 260 if(scrA<DBL_MIN)scrA=DBL_MIN; 261 if(llambda1<DBL_MIN)llambda1=DBL_MIN; 262 q1 = 2.*scrA*((1. + scrA)*log(1. + 1./scrA) - 1.); 263 if(q1<DBL_MIN)q1=DBL_MIN; 264 Gamma = 6.*scrA*(1. + scrA)*((1. + 2.*scrA)*log(1. + 1./scrA) - 2.)/q1; 265 Eta = atomPerVolume*tPathLength/llambda1; 266 delta = 0.90824829 - Eta*(0.102062073-Gamma*0.026374715); 267 268 nu = G4UniformRand(); 269 nu = std::sqrt(nu); 270 nu0 = (1.0 - nu)/2.; 271 nu1 = nu*delta; 272 nu2 = nu*(1.0-delta); 273 x_coord=(nu1*sinTheta1*cosPhi1+nu2*sinTheta2*(cosPhi1*cosPhi2-cosTheta1*sinPhi1*sinPhi2)+nu0*us); 274 y_coord=(nu1*sinTheta1*sinPhi1+nu2*sinTheta2*(sinPhi1*cosPhi2+cosTheta1*cosPhi1*sinPhi2)+nu0*vs); 275 z_coord=(nu0+nu1*cosTheta1+nu2*cosTheta2+ nu0*ws) ; 276 } 277 278 // displacement is computed relatively to the end point 279 if(!noscatt)z_coord -= 1.0;//for noscatt zcoord z_coord !=0. 280 G4double rr = sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord); 281 G4double r = rr*zPathLength; 282 /* 283 G4cout << "G4GS::SampleSecondaries: e(MeV)= " << kineticEnergy 284 << " sinTheta= " << sqrt(1.0 - ws*ws) << " r(mm)= " << r 285 << " trueStep(mm)= " << tPathLength 286 << " geomStep(mm)= " << zPathLength 287 << G4endl; 288 */ 268 289 if(r > tlimitminfix) { 269 290 270 G4ThreeVector latDirection = G4ThreeVector(x_coord*tPathLength,y_coord*tPathLength,z_coord*tPathLength);271 272 291 G4ThreeVector latDirection(x_coord/rr,y_coord/rr,z_coord/rr); 292 latDirection.rotateUz(oldDirection); 293 273 294 ComputeDisplacement(fParticleChange, latDirection, r, safety); 274 295 } … … 277 298 278 299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 279 void G4GoudsmitSaundersonMscModel::SampleCosineTheta(G4double lambdan, G4double scrA, 280 G4double &cost, G4double &sint) 300 301 void 302 G4GoudsmitSaundersonMscModel::SampleCosineTheta(G4double lambdan, G4double scrA, 303 G4double &cost, G4double &sint) 281 304 { 282 305 G4double u,Qn1,r1,tet; … … 284 307 Qn1=2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.); 285 308 286 if((lambdan<1.)||(Qn1<0.001))//plural scatt. or small angle scatt. 287 { 288 G4double xi1,lambdai=0.; 289 G4int i=0; 290 do {xi1=G4UniformRand(); 291 lambdai -=std::log(xi1); 292 xi +=2.*scrA*xi1/(1.-xi1 + scrA); 293 i++; 294 }while((lambdai<lambdan)&&(i<30)); 295 296 } 297 else { 298 if(Qn1>0.5)xi=2.*G4UniformRand();//isotropic distribution 299 else{// procedure described by Benedito in Ref.1 300 do{r1=G4UniformRand(); 301 u=GSTable->SampleTheta(lambdan,scrA,G4UniformRand()); 302 xi = 2.*u; 303 tet=acos(1.-xi); 304 }while(tet*r1*r1>sin(tet)); 305 } 306 } 307 309 if (Qn1<0.001)xi=-0.5*Qn1*log(G4UniformRand()); 310 else if(Qn1>0.5)xi=2.*G4UniformRand();//isotropic distribution 311 else 312 { 313 // procedure described by Benedito in Ref.1 314 do{ 315 r1=G4UniformRand(); 316 u=GSTable->SampleTheta(lambdan,scrA,G4UniformRand()); 317 xi = 2.*u; 318 tet=acos(1.-xi); 319 }while(tet*r1*r1>sin(tet)); 320 } 308 321 309 322 if(xi<0.)xi=0.; … … 315 328 316 329 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 317 // Cubic spline log-log interpolation of Lambda0 and Lambda1 318 // Screening parameter calculated according to Eq. 37 of Ref.1 319 void G4GoudsmitSaundersonMscModel::CalculateIntegrals(const G4ParticleDefinition* p,G4double Z, 320 G4double kinEnergy,G4double &Lam0, 321 G4double &Lam1) 330 // Polynomial log-log interpolation of Lambda0 and Lambda1 between 100 eV - 1 GeV 331 // linear log-log extrapolation between 1 GeV - 100 TeV 332 333 void 334 G4GoudsmitSaundersonMscModel::CalculateIntegrals(const G4ParticleDefinition* p,G4double Z, 335 G4double kinEnergy,G4double &Lam0, 336 G4double &Lam1) 322 337 { 323 //////// BEGIN OF: LAMBDA CALCULATION ////////////////////////////////324 G4double TCSEForThisAtom[106],FTCSEForThisAtom[106],TCSPForThisAtom[106],FTCSPForThisAtom[106];325 338 G4double summ00=0.0; 326 339 G4double summ10=0.0; 327 G4double InterpolatedValue=0.0; 328 329 330 G4int iZ = G4int(Z); 331 if(iZ > 103) iZ = 103; 332 for(G4int i=0;i<106;i++) 333 { 334 TCSEForThisAtom[i]=TCSE[iZ-1][i];FTCSEForThisAtom[i]=FTCSE[iZ-1][i]; 335 TCSPForThisAtom[i]=TCSP[iZ-1][i];FTCSPForThisAtom[i]=FTCSP[iZ-1][i]; 336 } 337 340 G4double x1,x2,y1,y2,acoeff,bcoeff; 338 341 G4double kineticE = kinEnergy; 339 342 if(kineticE<lowKEnergy)kineticE=lowKEnergy; 340 343 if(kineticE>highKEnergy)kineticE=highKEnergy; 341 344 kineticE /= eV; 345 G4double logE=log(kineticE); 346 347 G4int iZ = G4int(Z); 348 if(iZ > 103) iZ = 103; 349 350 G4int enerInd=0; 351 for(G4int i=1;i<106;i++) 352 { 353 if((logE>=ener[i])&&(logE<ener[i+1])){enerInd=i;break;} 354 } 355 356 if(p==G4Electron::Electron()) 357 { 358 if(kineticE<=1.0e+9)//Interpolation of the form y=ax²+b 359 { 360 x1=ener[enerInd]; 361 x2=ener[enerInd+1]; 362 y1=TCSE[iZ-1][enerInd]; 363 y2=TCSE[iZ-1][enerInd+1]; 364 acoeff=(y2-y1)/(x2*x2-x1*x1); 365 bcoeff=y2-acoeff*x2*x2; 366 summ00=acoeff*logE*logE+bcoeff; 367 summ00 =exp(summ00); 368 y1=FTCSE[iZ-1][enerInd]; 369 y2=FTCSE[iZ-1][enerInd+1]; 370 acoeff=(y2-y1)/(x2*x2-x1*x1); 371 bcoeff=y2-acoeff*x2*x2; 372 summ10=acoeff*logE*logE+bcoeff; 373 summ10 =exp(summ10); 374 } 375 else //Interpolation of the form y=ax+b 376 { 377 x1=ener[104]; 378 x2=ener[105]; 379 y1=TCSE[iZ-1][104]; 380 y2=TCSE[iZ-1][105]; 381 summ00=(y2-y1)*(logE-x1)/(x2-x1)+y1; 382 summ00 =exp(summ00); 383 y1=FTCSE[iZ-1][104]; 384 y2=FTCSE[iZ-1][105]; 385 summ10=(y2-y1)*(logE-x1)/(x2-x1)+y1; 386 summ10 =exp(summ10); 387 } 388 } 389 if(p==G4Positron::Positron()) 390 { 391 if(kinEnergy<=1.0e+9) 392 { 393 x1=ener[enerInd]; 394 x2=ener[enerInd+1]; 395 y1=TCSP[iZ-1][enerInd]; 396 y2=TCSP[iZ-1][enerInd+1]; 397 acoeff=(y2-y1)/(x2*x2-x1*x1); 398 bcoeff=y2-acoeff*x2*x2; 399 summ00=acoeff*logE*logE+bcoeff; 400 summ00 =exp(summ00); 401 y1=FTCSP[iZ-1][enerInd]; 402 y2=FTCSP[iZ-1][enerInd+1]; 403 acoeff=(y2-y1)/(x2*x2-x1*x1); 404 bcoeff=y2-acoeff*x2*x2; 405 summ10=acoeff*logE*logE+bcoeff; 406 summ10 =exp(summ10); 407 } 408 else 409 { 410 x1=ener[104]; 411 x2=ener[105]; 412 y1=TCSP[iZ-1][104]; 413 y2=TCSP[iZ-1][105]; 414 summ00=(y2-y1)*(logE-x1)/(x2-x1)+y1; 415 summ00 =exp(summ00); 416 y1=FTCSP[iZ-1][104]; 417 y2=FTCSP[iZ-1][105]; 418 summ10=(y2-y1)*(logE-x1)/(x2-x1)+y1; 419 summ10 =exp(summ10); 420 } 421 } 342 422 343 if(p==G4Electron::Electron())344 {345 MyValue= new G4DataInterpolation(ener,TCSEForThisAtom,106,0.0,0);346 InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));347 delete MyValue;348 summ00 = std::exp(InterpolatedValue);349 MyValue= new G4DataInterpolation(ener,FTCSEForThisAtom,106,0.0,0);350 InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));351 delete MyValue;352 summ10 = std::exp(InterpolatedValue);353 }354 if(p==G4Positron::Positron())355 {356 MyValue= new G4DataInterpolation(ener,TCSPForThisAtom,106,0.0,0);357 InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));358 delete MyValue;359 summ00 = std::exp(InterpolatedValue);360 MyValue= new G4DataInterpolation(ener,FTCSPForThisAtom,106,0.0,0);361 InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));362 delete MyValue;363 summ10 = std::exp(InterpolatedValue);364 }365 366 423 summ00 *=barn; 367 424 summ10 *=barn; … … 374 431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 375 432 //t->g->t step transformations taken from Ref.6 376 G4double G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit(const G4Track& track, 377 G4PhysicsTable* theTable, 378 G4double currentMinimalStep) 433 434 G4double 435 G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit(const G4Track& track, 436 G4PhysicsTable* theTable, 437 G4double currentMinimalStep) 379 438 { 380 439 tPathLength = currentMinimalStep; … … 406 465 407 466 G4double presafety = sp->GetSafety(); 467 468 //G4cout << "G4GS::StepLimit tPathLength= " 469 // <<tPathLength<<" safety= " << presafety 470 // << " range= " <<currentRange<< " lambda= "<<lambda1 471 // << " Alg: " << steppingAlgorithm <<G4endl; 408 472 409 473 // far from geometry boundary … … 437 501 else smallstep = 1.; 438 502 503 //define stepmin here (it depends on lambda!) 504 //rough estimation of lambda_elastic/lambda_transport 505 G4double rat = currentKinEnergy/MeV ; 506 rat = 1.e-3/(rat*(10.+rat)) ; 507 //stepmin ~ lambda_elastic 508 stepmin = rat*lambda1; 509 skindepth = skin*stepmin; 510 //define tlimitmin 511 tlimitmin = 10.*stepmin; 512 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix; 513 514 //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin 515 // << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl; 439 516 // constraint from the geometry 440 517 if((geomlimit < geombig) && (geomlimit > geommin)) … … 448 525 tgeom = geombig; 449 526 450 //define stepmin here (it depends on lambda!)451 //rough estimation of lambda_elastic/lambda_transport452 G4double rat = currentKinEnergy/MeV ;453 rat = 1.e-3/(rat*(10.+rat)) ;454 //stepmin ~ lambda_elastic455 stepmin = rat*lambda1;456 skindepth = skin*stepmin;457 458 //define tlimitmin459 tlimitmin = 10.*stepmin;460 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;461 462 527 } 463 528 … … 471 536 472 537 if(tlimit > tgeom) tlimit = tgeom; 538 539 //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit 540 // << " tlimit= " << tlimit << " presafety= " << presafety << G4endl; 473 541 474 542 // shortcut … … 499 567 if(tlimit < stepmin) tlimit = stepmin; 500 568 501 if(tPathLength > tlimit) tPathLength = tlimit 569 if(tPathLength > tlimit) tPathLength = tlimit; 502 570 503 571 } … … 519 587 520 588 if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined)) 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 589 { 590 rangeinit = currentRange; 591 fr = facrange; 592 // 9.1 like stepping for e+/e- only (not for muons,hadrons) 593 if(mass < masslimite) 594 { 595 if(lambda1 > currentRange) 596 rangeinit = lambda1; 597 if(lambda1 > lambdalimit) 598 fr *= 0.75+0.25*lambda1/lambdalimit; 599 } 600 601 //lower limit for tlimit 602 G4double rat = currentKinEnergy/MeV ; 603 rat = 1.e-3/(rat*(10.+rat)) ; 604 tlimitmin = 10.*lambda1*rat; 605 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix; 606 } 539 607 //step limit 540 608 tlimit = fr*rangeinit; … … 561 629 } 562 630 } 631 //G4cout << "tPathLength= " << tPathLength 632 // << " currentMinimalStep= " << currentMinimalStep << G4endl; 563 633 return tPathLength ; 564 634 } 565 635 566 636 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 637 567 638 G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength(G4double) 568 639 { … … 602 673 zmean = 1./(par1*par3) ; 603 674 } else { 604 G4double T1 = theManager->GetEnergy(particle,currentRange-tPathLength,currentCouple); 675 G4double T1 = theManager->GetEnergy(particle,currentRange-tPathLength, 676 currentCouple); 605 677 606 678 lambda11 = GetLambda(T1); … … 614 686 zPathLength = zmean ; 615 687 // sample z 616 if(samplez) 617 { 688 if(samplez) { 689 618 690 const G4double ztmax = 0.99, onethird = 1./3. ; 619 691 G4double zt = zmean/tPathLength ; 620 692 621 if (tPathLength > stepmin && zt < ztmax) 622 { 693 if (tPathLength > stepmin && zt < ztmax) { 694 623 695 G4double u,cz1; 624 if(zt >= onethird) 625 { 696 if(zt >= onethird) { 697 626 698 G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ; 627 699 cz1 = 1.+cz ; … … 629 701 G4double grej ; 630 702 do { 631 u = exp(log(G4UniformRand())/cz1) ; 632 grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ; 633 } while (grej < G4UniformRand()) ; 634 } 635 else 636 { 703 u = exp(log(G4UniformRand())/cz1) ; 704 grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ; 705 } while (grej < G4UniformRand()) ; 706 707 } else { 637 708 cz1 = 1./zt-1.; 638 709 u = 1.-exp(log(G4UniformRand())/cz1) ; … … 642 713 } 643 714 if(zPathLength > lambda1) zPathLength = lambda1; 644 645 return zPathLength; 646 } 647 648 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 649 650 G4double G4GoudsmitSaundersonMscModel::ComputeTrueStepLength(G4double geomStepLength) 715 //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda1 << G4endl; 716 717 return zPathLength; 718 } 719 720 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 721 722 G4double 723 G4GoudsmitSaundersonMscModel::ComputeTrueStepLength(G4double geomStepLength) 651 724 { 652 725 // step defined other than transportation … … 673 746 } 674 747 if(tPathLength < geomStepLength) tPathLength = geomStepLength; 748 //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl; 675 749 676 750 return tPathLength; … … 678 752 679 753 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 754 //Total & first transport x sections for e-/e+ generated from ELSEPA code 755 680 756 void G4GoudsmitSaundersonMscModel::LoadELSEPAXSections() 681 757 { 682 ///////////////////////////////////////683 //Total & first transport x sections of e-/e+ from ELSEPA code684 758 G4String filename = "XSECTIONS.dat"; 685 759 686 char* path = getenv("G4LEDATA"); 687 if (!path) 688 { 689 G4String excep = "G4GoudsmitSaundersonTable: G4LEDATA environment variable not set"; 690 G4Exception(excep); 691 } 692 693 G4String pathString(path); 694 G4String dirFile = pathString + "/msc_GS/" + filename; 695 FILE *infile; 696 infile = fopen(dirFile,"r"); 697 if (infile == 0) 698 { 699 G4String excep = "G4GoudsmitSaunderson - data files: " + dirFile + " not found"; 700 G4Exception(excep); 701 } 702 703 // Read parameters from tables and take logarithms 704 G4float aRead; 705 for(G4int i=0 ; i<106 ;i++){ 706 fscanf(infile,"%f\t",&aRead); 707 if(aRead > 0.0) aRead = std::log(aRead); 708 else aRead = 0.0; 709 ener[i]=aRead; 760 char* path = getenv("G4LEDATA"); 761 if (!path) 762 { 763 G4String excep = "G4GoudsmitSaundersonTable: G4LEDATA environment variable not set properly"; 764 G4Exception(excep); 765 } 766 767 G4String pathString(path); 768 G4String dirFile = pathString + "/msc_GS/" + filename; 769 FILE *infile; 770 infile = fopen(dirFile,"r"); 771 if (infile == 0) 772 { 773 G4String excep = "G4GoudsmitSaunderson - data files: " + dirFile + " not found"; 774 G4Exception(excep); 775 } 776 777 // Read parameters from tables and take logarithms 778 G4float aRead; 779 for(G4int i=0 ; i<106 ;i++){ 780 fscanf(infile,"%f\t",&aRead); 781 if(aRead > 0.0) aRead = log(aRead); 782 else aRead = 0.0; 783 ener[i]=aRead; 784 } 785 for(G4int j=0;j<103;j++){ 786 for(G4int i=0;i<106;i++){ 787 fscanf(infile,"%f\t",&aRead); 788 if(aRead > 0.0) aRead = log(aRead); 789 else aRead = 0.0; 790 TCSE[j][i]=aRead; 710 791 } 711 for(G4int j=0;j<103;j++){ 712 for(G4int i=0;i<106;i++){ 713 fscanf(infile,"%f\t",&aRead); 714 if(aRead > 0.0) aRead = std::log(aRead); 715 else aRead = 0.0; 716 TCSE[j][i]=aRead; 717 } 718 } 719 for(G4int j=0;j<103;j++){ 720 for(G4int i=0;i<106;i++){ 721 fscanf(infile,"%f\t",&aRead); 722 if(aRead > 0.0) aRead = std::log(aRead); 723 else aRead = 0.0; 724 FTCSE[j][i]=aRead; 725 } 726 } 727 for(G4int j=0;j<103;j++){ 728 for(G4int i=0;i<106;i++){ 729 fscanf(infile,"%f\t",&aRead); 730 if(aRead > 0.0) aRead = std::log(aRead); 731 else aRead = 0.0; 732 TCSP[j][i]=aRead; 733 } 734 } 735 for(G4int j=0;j<103;j++){ 736 for(G4int i=0;i<106;i++){ 737 fscanf(infile,"%f\t",&aRead); 738 if(aRead > 0.0) aRead = std::log(aRead); 739 else aRead = 0.0; 740 FTCSP[j][i]=aRead; 741 } 742 } 743 744 fclose(infile); 745 //End loading XSections and Energies 746 747 } 748 749 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 792 } 793 for(G4int j=0;j<103;j++){ 794 for(G4int i=0;i<106;i++){ 795 fscanf(infile,"%f\t",&aRead); 796 if(aRead > 0.0) aRead = log(aRead); 797 else aRead = 0.0; 798 FTCSE[j][i]=aRead; 799 } 800 } 801 for(G4int j=0;j<103;j++){ 802 for(G4int i=0;i<106;i++){ 803 fscanf(infile,"%f\t",&aRead); 804 if(aRead > 0.0) aRead = log(aRead); 805 else aRead = 0.0; 806 TCSP[j][i]=aRead; 807 } 808 } 809 for(G4int j=0;j<103;j++){ 810 for(G4int i=0;i<106;i++){ 811 fscanf(infile,"%f\t",&aRead); 812 if(aRead > 0.0) aRead = log(aRead); 813 else aRead = 0.0; 814 FTCSP[j][i]=aRead; 815 } 816 } 817 818 fclose(infile); 819 820 } 821 822 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -
trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonTable.cc
r1058 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4GoudsmitSaundersonTable.cc,v 1. 1 2009/03/05 18:48:30vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4GoudsmitSaundersonTable.cc,v 1.4 2009/08/28 16:36:52 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 39 39 // Modifications: 40 40 // 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style 41 // 26.08.2009 O.Kadri: avoiding unuseful calculations and optimizing the root finding parameter error's 42 // within SampleTheta method 41 43 // 42 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... … … 167 169 0.996273,0.997261,0.998097,0.998782,0.999315,0.999695,0.999924,1.0}; 168 170 169 G4double G4GoudsmitSaundersonTable::PDF[] = {-1.};170 G4double G4GoudsmitSaundersonTable::CPDF[] = {-1.};171 G4double* G4GoudsmitSaundersonTable::PDF = 0; 172 G4double* G4GoudsmitSaundersonTable::CPDF = 0; 171 173 172 174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 173 175 G4GoudsmitSaundersonTable::G4GoudsmitSaundersonTable() 174 176 { 175 if(PDF [0] < 0.0){177 if(PDF == 0){ 176 178 G4cout << "### G4GoudsmitSaundersonTable loading PDF data" << G4endl; 179 180 PDF = new G4double [76*11*320]; 181 CPDF= new G4double [76*11*320]; 182 177 183 LoadPDFandCPDFdata(); 178 184 } … … 181 187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 182 188 G4GoudsmitSaundersonTable::~G4GoudsmitSaundersonTable() 183 {} 189 { 190 if(PDF) { 191 delete [] PDF; 192 delete [] CPDF; 193 PDF = 0; 194 } 195 } 184 196 185 197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... … … 194 206 /////////////////////////////////////////////////////////////////////////// 195 207 // Find Lambda and Chia2 Index 196 for(G4int k=0;k<75;k++){if((lambda>=LAMBDAN[k])&&(lambda<LAMBDAN[k+1])) KIndex=k;}208 for(G4int k=0;k<75;k++){if((lambda>=LAMBDAN[k])&&(lambda<LAMBDAN[k+1])){KIndex=k;break;}} 197 209 for(G4int j=0;j<11;j++){A[j]=scrA[11*KIndex+j];} 198 for(G4int j=0;j<10;j++){if((Chia2>=A[j])&&(Chia2<A[j+1])) JIndex=j;}210 for(G4int j=0;j<10;j++){if((Chia2>=A[j])&&(Chia2<A[j+1])){JIndex=j;break;}} 199 211 200 212 /////////////////////////////////////////////////////////////////////////// … … 218 230 ThisPDF[i]=Ckj*PDF1[i]+CkjPlus1*PDF2[i]+CkPlus1j*PDF3[i]+CkPlus1jPlus1*PDF4[i]; 219 231 ThisCPDF[i]=Ckj*CPDF1[i]+CkjPlus1*CPDF2[i]+CkPlus1j*CPDF3[i]+CkPlus1jPlus1*CPDF4[i]; 232 // Find u Index using secant method 233 if((i!=0)&&((rndm>=ThisCPDF[i-1])&&(rndm<ThisCPDF[i]))) {IIndex=i-1;break;} 220 234 } 221 235 222 236 /////////////////////////////////////////////////////////////////////////// 223 // Find u Index using secant method224 for(G4int i=0;i<319;i++){if((rndm>=ThisCPDF[i])&&(rndm<ThisCPDF[i+1]))IIndex=i;}225 237 //CPDF^-1(rndm)=x ==> CPDF(x)=rndm; 226 238 aa=uvalues[IIndex]; 227 239 b=uvalues[IIndex+1]; 240 228 241 do{ 229 242 m=0.5*(aa+b); … … 233 246 if(F>0.)b=m; 234 247 else aa=m; 235 } while(sqrt((b-aa)*(b-aa))>1.0e- 10);248 } while(sqrt((b-aa)*(b-aa))>1.0e-6); 236 249 237 250 return m; -
trunk/source/processes/electromagnetic/standard/src/G4HeatedKleinNishinaCompton.cc
r1058 r1196 25 25 // 26 26 // $Id: G4HeatedKleinNishinaCompton.cc,v 1.5 2009/04/12 17:09:57 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4InitXscPAI.cc
r1007 r1196 26 26 // 27 27 // $Id: G4InitXscPAI.cc,v 1.9 2006/06/29 19:53:00 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // -
trunk/source/processes/electromagnetic/standard/src/G4IonFluctuations.cc
r1055 r1196 25 25 // 26 26 // $Id: G4IonFluctuations.cc,v 1.26 2009/03/31 13:24:40 toshito Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4KleinNishinaCompton.cc
r1055 r1196 25 25 // 26 26 // $Id: G4KleinNishinaCompton.cc,v 1.10 2009/05/15 17:12:33 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4MollerBhabhaModel.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4MollerBhabhaModel.cc,v 1.3 4 2009/04/24 17:15:46vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4MollerBhabhaModel.cc,v 1.35 2009/11/09 19:16:13 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 258 258 259 259 //density correction 260 G4double cden = material->GetIonisation()->GetCdensity();261 G4double mden = material->GetIonisation()->GetMdensity();262 G4double aden = material->GetIonisation()->GetAdensity();263 G4double x0den = material->GetIonisation()->GetX0density();264 G4double x1den = material->GetIonisation()->GetX1density();260 //G4double cden = material->GetIonisation()->GetCdensity(); 261 //G4double mden = material->GetIonisation()->GetMdensity(); 262 //G4double aden = material->GetIonisation()->GetAdensity(); 263 //G4double x0den = material->GetIonisation()->GetX0density(); 264 //G4double x1den = material->GetIonisation()->GetX1density(); 265 265 G4double x = log(bg2)/twoln10; 266 266 267 if (x >= x0den) { 268 dedx -= twoln10*x - cden; 269 if (x < x1den) dedx -= aden*pow(x1den-x, mden); 270 } 267 //if (x >= x0den) { 268 // dedx -= twoln10*x - cden; 269 // if (x < x1den) dedx -= aden*pow(x1den-x, mden); 270 //} 271 dedx -= material->GetIonisation()->DensityCorrection(x); 271 272 272 273 // now you can compute the total ionization loss -
trunk/source/processes/electromagnetic/standard/src/G4MscModel71.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4MscModel71.cc,v 1. 7 2009/04/09 18:41:18vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4MscModel71.cc,v 1.8 2009/11/01 13:05:01 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 104 104 G4MscModel71::G4MscModel71(G4double& m_dtrl, G4double& m_NuclCorrPar, 105 105 G4double& m_FactPar, G4double& m_factail, 106 106 G4bool& m_samplez, const G4String& nam) 107 107 : G4VEmModel(nam), 108 108 taubig(8.0), … … 117 117 { 118 118 stepmin = 1.e-6*mm; 119 currentRange = 0. ; 119 currentRange = 0.; 120 G4cout << G4endl; 121 G4cout << "!!! G4MscModel71 class is obsolete and will be removed for the next major Geant4 release !!!" << G4endl; 122 G4cout << "!!! Please use other models (G4UrbanMscModel90, 92, 93) !!!" << G4endl; 123 G4cout << G4endl; 120 124 } 121 125 -
trunk/source/processes/electromagnetic/standard/src/G4MultipleScattering.cc
r1007 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4MultipleScattering.cc,v 1.7 5 2008/10/15 17:53:44vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4MultipleScattering.cc,v 1.77 2009/11/01 13:05:01 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ----------------------------------------------------------------------------- … … 129 129 130 130 #include "G4MultipleScattering.hh" 131 #include "G4UrbanMscModel .hh"131 #include "G4UrbanMscModel92.hh" 132 132 #include "G4MscStepLimitType.hh" 133 133 #include "G4UrbanMscModel.hh" … … 141 141 { 142 142 isInitialized = false; 143 G4cout << G4endl; 144 G4cout << "!!! G4MultipleScattering class is obsolete and will be removed for the next major Geant4 release !!!" << G4endl; 145 G4cout << "!!! Please use G4eMultipleScattering for e+ and e- !!!" << G4endl; 146 G4cout << "!!! Please use G4MuMultipleScattering for mu+ and mu- !!!" << G4endl; 147 G4cout << "!!! Please use G4hMultipleScattering for hadrons and ions !!!" << G4endl; 148 G4cout << G4endl; 143 149 } 144 150 … … 157 163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 158 164 159 void G4MultipleScattering::InitialiseProcess(const G4ParticleDefinition* p)165 void G4MultipleScattering::InitialiseProcess(const G4ParticleDefinition*) 160 166 { 161 // Modification of parameters between runs 162 if(isInitialized) { 163 if (p->GetParticleType() != "nucleus") { 164 mscUrban->SetStepLimitType(StepLimitType()); 165 mscUrban->SetLateralDisplasmentFlag(LateralDisplasmentFlag()); 166 mscUrban->SetSkin(Skin()); 167 mscUrban->SetRangeFactor(RangeFactor()); 168 mscUrban->SetGeomFactor(GeomFactor()); 169 } 170 return; 171 } 172 173 // defaults for ions, which cannot be overwritten 174 if (p->GetParticleType() == "nucleus") { 175 SetStepLimitType(fMinimal); 176 SetLateralDisplasmentFlag(false); 177 SetBuildLambdaTable(false); 178 } 167 if(isInitialized) return; 179 168 180 169 // initialisation of parameters - defaults for particles other 181 170 // than ions can be overwritten by users 182 mscUrban = new G4UrbanMscModel(); 183 mscUrban->SetStepLimitType(StepLimitType()); 184 mscUrban->SetLateralDisplasmentFlag(LateralDisplasmentFlag()); 185 mscUrban->SetSkin(Skin()); 186 mscUrban->SetRangeFactor(RangeFactor()); 187 mscUrban->SetGeomFactor(GeomFactor()); 188 171 mscUrban = new G4UrbanMscModel92(); 189 172 AddEmModel(1,mscUrban); 190 173 isInitialized = true; -
trunk/source/processes/electromagnetic/standard/src/G4MultipleScattering71.cc
r1007 r1196 25 25 // 26 26 // $Id: G4MultipleScattering71.cc,v 1.5 2008/07/16 11:27:41 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-0 2$27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ----------------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PAIModel.cc,v 1. 47 2009/04/09 18:41:18vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4PAIModel.cc,v 1.51 2009/08/12 21:28:50 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 41 41 // 16.08.04 V.Grichine, bug fixed in massRatio for DEDX, CrossSection, SampleSecondary 42 42 // 08.04.05 Major optimisation of internal interfaces (V.Ivantchenko) 43 // 26.07.09 Fixed logic to work with several materials (V.Ivantchenko) 43 44 // 44 45 … … 63 64 #include "G4ParticleDefinition.hh" 64 65 #include "G4ParticleChangeForLoss.hh" 65 #include "G4GeometryTolerance.hh"66 66 67 67 //////////////////////////////////////////////////////////////////////// … … 162 162 163 163 // Prepare initialization 164 165 fPAItransferTable = new G4PhysicsTable(fTotBin);166 fPAIdEdxTable = new G4PhysicsTable(fTotBin);167 168 164 const G4ProductionCutsTable* theCoupleTable = 169 165 G4ProductionCutsTable::GetProductionCutsTable(); … … 182 178 //G4cout << "Reg <" <<curReg->GetName() << "> mat <" 183 179 // << fMaterial->GetName() << "> fCouple= " 184 // << fCutCouple<< G4endl;180 // << fCutCouple<<" " << p->GetParticleName() <<G4endl; 185 181 if( fCutCouple ) { 186 182 fMaterialCutsCoupleVector.push_back(fCutCouple); 183 184 fPAItransferTable = new G4PhysicsTable(fTotBin+1); 185 fPAIdEdxTable = new G4PhysicsTable(fTotBin+1); 187 186 188 187 fDeltaCutInKinEnergy = … … 263 262 G4double tau, Tmax, Tmin, Tkin, deltaLow, gamma, bg2 ; 264 263 265 if(fdEdxVector) delete fdEdxVector;266 264 fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy, 267 265 fHighestKineticEnergy, … … 272 270 deltaLow = 100.*eV; // 0.5*eV ; 273 271 274 for (G4int i = 0 ; i < fTotBin ; i++) //The loop for the kinetic energy272 for (G4int i = 0 ; i <= fTotBin ; i++) //The loop for the kinetic energy 275 273 { 276 274 LowEdgeEnergy = fParticleEnergyVector->GetLowEdgeEnergy(i) ; … … 283 281 Tkin = Tmax ; 284 282 285 // G4cout<<"proton Tkin = "<<LowEdgeEnergy/MeV<<" MeV"286 // <<" Tmax = "<<Tmax/MeV<<" MeV"<<G4endl;287 288 283 if ( Tmax < Tmin + deltaLow ) // low energy safety 289 284 Tkin = Tmin + deltaLow ; 290 285 291 /*292 G4PAIxSection protonPAI( fMatIndex,293 Tkin,294 bg2,295 fSandiaPhotoAbsCof,296 fSandiaIntervalNumber ) ;297 */298 286 fPAIySection.Initialize(fMaterial, Tkin, bg2); 299 287 … … 325 313 326 314 } // end of Tkin loop 327 // theLossTable->insert(fdEdxVector);328 // end of material loop329 // G4cout<<"G4PAIonisation::BuildPAIonisationTable() have been called"<<G4endl ;330 // G4cout<<"G4PAIonisation::BuildLossTable() have been called"<<G4endl ;331 315 } 332 316 … … 339 323 void G4PAIModel::BuildLambdaVector() 340 324 { 341 //G4double kCarTolerance = G4GeometryTolerance::GetInstance()342 // ->GetSurfaceTolerance();343 344 if (fLambdaVector) delete fLambdaVector;345 if (fdNdxCutVector) delete fdNdxCutVector;346 347 325 fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy, 348 326 fHighestKineticEnergy, … … 356 334 <<fDeltaCutInKinEnergy/keV<<" keV"<<G4endl; 357 335 } 358 for (G4int i = 0 ; i < fTotBin ; i++ )336 for (G4int i = 0 ; i <= fTotBin ; i++ ) 359 337 { 360 338 G4double dNdxCut = GetdNdxCut(i,fDeltaCutInKinEnergy) ; 361 G4double lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut ; 362 363 // if (lambda <= 1000*kCarTolerance) lambda = 1000*kCarTolerance ; // Mmm ??? 339 //G4cout << "i= " << i << " x= " << dNdxCut << G4endl; 340 if(dNdxCut < 0.0) dNdxCut = 0.0; 341 // G4double lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut ; 342 G4double lambda = DBL_MAX; 343 if(dNdxCut > 0.0) lambda = 1.0/dNdxCut; 364 344 365 345 fLambdaVector->PutValue(i, lambda) ; … … 395 375 y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1) ; 396 376 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer) ; 377 if(y1 < 0.0) y1 = 0.0; 378 if(y2 < 0.0) y2 = 0.0; 397 379 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl; 398 380 x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ; … … 409 391 } 410 392 // G4cout<<""<<dNdxCut<<G4endl; 393 if(dNdxCut < 0.0) dNdxCut = 0.0; 411 394 return dNdxCut ; 412 395 } … … 420 403 G4int iTransfer; 421 404 G4double x1, x2, y1, y2, dEdxCut; 422 // 405 //G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl; 423 406 // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) 424 407 // <<G4endl; … … 438 421 y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1) ; 439 422 y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer) ; 440 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl; 423 if(y1 < 0.0) y1 = 0.0; 424 if(y2 < 0.0) y2 = 0.0; 425 //G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl; 441 426 x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ; 442 427 x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ; 443 // 428 //G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; 444 429 445 430 if ( y1 == y2 ) dEdxCut = y2 ; … … 451 436 else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; 452 437 } 453 // G4cout<<""<<dEdxCut<<G4endl; 438 //G4cout<<""<<dEdxCut<<G4endl; 439 if(dEdxCut < 0.0) dEdxCut = 0.0; 454 440 return dEdxCut ; 455 441 } … … 463 449 { 464 450 G4int iTkin,iPlace; 465 size_t jMat;466 451 467 452 //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy); … … 474 459 const G4MaterialCutsCouple* matCC = CurrentCouple(); 475 460 476 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat ) 461 size_t jMat = 0; 462 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat ) 477 463 { 478 464 if( matCC == fMaterialCutsCoupleVector[jMat] ) break; 479 465 } 480 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--; 466 //G4cout << "jMat= " << jMat 467 // << " jMax= " << fMaterialCutsCoupleVector.size() 468 // << " matCC: " << matCC; 469 //if(matCC) G4cout << " mat: " << matCC->GetMaterial()->GetName(); 470 // G4cout << G4endl; 471 if(jMat == fMaterialCutsCoupleVector.size()) return 0.0; 481 472 482 473 fPAIdEdxTable = fPAIdEdxBank[jMat]; 483 474 fdEdxVector = fdEdxTable[jMat]; 484 for(iTkin = 0 ; iTkin < fTotBin ; iTkin++)475 for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++) 485 476 { 486 477 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ; 487 478 } 479 //G4cout << "E= " << scaledTkin << " iTkin= " << iTkin 480 // << " " << fdEdxVector->GetVectorLength()<<G4endl; 488 481 iPlace = iTkin - 1; 489 482 if(iPlace < 0) iPlace = 0; 483 else if(iPlace >= fTotBin) iPlace = fTotBin-1; 490 484 G4double dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) ); 491 485 if( dEdx < 0.) dEdx = 0.; … … 502 496 { 503 497 G4int iTkin,iPlace; 504 size_t jMat;505 498 G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy); 506 499 if(tmax <= cutEnergy) return 0.0; … … 511 504 const G4MaterialCutsCouple* matCC = CurrentCouple(); 512 505 513 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat ) 506 size_t jMat = 0; 507 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat ) 514 508 { 515 509 if( matCC == fMaterialCutsCoupleVector[jMat] ) break; 516 510 } 517 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;511 if(jMat == fMaterialCutsCoupleVector.size()) return 0.0; 518 512 519 513 fPAItransferTable = fPAIxscBank[jMat]; 520 514 521 for(iTkin = 0 ; iTkin < fTotBin ; iTkin++)515 for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++) 522 516 { 523 517 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ; … … 525 519 iPlace = iTkin - 1; 526 520 if(iPlace < 0) iPlace = 0; 527 528 // G4cout<<"iPlace = "<<iPlace<<"; tmax = " 521 else if(iPlace >= fTotBin) iPlace = fTotBin-1; 522 523 //G4cout<<"iPlace = "<<iPlace<<"; tmax = " 529 524 // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl; 530 525 cross1 = GetdNdxCut(iPlace,tmax) ; 531 // 526 //G4cout<<"cross1 = "<<cross1<<G4endl; 532 527 cross2 = GetdNdxCut(iPlace,cutEnergy) ; 533 // 528 //G4cout<<"cross2 = "<<cross2<<G4endl; 534 529 cross = (cross2-cross1)*charge2; 535 // 536 if( cross < DBL_MIN) cross = DBL_MIN;530 //G4cout<<"cross = "<<cross<<G4endl; 531 if( cross < DBL_MIN) cross = 0.0; 537 532 // if( cross2 < DBL_MIN) cross2 = DBL_MIN; 538 533 … … 552 547 G4double maxEnergy) 553 548 { 554 size_t jMat ;555 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size(); ++jMat )549 size_t jMat = 0; 550 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat ) 556 551 { 557 552 if( matCC == fMaterialCutsCoupleVector[jMat] ) break; 558 553 } 559 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;554 if(jMat == fMaterialCutsCoupleVector.size()) return; 560 555 561 556 fPAItransferTable = fPAIxscBank[jMat]; … … 636 631 G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W ; 637 632 638 for(iTkin=0;iTkin< fTotBin;iTkin++)633 for(iTkin=0;iTkin<=fTotBin;iTkin++) 639 634 { 640 635 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ; … … 643 638 // G4cout<<"from search, iPlace = "<<iPlace<<G4endl ; 644 639 if(iPlace < 0) iPlace = 0; 640 else if(iPlace >= fTotBin) iPlace = fTotBin-1; 645 641 dNdxCut1 = (*fdNdxCutVector)(iPlace) ; 646 642 // G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ; … … 758 754 G4double&) 759 755 { 760 size_t jMat ;761 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size(); ++jMat )756 size_t jMat = 0; 757 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat ) 762 758 { 763 759 if( material == fMaterialCutsCoupleVector[jMat]->GetMaterial() ) break; 764 760 } 765 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;761 if(jMat == fMaterialCutsCoupleVector.size()) return 0.0; 766 762 767 763 fPAItransferTable = fPAIxscBank[jMat]; 768 764 fdNdxCutVector = fdNdxCutTable[jMat]; 769 765 770 G4int iTkin, iTransfer, iPlace 766 G4int iTkin, iTransfer, iPlace; 771 767 G4long numOfCollisions=0; 772 768 773 // G4cout<<"G4PAIModel::SampleFluctuations"<<G4endl ; 774 //G4cout<<"in: "<<fMaterialCutsCoupleVector[jMat]->GetMaterial()->GetName()<<G4endl ; 775 769 //G4cout<<"G4PAIModel::SampleFluctuations"<<G4endl ; 770 //G4cout<<"in: "<<fMaterialCutsCoupleVector[jMat]->GetMaterial()->GetName()<<G4endl; 776 771 G4double loss = 0.0, charge2 ; 777 772 G4double stepSum = 0., stepDelta, lambda, omega; … … 784 779 G4double TkinScaled = Tkin*MassRatio ; 785 780 786 for(iTkin=0;iTkin< fTotBin;iTkin++)781 for(iTkin=0;iTkin<=fTotBin;iTkin++) 787 782 { 788 783 if(TkinScaled < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ; … … 790 785 iPlace = iTkin - 1 ; 791 786 if(iPlace < 0) iPlace = 0; 792 // G4cout<<"from search, iPlace = "<<iPlace<<G4endl ; 787 else if(iPlace >= fTotBin) iPlace = fTotBin - 1; 788 //G4cout<<"from search, iPlace = "<<iPlace<<G4endl ; 793 789 dNdxCut1 = (*fdNdxCutVector)(iPlace) ; 794 // G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ; 795 790 //G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ; 796 791 797 792 if(iTkin == fTotBin) // Fermi plato, try from left … … 810 805 numOfCollisions++; 811 806 } 812 // 807 //G4cout<<"##1 numOfCollisions = "<<numOfCollisions<<G4endl ; 813 808 814 809 while(numOfCollisions) … … 823 818 } 824 819 omega = GetEnergyTransfer(iPlace,position,iTransfer); 825 // 820 //G4cout<<"G4PAIModel::SampleFluctuations, omega = "<<omega/keV<<" keV; "<<"\t"; 826 821 loss += omega; 827 822 numOfCollisions-- ; … … 831 826 { 832 827 dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ; 833 // G4cout<<"dNdxCut2 = "<<dNdxCut2<<G4endl;828 //G4cout<<"dNdxCut2 = "<<dNdxCut2<< " iTkin= "<<iTkin<<" iPlace= "<<iPlace<<G4endl; 834 829 835 830 if(iTkin == 0) // Tkin is too small, trying from right only … … 839 834 // numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ; 840 835 // numOfCollisions = G4Poisson(meanNumber) ; 841 if( meanNumber > 0.) lambda = step/meanNumber;842 else lambda = DBL_MAX;843 while(numb)844 845 846 847 848 849 850 851 //G4cout<<"##2 numOfCollisions = "<<numOfCollisions<<G4endl ;836 if( meanNumber > 0.) lambda = step/meanNumber; 837 else lambda = DBL_MAX; 838 while(numb) 839 { 840 stepDelta = CLHEP::RandExponential::shoot(lambda); 841 stepSum += stepDelta; 842 if(stepSum >= step) break; 843 numOfCollisions++; 844 } 845 846 //G4cout<<"##2 numOfCollisions = "<<numOfCollisions<<G4endl ; 852 847 853 848 while(numOfCollisions) … … 862 857 } 863 858 omega = GetEnergyTransfer(iPlace,position,iTransfer); 864 // 859 //G4cout<<omega/keV<<"\t"; 865 860 loss += omega; 866 861 numOfCollisions-- ; … … 875 870 W2 = (TkinScaled - E1)*W ; 876 871 877 // G4cout<<"(*(*fPAItransferTable)(iPlace))(0) = "<< 872 //G4cout << fPAItransferTable->size() << G4endl; 873 //G4cout<<"(*(*fPAItransferTable)(iPlace))(0) = "<< 878 874 // (*(*fPAItransferTable)(iPlace))(0)<<G4endl ; 879 // 875 //G4cout<<"(*(*fPAItransferTable)(iPlace+1))(0) = "<< 880 876 // (*(*fPAItransferTable)(iPlace+1))(0)<<G4endl ; 881 877 … … 885 881 // numOfCollisions = RandPoisson::shoot(meanNumber) ; 886 882 // numOfCollisions = G4Poisson(meanNumber) ; 887 if( meanNumber > 0.) lambda = step/meanNumber;888 else lambda = DBL_MAX;889 while(numb)890 891 892 893 894 895 896 897 //G4cout<<"##3 numOfCollisions = "<<numOfCollisions<<endl ;883 if( meanNumber > 0.) lambda = step/meanNumber; 884 else lambda = DBL_MAX; 885 while(numb) 886 { 887 stepDelta = CLHEP::RandExponential::shoot(lambda); 888 stepSum += stepDelta; 889 if(stepSum >= step) break; 890 numOfCollisions++; 891 } 892 893 //G4cout<<"##3 numOfCollisions = "<<numOfCollisions<<endl ; 898 894 899 895 while(numOfCollisions) … … 923 919 } 924 920 } 925 // 921 //G4cout<<"PAIModel AlongStepLoss = "<<loss/keV<<" keV, on step = " 926 922 // <<step/mm<<" mm"<<G4endl ; 927 923 if(loss > Tkin) loss=Tkin; -
trunk/source/processes/electromagnetic/standard/src/G4PAIPhotonModel.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4PAIPhotonModel.cc,v 1.2 2 2009/04/09 18:41:18vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4PAIPhotonModel.cc,v 1.23 2009/07/26 15:51:01 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 310 310 deltaLow = 100.*eV; // 0.5*eV ; 311 311 312 for (G4int i = 0 ; i < fTotBin ; i++) //The loop for the kinetic energy312 for (G4int i = 0 ; i <= fTotBin ; i++) //The loop for the kinetic energy 313 313 { 314 314 LowEdgeEnergy = fProtonEnergyVector->GetLowEdgeEnergy(i) ; … … 441 441 <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl; 442 442 } 443 for ( i = 0 ; i < fTotBin ; i++ )443 for ( i = 0 ; i <= fTotBin ; i++ ) 444 444 { 445 445 dNdxPhotonCut = GetdNdxPhotonCut(i,photonCutInKineticEnergyNow); … … 664 664 fPAIdEdxTable = fPAIdEdxBank[jMat]; 665 665 fdEdxVector = fdEdxTable[jMat]; 666 for(iTkin = 0 ; iTkin < fTotBin ; iTkin++)666 for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++) 667 667 { 668 668 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ; … … 722 722 fPAIplasmonTable = fPAIplasmonBank[jMat]; 723 723 724 for(iTkin = 0 ; iTkin < fTotBin ; iTkin++)724 for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++) 725 725 { 726 726 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ; … … 789 789 790 790 G4int iTkin; 791 for(iTkin=0;iTkin< fTotBin;iTkin++)791 for(iTkin=0;iTkin<=fTotBin;iTkin++) 792 792 { 793 793 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ; … … 1052 1052 G4double cof = step*charge2; 1053 1053 1054 for( iTkin = 0; iTkin < fTotBin; iTkin++)1054 for( iTkin = 0; iTkin <= fTotBin; iTkin++) 1055 1055 { 1056 1056 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ; -
trunk/source/processes/electromagnetic/standard/src/G4PAIxSection.cc
r1007 r1196 26 26 // 27 27 // $Id: G4PAIxSection.cc,v 1.24 2008/05/30 16:04:40 grichine Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // -
trunk/source/processes/electromagnetic/standard/src/G4PAIySection.cc
r1007 r1196 25 25 // 26 26 // 27 // $Id: G4PAIySection.cc,v 1. 3 2007/10/01 18:38:10vnivanch Exp $28 // GEANT4 tag $Name: geant4-09-0 2$27 // $Id: G4PAIySection.cc,v 1.4 2009/07/26 15:51:01 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // … … 41 41 // 42 42 // 01.10.07, V.Ivanchenko create using V.Grichine G4PAIxSection class 43 // 26.07.09, V.Ivanchenko added protection for mumerical exceptions for 44 // low-density materials 43 45 // 44 46 … … 128 130 // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t" 129 131 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl ; 130 if(fEnergyInterval[i+1]-fEnergyInterval[i] >132 if(fEnergyInterval[i+1]-fEnergyInterval[i] < 131 133 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i])) 132 {133 continue ;134 }135 else136 134 { 137 135 for(j=i;j<fIntervalNumber;j++) … … 529 527 // result *= (1-exp(-be2/betaBohr2)) ; 530 528 result *= (1-exp(-be4/betaBohr4)) ; 531 if(fDensity >= 0.1) 529 // if(fDensity >= 0.1) 530 if(x8 > 0.) 532 531 { 533 532 result /= x8 ; … … 586 585 dNdxC *= (1-exp(-be4/betaBohr4)) ; 587 586 588 if(fDensity >= 0.1)589 {590 587 // if(fDensity >= 0.1) 588 // { 589 modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) + 591 590 fImPartDielectricConst[i]*fImPartDielectricConst[i] ; 592 dNdxC /= modul2 ; 593 } 591 if(modul2 > 0.) 592 { 593 dNdxC /= modul2 ; 594 } 594 595 return dNdxC ; 595 596 … … 626 627 dNdxP *= (1-exp(-be4/betaBohr4)) ; 627 628 628 if( fDensity >= 0.1 ) 629 { 630 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 631 fImPartDielectricConst[i]*fImPartDielectricConst[i] ; 632 dNdxP /= modul2 ; 633 } 629 // if( fDensity >= 0.1 ) 630 // { 631 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 632 fImPartDielectricConst[i]*fImPartDielectricConst[i] ; 633 if(modul2 > 0.) 634 { 635 dNdxP /= modul2 ; 636 } 634 637 return dNdxP ; 635 638 … … 800 803 G4double G4PAIySection::SumOverInterCerenkov( G4int i ) 801 804 { 802 G4double x0,x1,y0,yy1,a, b,c,result ;805 G4double x0,x1,y0,yy1,a,c,result ; 803 806 804 807 x0 = fSplineEnergy[i] ; … … 811 814 c = x1/x0; 812 815 a = log10(yy1/y0)/log10(c) ; 813 b = y0/pow(x0,a) ; 816 G4double b = 0.0; 817 if(a < 20.) b = y0/pow(x0,a) ; 814 818 815 819 a += 1.0 ; … … 833 837 G4double G4PAIySection::SumOverInterPlasmon( G4int i ) 834 838 { 835 G4double x0,x1,y0,yy1,a, b,c,result ;839 G4double x0,x1,y0,yy1,a,c,result ; 836 840 837 841 x0 = fSplineEnergy[i] ; … … 841 845 c =x1/x0; 842 846 a = log10(yy1/y0)/log10(c) ; 843 // b = log10(y0) - a*log10(x0) ; 844 b = y0/pow(x0,a) ; 847 848 G4double b = 0.0; 849 if(a < 20.) b = y0/pow(x0,a) ; 845 850 846 851 a += 1.0 ; … … 864 869 G4double en0 ) 865 870 { 866 G4double x0,x1,y0,yy1,a, b,c,d,e0,result ;871 G4double x0,x1,y0,yy1,a,c,d,e0,result ; 867 872 868 873 e0 = en0 ; … … 875 880 d = e0/x0; 876 881 a = log10(yy1/y0)/log10(x1/x0) ; 877 // b0 = log10(y0) - a*log10(x0) ; 878 b = y0/pow(x0,a); // pow(10.,b) ; 882 883 G4double b = 0.0; 884 if(a < 20.) b = y0/pow(x0,a) ; 879 885 880 886 a += 1 ; … … 933 939 G4double en0 ) 934 940 { 935 G4double x0,x1,y0,yy1,a, b,c,d,e0,result ;941 G4double x0,x1,y0,yy1,a,c,d,e0,result ; 936 942 937 943 e0 = en0 ; … … 944 950 d = e0/x0; 945 951 a = log10(yy1/y0)/log10(x1/x0) ; 946 // b0 = log10(y0) - a*log10(x0) ; 947 b = y0/pow(x0,a); // pow(10.,b) ; 952 953 G4double b = 0.0; 954 if(a < 20.) b = y0/pow(x0,a) ; 948 955 949 956 a += 2 ; … … 964 971 d = e0/x0; 965 972 a = log10(yy1/y0)/log10(x1/x0) ; 966 // b0 = log10(y0) - a*log10(x0) ; 967 b = y0/pow(x0,a) ; 973 974 if(a < 20.) b = y0/pow(x0,a) ; 975 968 976 a += 2 ; 969 977 if(a == 0) … … 987 995 G4double en0 ) 988 996 { 989 G4double x0,x1,y0,yy1,a, b,e0,c,d,result ;997 G4double x0,x1,y0,yy1,a,e0,c,d,result ; 990 998 991 999 e0 = en0 ; … … 996 1004 997 1005 // G4cout<<G4endl; 998 // 1006 //G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1 999 1007 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl; 1000 1008 c = x1/x0 ; 1001 1009 d = e0/x0 ; 1002 1010 a = log10(yy1/y0)/log10(c) ; 1003 // b0 = log10(y0) - a*log10(x0) ; 1004 b = y0/pow(x0,a); // pow(10.,b0) ; 1011 1012 G4double b = 0.0; 1013 if(a < 20.) b = y0/pow(x0,a) ; 1005 1014 1006 1015 a += 1.0 ; … … 1012 1021 else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ; 1013 1022 1014 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;1023 //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl; 1015 1024 1016 1025 x0 = fSplineEnergy[i - 1] ; … … 1019 1028 yy1 = fdNdxCerenkov[i - 2] ; 1020 1029 1021 // 1030 //G4cout<<"x0 ="<<x0<<"; x1 = "<<x1 1022 1031 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl; 1023 1032 … … 1025 1034 d = e0/x0 ; 1026 1035 a = log10(yy1/y0)/log10(x1/x0) ; 1027 // b0 = log10(y0) - a*log10(x0) ; 1028 b = y0/pow(x0,a); // pow(10.,b0) ; 1036 1037 // G4cout << "a= " << a << G4endl; 1038 if(a < 20.) b = y0/pow(x0,a) ; 1039 1040 if(a > 20.0) b = 0.0; 1041 else b = y0/pow(x0,a); // pow(10.,b0) ; 1042 1043 //G4cout << "b= " << b << G4endl; 1029 1044 1030 1045 a += 1.0 ; … … 1032 1047 else result += y0*(e0*pow(d,a-1) - x0 )/a ; 1033 1048 a += 1.0 ; 1049 //G4cout << "result= " << result << G4endl; 1034 1050 1035 1051 if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0) ; 1036 1052 else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ; 1037 1053 1038 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = " 1039 // <<b<<"; result = "<<result<<G4endl; 1054 //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl; 1040 1055 1041 1056 return result ; … … 1051 1066 G4double en0 ) 1052 1067 { 1053 G4double x0,x1,y0,yy1,a, b,c,d,e0,result ;1068 G4double x0,x1,y0,yy1,a,c,d,e0,result ; 1054 1069 1055 1070 e0 = en0 ; … … 1062 1077 d = e0/x0 ; 1063 1078 a = log10(yy1/y0)/log10(c) ; 1064 // b0 = log10(y0) - a*log10(x0) ; 1065 b = y0/pow(x0,a); //pow(10.,b) ; 1079 1080 G4double b = 0.0; 1081 if(a < 20.) b = y0/pow(x0,a) ; 1066 1082 1067 1083 a += 1.0 ; … … 1081 1097 d = e0/x0 ; 1082 1098 a = log10(yy1/y0)/log10(c) ; 1083 // b0 = log10(y0) - a*log10(x0) ;1084 b = y0/pow(x0,a);// pow(10.,b0) ;1099 1100 if(a < 20.) b = y0/pow(x0,a) ; 1085 1101 1086 1102 a += 1.0 ; -
trunk/source/processes/electromagnetic/standard/src/G4PEEffectModel.cc
r1055 r1196 25 25 // 26 26 // $Id: G4PEEffectModel.cc,v 1.8 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4PSTARStopping.cc
r1007 r1196 25 25 // 26 26 // $Id: G4PSTARStopping.cc,v 1.8 2008/11/24 18:28:09 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-0 2$27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 29 29 //--------------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4PairProductionRelModel.cc
r1058 r1196 25 25 // 26 26 // $Id: G4PairProductionRelModel.cc,v 1.3 2009/05/15 17:12:33 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4PhotoElectricEffect.cc
r1055 r1196 25 25 // 26 26 // $Id: G4PhotoElectricEffect.cc,v 1.42 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // -
trunk/source/processes/electromagnetic/standard/src/G4PolarizedComptonScattering.cc
r1007 r1196 26 26 // 27 27 // $Id: G4PolarizedComptonScattering.cc,v 1.18 2008/10/15 17:53:44 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-0 2$28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // -
trunk/source/processes/electromagnetic/standard/src/G4UniversalFluctuation.cc
r1055 r1196 25 25 // 26 26 // $Id: G4UniversalFluctuation.cc,v 1.22 2009/03/20 18:11:23 urban Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel.cc
r1055 r1196 25 25 // 26 26 // 27 // $Id: G4UrbanMscModel.cc,v 1.9 0 2009/04/29 13:30:22vnivanch Exp $28 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // $Id: G4UrbanMscModel.cc,v 1.91 2009/07/20 18:41:34 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------------- … … 122 122 lambdalimit = 1.*mm; 123 123 fr = 0.02; 124 facsafety = 0.3;124 // facsafety = 0.3; 125 125 taubig = 8.0; 126 126 tausmall = 1.e-16; -
trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel2.cc
r1055 r1196 25 25 // 26 26 // 27 // $Id: G4UrbanMscModel2.cc,v 1.2 6 2009/05/19 06:26:10 urbanExp $28 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // $Id: G4UrbanMscModel2.cc,v 1.27 2009/07/20 18:41:34 vnivanch Exp $ 28 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 29 29 // 30 30 // ------------------------------------------------------------------- … … 119 119 lambdalimit = 1.*mm; 120 120 fr = 0.02; 121 facsafety = 0.3;121 //facsafety = 0.3; 122 122 taubig = 8.0; 123 123 tausmall = 1.e-16; … … 462 462 presafety = sp->GetSafety(); 463 463 464 // G4cout << "G4UrbanMscModel2::ComputeTruePathLengthLimit tPathLength= " 465 // <<tPathLength<<" safety= " << presafety 466 // << " range= " <<currentRange<<G4endl; 464 //G4cout << "G4Urban2::StepLimit tPathLength= " 465 // <<tPathLength<<" safety= " << presafety 466 // << " range= " <<currentRange<< " lambda= "<<lambda0 467 // << " Alg: " << steppingAlgorithm <<G4endl; 467 468 468 469 // far from geometry boundary … … 506 507 tlimitmin = 10.*stepmin; 507 508 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix; 508 509 //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin 510 // << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl; 509 511 // constraint from the geometry 510 512 if((geomlimit < geombig) && (geomlimit > geommin)) … … 535 537 if(tlimit > tgeom) tlimit = tgeom; 536 538 537 // 538 // << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;539 //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit 540 // << " tlimit= " << tlimit << " presafety= " << presafety << G4endl; 539 541 540 542 // shortcut … … 627 629 } 628 630 } 629 // G4cout << "tPathLength= " << tPathLength << " geomlimit= " << geomlimit631 //G4cout << "tPathLength= " << tPathLength 630 632 // << " currentMinimalStep= " << currentMinimalStep << G4endl; 631 633 return tPathLength ; … … 713 715 714 716 if(zPathLength > lambda0) zPathLength = lambda0; 715 717 //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda0 << G4endl; 716 718 return zPathLength; 717 719 } … … 744 746 } 745 747 if(tPathLength < geomStepLength) tPathLength = geomStepLength; 748 //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl; 746 749 747 750 return tPathLength; … … 799 802 800 803 G4double r = SampleDisplacement(); 801 /* 804 /* 802 805 G4cout << "G4UrbanMscModel2::SampleSecondaries: e(MeV)= " << kineticEnergy 803 806 << " sinTheta= " << sth << " r(mm)= " << r -
trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel90.cc
r1055 r1196 25 25 // 26 26 // $Id: G4UrbanMscModel90.cc,v 1.13 2009/04/10 18:10:58 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4WaterStopping.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4WaterStopping.cc,v 1.1 6 2009/05/15 19:04:21vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4WaterStopping.cc,v 1.18 2009/06/19 10:39:48 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 29 29 //--------------------------------------------------------------------------- … … 67 67 { 68 68 G4double res = 0.0; 69 if( iz < 3 || iz > 18) return res;69 if((iz > 26) || (iz < 3) || (iz > 18 && iz < 26)) return res; 70 70 G4bool b; 71 71 G4int idx = iz - 3; … … 99 99 G4int i; 100 100 //..List of ions 101 G4int zz[1 6] = {3, 4, 5, 6, 7, 8, 9, 10,11,12,13,14,15,16,17,18};102 G4int aa[1 6] = {7, 9, 11, 12, 14, 16, 19, 20, 23, 24, 27, 28,31,32, 35,40};103 // G4double A_Ion[17] = {6.941,9.0122,10.811,12.011,14.007,15.999,18.998,20.180,22.990,24.305,26.982,28.086,30.974,32.065,35.453,39.948 };104 for(i=0; i<1 6; i++) {101 G4int zz[17] = {3, 4, 5, 6, 7, 8, 9, 10,11,12,13,14,15,16,17,18,26}; 102 G4int aa[17] = {7, 9, 11, 12, 14, 16, 19, 20, 23, 24, 27, 28,31,32, 35,40,56}; 103 // G4double A_Ion[17] = {6.941,9.0122,10.811,12.011,14.007,15.999,18.998,20.180,22.990,24.305,26.982,28.086,30.974,32.065,35.453,39.948,55.845}; 104 for(i=0; i<17; i++) { 105 105 Z[i] = zz[i]; 106 106 A[i] = G4double(aa[i]); … … 145 145 AddData(E,G4_WATER_Ar,factor); 146 146 147 G4double G4_WATER_Fe [53]={6.5394, 7.3060, 8.7367, 10.0690, 11.3310, 12.5470, 13.7280, 14.8780, 15.9980, 21.1160, 25.4850, 29.1260, 32.0640, 36.1770, 38.6920, 40.2760, 41.2950, 41.9530, 42.3710, 42.6210, 42.5910, 41.6900, 40.5190, 39.2690, 36.8000, 34.5040, 32.4190, 30.5410, 28.8480, 27.3170, 25.9310, 20.6170, 17.0680, 14.5540, 12.6930, 10.1410, 8.4892, 7.3402, 6.4976, 5.8545, 5.3479, 4.9387, 3.6892, 3.0503, 2.6620, 2.4014, 2.0756, 1.8825, 1.7569, 1.6702, 1.6079, 1.5619, 1.5267}; 148 AddData(E,G4_WATER_Fe,factor); 149 147 150 if(corr) { 148 for(i=0; i<1 6; i++) {corr->AddStoppingData(Z[i], aa[i], "G4_WATER", dedx[i]);}151 for(i=0; i<17; i++) {corr->AddStoppingData(Z[i], aa[i], "G4_WATER", dedx[i]);} 149 152 } 150 153 } -
trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4WentzelVIModel.cc,v 1.3 2 2009/05/10 16:09:29vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4WentzelVIModel.cc,v 1.37 2009/10/28 10:14:13 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 135 135 // reset parameters 136 136 SetupParticle(p); 137 tkin = targetZ = mom2 = DBL_MIN;137 tkin = targetZ = mom2 = 0.0; 138 138 ecut = etag = DBL_MAX; 139 139 currentRange = 0.0; … … 158 158 { 159 159 SetupParticle(p); 160 G4double ekin = std::max(lowEnergyLimit, kinEnergy);161 SetupKinematic( ekin, cutEnergy);162 SetupTarget(Z, ekin);160 if(kinEnergy < lowEnergyLimit) return 0.0; 161 SetupKinematic(kinEnergy, cutEnergy); 162 SetupTarget(Z, kinEnergy); 163 163 G4double xsec = ComputeTransportXSectionPerAtom(); 164 /* 164 /* 165 165 G4cout << "CS: e= " << tkin << " cosEl= " << cosTetMaxElec2 166 166 << " cosN= " << cosTetMaxNuc2 << " xsec(bn)= " << xsec/barn … … 392 392 // << particle->GetParticleName() << G4endl; 393 393 G4double kinEnergy = dynParticle->GetKineticEnergy(); 394 if(kinEnergy <= DBL_MIN || tPathLength <= DBL_MIN) return; 394 395 // ignore scattering for zero step length and enegy below the limit 396 if(kinEnergy < lowEnergyLimit || tPathLength <= DBL_MIN) return; 395 397 396 398 G4double ekin = preKinEnergy; … … 423 425 424 426 // for low-energy e-,e+ no limit 425 ekin = std::max(ekin, lowEnergyLimit);426 427 SetupKinematic(ekin, cut); 427 428 … … 431 432 xsec = ComputeXSectionPerVolume(); 432 433 433 if(xtsec > DBL_MIN) x1 = 0.5*tPathLength*xtsec;434 else 434 if(xtsec > 0.0) x1 = 0.5*tPathLength*xtsec; 435 else x1 = 0.0; 435 436 436 437 /* … … 452 453 // cost is sampled ------------------------------ 453 454 G4double cost = 1.0 - 2.0*z; 454 if(cost < -1.0) cost = -1.0;455 else if(cost > 1.0) cost = 1.0;455 // if(cost < -1.0) cost = -1.0; 456 // else if(cost > 1.0) cost = 1.0; 456 457 G4double sint = sqrt((1.0 - cost)*(1.0 + cost)); 457 458 … … 472 473 473 474 // sample MSC scattering for large angle 474 // extra central scattering for h olf step475 // extra central scattering for half step 475 476 if(largeAng) { 476 477 isscat = true; … … 480 481 } while (z > 1.0); 481 482 cost = 1.0 - 2.0*z; 482 if(std::abs(cost) > 1.0) cost = 1.0;483 //if(std::fabs(cost) > 1.0) cost = 1.0; 483 484 484 485 sint = sqrt((1.0 - cost)*(1.0 + cost)); … … 491 492 } 492 493 493 // sample R eserford scattering for large angle494 // sample Rutherford scattering for large angle 494 495 if(xsec > DBL_MIN) { 495 496 G4double t = tPathLength; … … 536 537 if(zz1 < 1.0) { 537 538 isscat = true; 538 //G4cout << "R eserford zz1= " << zz1 << " t= " << t << G4endl;539 //G4cout << "Rutherford zz1= " << zz1 << " t= " << t << G4endl; 539 540 sint = sqrt((1.0 - zz1)*(1.0 + zz1)); 540 541 //G4cout << "sint= " << sint << G4endl; … … 637 638 if(cosThetaMin > cosnm) { 638 639 639 // R eserford part640 // Rutherford part 640 641 G4double s = screenZ*formfactA; 641 642 G4double z1 = 1.0 - cosnm + screenZ; … … 711 712 zcorr += 0.5*x*s; 712 713 713 // R eserford part714 // Rutherford part 714 715 G4double z1 = 1.0 - cosTetMaxElec + screenZ; 715 716 G4double z2 = (cosThetaMin - cosTetMaxElec)/x1; … … 730 731 zcorr += 0.5*x*s; 731 732 732 // R eserford part733 // Rutherford part 733 734 s = screenZ*formfactA; 734 735 G4double w = 1.0 + 2.0*s; -
trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlung.cc
r1055 r1196 25 25 // 26 26 // $Id: G4eBremsstrahlung.cc,v 1.56 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungModel.cc
r1055 r1196 25 25 // 26 26 // $Id: G4eBremsstrahlungModel.cc,v 1.44 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc
r1055 r1196 25 25 // 26 26 // $Id: G4eBremsstrahlungRelModel.cc,v 1.14 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4eCoulombScatteringModel.cc
r1055 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eCoulombScatteringModel.cc,v 1. 69 2009/05/10 16:09:29vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $26 // $Id: G4eCoulombScatteringModel.cc,v 1.78 2009/10/28 10:14:13 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 39 39 // 40 40 // Modifications: 41 // 41 42 // 01.08.06 V.Ivanchenko extend upper limit of table to TeV and review the 42 43 // logic of building - only elements from G4ElementTable … … 45 46 // 09.10.07 V.Ivanchenko reorganized methods, add cut dependence in scattering off e- 46 47 // 09.06.08 V.Ivanchenko add SelectIsotope and sampling of the recoil ion 48 // 16.06.09 C.Consolandi fixed computation of effective mass 49 // 47 50 // 48 51 // Class Description: … … 63 66 #include "G4Proton.hh" 64 67 #include "G4ParticleTable.hh" 68 #include "G4ProductionCutsTable.hh" 69 #include "G4NucleiProperties.hh" 65 70 66 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... … … 87 92 currentMaterial = 0; 88 93 currentElement = 0; 89 lowEnergyLimit =keV;94 lowEnergyLimit = 0.1*keV; 90 95 G4double p0 = electron_mass_c2*classic_electr_radius; 91 96 coeff = twopi*p0*p0; 92 97 tkin = targetZ = mom2 = DBL_MIN; 93 98 elecXSection = nucXSection = 0.0; 94 recoilThreshold = 100.*keV;99 recoilThreshold = 0.*keV; 95 100 ecut = DBL_MAX; 96 101 particle = 0; … … 130 135 ecut = etag = DBL_MAX; 131 136 cosThetaMin = cos(PolarAngleLimit()); 132 currentCuts = &cuts;137 pCuts = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3); 133 138 //G4cout << "!!! G4eCoulombScatteringModel::Initialise for " 134 139 // << p->GetParticleName() << " cos(TetMin)= " << cosThetaMin 135 140 // << " cos(TetMax)= " << cosThetaMax <<G4endl; 141 // G4cout << "cut0= " << cuts[0] << " cut1= " << cuts[1] << G4endl; 136 142 if(!isInitialised) { 137 143 isInitialised = true; … … 185 191 G4double xsec = 0.0; 186 192 SetupParticle(p); 187 G4double ekin = std::max(lowEnergyLimit, kinEnergy);188 SetupKinematic( ekin, cutEnergy);193 if(kinEnergy < lowEnergyLimit) return xsec; 194 SetupKinematic(kinEnergy, cutEnergy); 189 195 if(cosTetMaxNuc < cosTetMinNuc) { 190 SetupTarget(Z, ekin);196 SetupTarget(Z, kinEnergy); 191 197 xsec = CrossSectionPerAtom(); 192 198 } … … 196 202 << " cosTetMaxElec= " << cosTetMaxElec 197 203 << " screenZ= " << screenZ 198 << " formfactA= " << formfactA 199 << " cosTetMaxHad= " << cosTetMaxHad << G4endl; 204 << " formfactA= " << formfactA << G4endl; 200 205 */ 201 206 return xsec; … … 207 212 { 208 213 // This method needs initialisation before be called 209 G4double fac = coeff*targetZ*chargeSquare*kinFactor; 214 //G4double fac = coeff*targetZ*chargeSquare*invbeta2/mom2; 215 216 G4double meff = targetMass/(mass+targetMass); 217 G4double fac = coeff*targetZ*chargeSquare*invbeta2/(mom2*meff*meff); 218 210 219 elecXSection = 0.0; 211 220 nucXSection = 0.0; … … 257 266 { 258 267 G4double kinEnergy = dp->GetKineticEnergy(); 259 if(kinEnergy < = DBL_MIN) return;268 if(kinEnergy < lowEnergyLimit) return; 260 269 DefineMaterial(couple); 261 270 SetupParticle(dp->GetDefinition()); 262 G4double ekin = std::max(lowEnergyLimit, kinEnergy); 263 SetupKinematic( ekin, cutEnergy);271 272 SetupKinematic(kinEnergy, cutEnergy); 264 273 //G4cout << "G4eCoulombScatteringModel::SampleSecondaries e(MeV)= " 265 // << kinEnergy << " " << particle->GetParticleName() << G4endl; 274 // << kinEnergy << " " << particle->GetParticleName() 275 // << " cut= " << cutEnergy<< G4endl; 266 276 267 277 // Choose nucleus 268 currentElement = SelectRandomAtom(couple,particle,ekin,cutEnergy,ekin); 269 270 SetupTarget(currentElement->GetZ(),ekin); 278 currentElement = SelectRandomAtom(couple,particle, 279 kinEnergy,cutEnergy,kinEnergy); 280 281 SetupTarget(currentElement->GetZ(),kinEnergy); 282 283 G4int ia = SelectIsotopeNumber(currentElement); 284 targetMass = G4NucleiProperties::GetNuclearMass(ia, iz); 271 285 272 286 G4double cost = SampleCosineTheta(); … … 276 290 G4double sint = sqrt(z1*(1.0 + cost)); 277 291 278 //G4cout<<"## Sampled sint= " << sint << " Z= " << targetZ 292 //G4cout<<"## Sampled sint= " << sint << " Z= " << targetZ << " A= " << ia 279 293 // << " screenZ= " << screenZ << " cn= " << formfactA << G4endl; 280 294 … … 289 303 // recoil sampling assuming a small recoil 290 304 // and first order correction to primary 4-momentum 291 if(lowEnergyLimit < kinEnergy) { 292 G4int ia = SelectIsotopeNumber(currentElement); 293 G4double Trec = z1*mom2/(amu_c2*G4double(ia)); 294 295 if(Trec > recoilThreshold) { 296 G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz); 297 Trec = z1*mom2/ion->GetPDGMass(); 298 if(Trec < kinEnergy) { 299 G4ThreeVector dir = (direction - newDirection).unit(); 300 G4DynamicParticle* newdp = new G4DynamicParticle(ion, dir, Trec); 301 fvect->push_back(newdp); 302 fParticleChange->SetProposedKineticEnergy(kinEnergy - Trec); 303 } 304 } 305 G4double q2 = 2*z1*mom2; 306 G4double trec = q2/(sqrt(targetMass*targetMass + q2) + targetMass); 307 G4double finalT = kinEnergy - trec; 308 //G4cout<<"G4eCoulombScatteringModel: finalT= "<<finalT<<" Trec= "<<trec<<G4endl; 309 if(finalT <= lowEnergyLimit) { 310 trec = kinEnergy; 311 finalT = 0.0; 312 } 313 314 fParticleChange->SetProposedKineticEnergy(finalT); 315 G4double tcut = recoilThreshold; 316 if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); } 317 318 if(trec > tcut) { 319 G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz); 320 G4ThreeVector dir = (direction*sqrt(mom2) - 321 newDirection*sqrt(finalT*(2*mass + finalT))).unit(); 322 G4DynamicParticle* newdp = new G4DynamicParticle(ion, dir, trec); 323 fvect->push_back(newdp); 324 } else { 325 fParticleChange->ProposeLocalEnergyDeposit(trec); 326 fParticleChange->ProposeNonIonizingEnergyDeposit(trec); 305 327 } 306 328 … … 324 346 } 325 347 326 /* 348 /* 327 349 G4cout << "SampleCost: e(MeV)= " << tkin 328 << " ctmin= " << cosThetaMin 329 << " ctmaxN= " << cosTetMaxNuc 330 << " ctmax= " << costm 331 << " Z= " << targetZ << " A= " << targetA 350 << " 1-ctmaxN= " << 1. - cosTetMinNuc 351 << " 1-ctmax= " << 1. - costm 352 << " Z= " << targetZ 332 353 << G4endl; 333 354 */ 355 334 356 if(costm >= cosTetMinNuc) return 2.0; 335 357 … … 343 365 } while ( G4UniformRand() > grej*grej ); 344 366 345 //G4cout << "z= " << z1 << " cross= " << nucXSection/barn 346 // << " crossE= " << elecXSection/barn << G4endl; 367 if(mass > MeV) { 368 if(G4UniformRand() > (1. - z1*0.5)/(1.0 + z1*sqrt(mom2)/targetMass)) { 369 return 2.0; 370 } 371 } 372 //G4cout << "z1= " << z1 << " cross= " << nucXSection/barn 373 // << " crossE= " << elecXSection/barn << G4endl; 347 374 348 375 return 1.0 - z1; -
trunk/source/processes/electromagnetic/standard/src/G4eIonisation.cc
r1055 r1196 25 25 // 26 26 // $Id: G4eIonisation.cc,v 1.57 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4eMultipleScattering.cc
r1007 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4eMultipleScattering.cc,v 1. 7 2008/10/23 17:55:20vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4eMultipleScattering.cc,v 1.10 2009/11/01 13:05:01 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ----------------------------------------------------------------------------- … … 45 45 46 46 #include "G4eMultipleScattering.hh" 47 #include "G4UrbanMscModel 2.hh"47 #include "G4UrbanMscModel92.hh" 48 48 #include "G4MscStepLimitType.hh" 49 49 #include "G4Electron.hh" … … 58 58 { 59 59 isInitialized = false; 60 SetRangeFactor(0.04);61 60 } 62 61 … … 76 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 77 76 78 void G4eMultipleScattering::InitialiseProcess(const G4ParticleDefinition* p)77 void G4eMultipleScattering::InitialiseProcess(const G4ParticleDefinition*) 79 78 { 80 // Modification of parameters between runs 81 if(isInitialized) { 82 if (p->GetParticleType() != "nucleus") { 83 mscUrban->SetStepLimitType(StepLimitType()); 84 mscUrban->SetLateralDisplasmentFlag(LateralDisplasmentFlag()); 85 mscUrban->SetSkin(Skin()); 86 mscUrban->SetRangeFactor(RangeFactor()); 87 mscUrban->SetGeomFactor(GeomFactor()); 88 } 89 return; 90 } 91 92 // defaults for ions, which cannot be overwritten 93 if (p->GetParticleType() == "nucleus") { 94 SetStepLimitType(fMinimal); 95 SetLateralDisplasmentFlag(false); 96 SetBuildLambdaTable(false); 97 } 79 if(isInitialized) return; 98 80 99 81 // initialisation of parameters - defaults for particles other 100 82 // than ions can be overwritten by users 101 mscUrban = new G4UrbanMscModel2(); 102 mscUrban->SetStepLimitType(StepLimitType()); 103 mscUrban->SetLateralDisplasmentFlag(LateralDisplasmentFlag()); 104 mscUrban->SetSkin(Skin()); 105 mscUrban->SetRangeFactor(RangeFactor()); 106 mscUrban->SetGeomFactor(GeomFactor()); 107 83 G4VMscModel* mscUrban = new G4UrbanMscModel92(); 108 84 AddEmModel(1,mscUrban); 109 85 isInitialized = true; -
trunk/source/processes/electromagnetic/standard/src/G4eeToTwoGammaModel.cc
r1055 r1196 25 25 // 26 26 // $Id: G4eeToTwoGammaModel.cc,v 1.15 2009/04/09 18:41:18 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4eplusAnnihilation.cc
r1055 r1196 25 25 // 26 26 // $Id: G4eplusAnnihilation.cc,v 1.30 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4hIonisation.cc
r1055 r1196 25 25 // 26 26 // $Id: G4hIonisation.cc,v 1.82 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- -
trunk/source/processes/electromagnetic/standard/src/G4hMultipleScattering.cc
r1007 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4hMultipleScattering.cc,v 1.1 3 2008/10/15 17:53:44vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4hMultipleScattering.cc,v 1.16 2009/11/01 13:05:01 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ----------------------------------------------------------------------------- … … 59 59 { 60 60 isInitialized = false; 61 isIon = false;62 61 SetStepLimitType(fMinimal); 63 62 } … … 77 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 78 77 79 void G4hMultipleScattering::InitialiseProcess(const G4ParticleDefinition* p)78 void G4hMultipleScattering::InitialiseProcess(const G4ParticleDefinition*) 80 79 { 81 // Modification of parameters between runs 82 if(isInitialized) { 83 if (p->GetParticleType() != "nucleus" && p->GetPDGMass() < GeV) { 84 mscUrban->SetStepLimitType(StepLimitType()); 85 mscUrban->SetLateralDisplasmentFlag(LateralDisplasmentFlag()); 86 mscUrban->SetSkin(Skin()); 87 mscUrban->SetRangeFactor(RangeFactor()); 88 mscUrban->SetGeomFactor(GeomFactor()); 89 } 90 return; 91 } 92 93 // defaults for ions, which cannot be overwritten 94 if (p->GetParticleType() == "nucleus" || p->GetPDGMass() > GeV) { 95 SetStepLimitType(fMinimal); 96 SetLateralDisplasmentFlag(false); 97 SetBuildLambdaTable(false); 98 if(p->GetParticleType() == "nucleus") isIon = true; 99 } 100 101 // initialisation of parameters 102 G4String part_name = p->GetParticleName(); 103 mscUrban = new G4UrbanMscModel90(); 104 105 mscUrban->SetStepLimitType(StepLimitType()); 106 mscUrban->SetLateralDisplasmentFlag(LateralDisplasmentFlag()); 107 mscUrban->SetSkin(Skin()); 108 mscUrban->SetRangeFactor(RangeFactor()); 109 mscUrban->SetGeomFactor(GeomFactor()); 110 80 if(isInitialized) return; 81 G4VMscModel* mscUrban = new G4UrbanMscModel90(); 111 82 AddEmModel(1,mscUrban); 112 83 isInitialized = true; … … 118 89 { 119 90 G4cout << " RangeFactor= " << RangeFactor() 120 << ", step limit type: " << StepLimitType()121 << ", lat eralDisplacement: " << LateralDisplasmentFlag()122 << ", skin= " << Skin() 123 //<< ", geomFactor= " << GeomFactor()91 << ", stepLimitType: " << StepLimitType() 92 << ", latDisplacement: " << LateralDisplasmentFlag() 93 << ", skin= " << Skin() 94 << ", geomFactor= " << GeomFactor() 124 95 << G4endl; 125 96 } … … 127 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 128 99 129 G4double G4hMultipleScattering::AlongStepGetPhysicalInteractionLength(130 const G4Track& track,131 G4double,132 G4double currentMinimalStep,133 G4double& currentSafety,134 G4GPILSelection* selection)135 {136 // get Step limit proposed by the process137 valueGPILSelectionMSC = NotCandidateForSelection;138 139 G4double escaled = track.GetKineticEnergy();140 if(isIon) escaled *= track.GetDynamicParticle()->GetMass()/proton_mass_c2;141 142 G4double steplength = GetMscContinuousStepLimit(track,143 escaled,144 currentMinimalStep,145 currentSafety);146 // G4cout << "StepLimit= " << steplength << G4endl;147 // set return value for G4GPILSelection148 *selection = valueGPILSelectionMSC;149 return steplength;150 }151 152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......153 -
trunk/source/processes/electromagnetic/standard/src/G4ionGasIonisation.cc
r1007 r1196 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4ionGasIonisation.cc,v 1.1 4 2008/09/12 16:26:34vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 2$26 // $Id: G4ionGasIonisation.cc,v 1.15 2009/11/10 11:50:30 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 58 58 : G4ionIonisation(name) 59 59 { 60 G4cout << G4endl; 61 G4cout << "!!! G4ionGasIonisation class is obsolete and may be removed for the next major Geant4 release !!!" << G4endl; 62 G4cout << "!!! Please use G4ionIonisation for ions !!!" << G4endl; 63 G4cout << G4endl; 60 64 } 61 65 -
trunk/source/processes/electromagnetic/standard/src/G4ionIonisation.cc
r1055 r1196 25 25 // 26 26 // $Id: G4ionIonisation.cc,v 1.66 2009/02/20 12:06:37 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-03- beta-cand-01 $27 // GEANT4 tag $Name: geant4-09-03-cand-01 $ 28 28 // 29 29 // -------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.