Changeset 1315 for trunk/source/processes/electromagnetic/utils/src/G4EmCorrections.cc

Timestamp:

Jun 18, 2010, 11:42:07 AM (14 years ago)

Author:

garnier

Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

File:

• trunk/source/processes/electromagnetic/utils/src/G4EmCorrections.cc (modified) (15 diffs)

trunk/source/processes/electromagnetic/utils/src/G4EmCorrections.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4EmCorrections.cc,v 1.54 2009/10/29 17:56:36 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4EmCorrections.cc,v 1.58 2010/06/04 09:28:46 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -68 +68 @@
 #include "G4ProductionCutsTable.hh"
 #include "G4MaterialCutsCouple.hh"
+#include "G4AtomicShells.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -123 +124 @@
   G4double sum = (2.0*(Barkas + Bloch) + Mott);
 
-  if(verbose > 1)
+  if(verbose > 1) {
     G4cout << "EmCorrections: E(MeV)= " << e/MeV << " Barkas= " << Barkas
            << " Bloch= " << Bloch << " Mott= " << Mott 
-           << " Sum= " << sum << G4endl; 
-
+           << " Sum= " << sum << " q2= " << q2 << G4endl; 
+  }
   sum *= material->GetElectronDensity() * q2 *  twopi_mc2_rcl2 /beta2;
   return sum;
@@ -165 +166 @@
 //   Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
   SetupKinematics(p, mat, e);
-  if(tau <= 0.0) return 0.0;
+  if(tau <= 0.0) { return 0.0; }
 
   G4double Barkas = BarkasCorrection (p, mat, e);
@@ -180 +181 @@
   sum *= material->GetElectronDensity() * q2 *  twopi_mc2_rcl2 /beta2;
 
-  if(verbose > 1) G4cout << " Sum= " << sum << G4endl; 
+  if(verbose > 1) { G4cout << " Sum= " << sum << G4endl; } 
   return sum;
 }
@@ -219 +220 @@
     sum = ComputeIonCorrections(p,couple->GetMaterial(),e) - rest/e;
 
-    if(verbose > 1) G4cout << " Sum= " << sum << " dSum= " << rest/e << G4endl; 
+    if(verbose > 1) { G4cout << " Sum= " << sum << " dSum= " << rest/e << G4endl; } 
   }
   return sum;
@@ -267 +268 @@
     }
     G4double e0= 13.6*eV*Z2;
-    term += f*atomDensity[i]*KShell(shells.GetBindingEnergy(iz,0)/e0,ba2/Z2)/Z;
+    term += f*atomDensity[i]*KShell(G4AtomicShells::GetBindingEnergy(iz,0)/e0,ba2/Z2)/Z;
   }
 
@@ -293 +294 @@
       G4double e0= 13.6*eV*Z2*0.25;
       G4double f = 0.125;
