Changeset 1340 for trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc

Timestamp:

Nov 5, 2010, 3:45:55 PM (14 years ago)

Author:

garnier

Message:

update ti head

File:

• trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc (modified) (18 diffs)

trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4eBremsstrahlungRelModel.cc,v 1.15 2010/04/06 17:02:23 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4eBremsstrahlungRelModel.cc,v 1.18 2010/11/04 17:30:32 vnivanch Exp $
+// GEANT4 tag $Name: emstand-V09-03-25 $
 //
 // -------------------------------------------------------------------
@@ -42 +42 @@
 // 13.11.08    add SetLPMflag and SetLPMconstant methods
 // 13.11.08    change default LPMconstant value
+// 13.10.10    add angular distributon interface (VI)
 //
 // Main References:
@@ -65 +66 @@
 #include "G4ParticleChangeForLoss.hh"
 #include "G4LossTableManager.hh"
-
+#include "G4ModifiedTsai.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -90 +91 @@
     use_completescreening(true),isInitialised(false)
 {
-  if(p) SetParticle(p);
   theGamma = G4Gamma::Gamma();
 
   minThreshold = 0.1*keV;
-  SetLowEnergyLimit(GeV);  
+  lowKinEnergy = GeV;
+  SetLowEnergyLimit(lowKinEnergy);  
 
   nist = G4NistManager::Instance();  
+
+  SetLPMFlag(true);
+  SetAngularDistribution(new G4ModifiedTsai());
+
+  particleMass = kinEnergy = totalEnergy = currentZ = z13 = z23 = lnZ = Fel 
+    = Finel = fCoulomb = fMax = densityFactor = densityCorr = lpmEnergy 
+    = xiLPM = phiLPM = gLPM = klpm = kp = 0.0;
+
+  energyThresholdLPM = 1.e39;
+
   InitialiseConstants();
-
-  SetLPMFlag(true);
+  if(p) { SetParticle(p); }
 }
 
@@ -125 +135 @@
   particle = p;
   particleMass = p->GetPDGMass();
-  if(p == G4Electron::Electron()) isElectron = true;
-  else                            isElectron = false;
+  if(p == G4Electron::Electron()) { isElectron = true; }
+  else                            { isElectron = false;}
 }
 
@@ -140 +150 @@
 
 void G4eBremsstrahlungRelModel::SetupForMaterial(const G4ParticleDefinition*,
-                                                 const G4Material* mat, G4double kineticEnergy)
+                                                 const G4Material* mat, 
+                                                 G4double kineticEnergy)
 {
   densityFactor = mat->GetElectronDensity()*fMigdalConstant;
@@ -146 +157 @@
 
   // Threshold for LPM effect (i.e. below which LPM hidden by density effect) 
-  if (LPMFlag()) 
+  if (LPMFlag()) {
     energyThresholdLPM=sqrt(densityFactor)*lpmEnergy;
-  else
+  } else {
      energyThresholdLPM=1.e39;   // i.e. do not use LPM effect
-
+  }
   // calculate threshold for density effect
   kinEnergy   = kineticEnergy;
@@ -168 +179 @@
                                            const G4DataVector& cuts)
 {
-  if(p) SetParticle(p);
-
-  highKinEnergy = HighEnergyLimit();
+  if(p) { SetParticle(p); }
+
   lowKinEnergy  = LowEnergyLimit();
 
@@ -177 +187 @@
   InitialiseElementSelectors(p, cuts);
 
-  if(isInitialised) return;
+  if(isInitialised) { return; }
   fParticleChange = GetParticleChangeForLoss();
   isInitialised = true;
@@ -190 +200 @@
                                                    G4double cutEnergy)
 {
-  if(!particle) SetParticle(p);
-  if(kineticEnergy < lowKinEnergy) return 0.0;
+  if(!particle) { SetParticle(p); }
+  if(kineticEnergy < lowKinEnergy) { return 0.0; }
   G4double cut = std::min(cutEnergy, kineticEnergy);
-  if(cut == 0.0) return 0.0;
+  if(cut == 0.0) { return 0.0; }
 
   SetupForMaterial(particle, material,kineticEnergy);
@@ -261 +271 @@
                                               G4double maxEnergy)
 {
-  if(!particle) SetParticle(p);
-  if(kineticEnergy < lowKinEnergy) return 0.0;
+  if(!particle) { SetParticle(p); }
+  if(kineticEnergy < lowKinEnergy) { return 0.0; }
   G4double cut  = std::min(cutEnergy, kineticEnergy);
   G4double tmax = std::min(maxEnergy, kineticEnergy);
 
-  if(cut >= tmax) return 0.0;
+  if(cut >= tmax) { return 0.0; }
 
   SetCurrentElement(Z);
@@ -273 +283 @@
 
   // allow partial integration
-  if(tmax < kinEnergy) cross -= ComputeXSectionPerAtom(tmax);
+  if(tmax < kinEnergy) { cross -= ComputeXSectionPerAtom(tmax); }
   
   cross *= Z*Z*bremFactor;
@@ -390 +400 @@
   // *** make sure suppression is smaller than 1 ***
   // *** caused by Migdal approximation in xi    ***
-  if (xiLPM*phiLPM>1. || s>0.57)  xiLPM=1./phiLPM;
+  if (xiLPM*phiLPM>1. || s>0.57)  { xiLPM=1./phiLPM; }
 }
 
@@ -401 +411 @@
 //    * complete screening
 {
-  if(gammaEnergy < 0.0) return 0.0;
+  if(gammaEnergy < 0.0) { return 0.0; }
 
   G4double y = gammaEnergy/totalEnergy;
@@ -429 +439 @@
 {
 
-  if(gammaEnergy < 0.0) return 0.0;
+  if(gammaEnergy < 0.0) { return 0.0; }
 
   G4double y = gammaEnergy/totalEnergy;
@@ -465 +475 @@
 {
   G4double kineticEnergy = dp->GetKineticEnergy();
-  if(kineticEnergy < lowKinEnergy) return;
+  if(kineticEnergy < lowKinEnergy) { return; }
   G4double cut  = std::min(cutEnergy, kineticEnergy);
   G4double emax = std::min(maxEnergy, kineticEnergy);
-  if(cut >= emax) return;
+  if(cut >= emax) { return; }
 
   SetupForMaterial(particle, couple->GetMaterial(),kineticEnergy);
@@ -483 +493 @@
   //  G4double fmax= fMax;
   G4bool highe = true;
-  if(totalEnergy < energyThresholdLPM) highe = false;
+  if(totalEnergy < energyThresholdLPM) { highe = false; }
  
   G4double xmin = log(cut*cut + densityCorr);
@@ -507 +517 @@
 
   //
-  //  angles of the emitted gamma. ( Z - axis along the parent particle)
+  // angles of the emitted gamma. ( Z - axis along the parent particle)
+  // use general interface
   //
-  //  universal distribution suggested by L. Urban 
-  // (Geant3 manual (1993) Phys211),
-  //  derived from Tsai distribution (Rev Mod Phys 49,421(1977))
-
-  G4double u;
-  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
-
-  if (9./(9.+d) > G4UniformRand()) u = - log(G4UniformRand()*G4UniformRand())/a1;
-     else                          u = - log(G4UniformRand()*G4UniformRand())/a2;
-
-  G4double theta = u*particleMass/totalEnergy;
+  G4double theta = GetAngularDistribution()->PolarAngle(totalEnergy,
+                                                        totalEnergy-gammaEnergy,
+                                                        (G4int)currentZ);
+
   G4double sint = sin(theta);
   G4double phi = twopi * G4UniformRand();