-      G4int nmax = std::min(4,shells.GetNumberOfShells(iz));
+      G4int nmax = std::min(4,G4AtomicShells::GetNumberOfShells(iz));
       for(G4int j=1; j<nmax; j++) {
-        G4double ne = G4double(shells.GetNumberOfElectrons(iz,j));
-        G4double e1 = shells.GetBindingEnergy(iz,j);
+        G4double ne = G4double(G4AtomicShells::GetNumberOfElectrons(iz,j));
+        G4double e1 = G4AtomicShells::GetBindingEnergy(iz,j);
         //   G4cout << "LShell: j= " << j << " ne= " << ne << " e(eV)= " << e/eV
         //  << " e0(eV)= " << e0/eV << G4endl;
@@ -417 +418 @@
   for (G4int i = 0; i<numberOfElements; i++) {
 
+    G4double res = 0.0;
     G4double Z = (*theElementVector)[i]->GetZ();
     G4int   iz = G4int(Z);
@@ -426 +428 @@
     }
     G4double e0= 13.6*eV*Z2;
-    term += f*atomDensity[i]*KShell(shells.GetBindingEnergy(iz,0)/e0,ba2/Z2)/Z;
+    res += f*KShell(G4AtomicShells::GetBindingEnergy(iz,0)/e0,ba2/Z2);
     if(2 < iz) {
       G4double Zeff = Z - ZD[10];
@@ -434 +436 @@
       f = 0.125;
       G4double eta = ba2/Z2;
-      G4int ntot = shells.GetNumberOfShells(iz);
+      G4int ntot = G4AtomicShells::GetNumberOfShells(iz);
       G4int nmax = std::min(4, ntot);
       G4double norm   = 0.0;
       G4double eshell = 0.0;
-      for(G4int j=1; j<nmax; j++) {
-        G4double x = G4double(shells.GetNumberOfElectrons(iz,j));
-        G4double e1 = shells.GetBindingEnergy(iz,j);
-        norm   += x;
-        eshell += e1*x;
-        term += f*x*atomDensity[i]*LShell(e1/e0,eta)/Z;
+      for(G4int j=1; j<nmax; ++j) {
+        G4double e1 = G4AtomicShells::GetBindingEnergy(iz,j);
+        G4int ne = G4AtomicShells::GetNumberOfElectrons(iz,j);
+        norm   += ne;
+        eshell += e1*ne;
+        res += f*ne*LShell(e1/e0,eta);
       }
-      if(10 < iz) {
+      if(ntot > nmax) {
         eshell /= norm;
         G4double eeff = eshell*eta;
-        for(G4int k=nmax; k<ntot; k++) {
-          G4double x = G4double(shells.GetNumberOfElectrons(iz,k));
-          G4double e1 = shells.GetBindingEnergy(iz,k);
-          term += f*x*atomDensity[i]*LShell(e1/e0,eeff/e1)/Z;
-          //          term += f*x*atomDensity[i]*LShell(eshell/e0,eeff/e1)/Z;
+        for(G4int k=nmax; k<ntot; ++k) {
+          G4int ne = G4AtomicShells::GetNumberOfElectrons(iz,k);
+          G4double e1 = G4AtomicShells::GetBindingEnergy(iz,k);
+          res += f*ne*LShell(e1/e0,eeff/e1);
+          //res += f*ne*LShell(e1/eshell,eeff/e1);
         }
+      }
         /*
         if(28 >= iz) {
-          term += f*(Z - 10.)*atomDensity[i]*LShell(eshell,HM[iz-11]*eta)/Z; 
+          res += f*(Z - 10.)*LShell(eshell,HM[iz-11]*eta); 
         } else if(32 >= iz) {
-          term += f*18.0*atomDensity[i]*LShell(eshell,HM[iz-11]*eta)/Z;
+          res += f*18.0*LShell(eshell,HM[iz-11]*eta);
         } else if(60 >= iz) {
-          term += f*18.0*atomDensity[i]*LShell(eshell,HM[iz-11]*eta)/Z; 
-          term += f*(Z - 28.)*atomDensity[i]*LShell(eshell,HN[iz-33]*eta)/Z; 
+          res += f*18.0*LShell(eshell,HM[iz-11]*eta); 
+          res += f*(Z - 28.)*LShell(eshell,HN[iz-33]*eta); 
         } else {
-          term += f*18.0*atomDensity[i]*LShell(eshell,HM[53]*eta)/Z;
-          term += f*32.0*atomDensity[i]*LShell(eshell,HN[30]*eta)/Z;
-          term += f*(Z - 60.)*atomDensity[i]*LShell(eshell,150.*eta)/Z; 
+          res += f*18.0*LShell(eshell,HM[50]*eta);
+          res += f*32.0*LShell(eshell,HN[30]*eta);
+          res += f*(Z - 60.)*LShell(eshell,150.*eta); 
         }
         */
-      }
-    }
-    //term += atomDensity[i]*MSH[iz]/(ba2*ba2);
+    }
+    //term += atomDensity[i]*(res/Z + MSH[iz]/(ba2*ba2));
+    term += res*atomDensity[i]/Z;
   }
 
   term /= material->GetTotNbOfAtomsPerVolume();
+  //if(charge < 0.0) { term = -term; }
   return term;
 }
@@ -552 +556 @@
 
   BarkasTerm *= 1.29*charge/material->GetTotNbOfAtomsPerVolume();
+
+  // temporary protection
+  if(charge < -7.0 ) { BarkasTerm *= (-7.0/charge); }
+
   return BarkasTerm;
 }
@@ -574 +582 @@
   } while (del > 0.01*term);
 
-  return -y2*term;
+  G4double res = -y2*term;
+  // temporary protection
+  if(q2 > 49. && res < -0.2) { res = -0.2; }
+
+  return res;
 }
 
@@ -808 +820 @@
   ionList[idx]  = ion;
   stopData[idx] = vv;
-  if(verbose>1) G4cout << "End data set " << G4endl;
+  if(verbose>1) { G4cout << "End data set " << G4endl; }
 }
 
@@ -1395 +1407 @@
     18.2
   };
-  for(i=0; i<53; i++) {HM[i] = hm[i];}
-  for(i=0; i<31; i++) {HN[i] = hn[i];}
+  for(i=0; i<53; ++i) {HM[i] = hm[i];}
+  for(i=0; i<31; ++i) {HN[i] = hn[i];}
 
   const G4double mm[93] = {