Changeset 1196 for trunk/source/processes/electromagnetic/utils/src

Timestamp:

Nov 25, 2009, 5:13:58 PM (15 years ago)

Author:

garnier

Message:

update CVS release candidate geant4.9.3.01

Location:

trunk/source/processes/electromagnetic/utils/src

Files:

• G4DummyModel.cc (modified) (1 diff)
• G4ElectronIonPair.cc (modified) (1 diff)
• G4EmCalculator.cc (modified) (19 diffs)
• G4EmConfigurator.cc (modified) (3 diffs)
• G4EmCorrections.cc (modified) (5 diffs)
• G4EmElementSelector.cc (modified) (6 diffs)
• G4EmModelManager.cc (modified) (30 diffs)
• G4EmMultiModel.cc (modified) (1 diff)
• G4EmProcessOptions.cc (modified) (2 diffs)
• G4EmSaturation.cc (modified) (8 diffs)
• G4EnergyLossMessenger.cc (modified) (8 diffs)
• G4EnergyLossTables.cc (modified) (1 diff)
• G4LossTableBuilder.cc (modified) (9 diffs)
• G4LossTableManager.cc (modified) (3 diffs)
• G4VEmFluctuationModel.cc (modified) (1 diff)
• G4VEmModel.cc (modified) (7 diffs)
• G4VEmProcess.cc (modified) (24 diffs)
• G4VEnergyLoss.cc (modified) (1 diff)
• G4VEnergyLossProcess.cc (modified) (45 diffs)
• G4VMscModel.cc (modified) (2 diffs)
• G4VMultipleScattering.cc (modified) (14 diffs)
• G4ionEffectiveCharge.cc (modified) (2 diffs)

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

@@ -25 +25 @@
 //
 // $Id: G4DummyModel.cc,v 1.4 2009/04/07 18:39:47 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------

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

@@ -25 +25 @@
 //
 // $Id: G4ElectronIonPair.cc,v 1.2 2008/10/17 14:46:16 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4EmCalculator.cc,v 1.46 2009/02/24 09:56:03 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4EmCalculator.cc,v 1.48 2009/11/11 23:59:48 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -127 +127 @@
   if(couple && UpdateParticle(p, kinEnergy) ) {
     res = manager->GetDEDX(p, kinEnergy, couple);
+    /*
     if(isIon) {
-      G4double eth = 2.0*MeV/massRatio;
-      if(kinEnergy > eth) {
-        G4double x1 = corr->ComputeIonCorrections(p,mat,kinEnergy);
-        G4double x2 = corr->ComputeIonCorrections(p,mat,eth);
-        res += x1 - x2*eth/kinEnergy;
-        /*      
-        G4cout << "### GetDEDX: E= " << kinEnergy << " res= " << res 
-               << " x1= " << x1 << " x2= " << x2 
-               << " del= " << x1 - x2*eth/kinEnergy << G4endl;;
-        */
-      } 
-    }
-
+      if(FindEmModel(p, currentProcessName, kinEnergy)) {
+        G4double length = 1.*um;
+        G4cout << "### GetDEDX: E= " << kinEnergy << " res= " << res; 
+        G4double eloss = res*length;
+        G4double niel  = 0.0;
+        currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length);
+        res = eloss/length; 
+     
+        //G4cout << "### GetDEDX: E= " << kinEnergy << " res= " << res;  
+        G4cout << " res1= " << res << G4endl;;
+      }
+    } 
+    */
     if(verbose>0) {
       G4cout << "G4EmCalculator::GetDEDX: E(MeV)= " << kinEnergy/MeV
@@ -147 +148 @@
              << "  " <<  p->GetParticleName()
              << " in " <<  mat->GetName()
+             << " isIon= " << isIon
              << G4endl;
     }
@@ -307 +309 @@
     FindLambdaTable(p, processName);
     if(currentLambda) {
-      G4bool b;
       G4double e = kinEnergy*massRatio;
-      res = (((*currentLambda)[idx])->GetValue(e,b))*chargeSquare;
+      res = (((*currentLambda)[idx])->Value(e))*chargeSquare;
       if(verbose>0) {
         G4cout << "E(MeV)= " << kinEnergy/MeV
@@ -436 +437 @@
       }
 
-      // emulate boundary region for different parameterisations
+      // emulate smoothing procedure
       G4double eth = currentModel->LowEnergyLimit();
       //      G4cout << "massRatio= " << massRatio << " eth= " << eth << G4endl;
-      if(eth > 0.05*MeV && eth < 10.*MeV && escaled > eth && 
-         loweModel != currentModel && loweModel) {
+      if(loweModel) {
         G4double res0 = 0.0;
         G4double res1 = 0.0;
@@ -462 +462 @@
                << res1 <<  "  q2= " << chargeSquare << G4endl;
         */
-        res *= (1.0 + (res0/res1 - 1.0)*eth/escaled);
-
-        if(isIon) {
-          G4double ethscaled = eth/massRatio;
-          if(kinEnergy > ethscaled) {
-            G4double x1 = corr->ComputeIonCorrections(p,mat,kinEnergy);
-            G4double x2 = corr->ComputeIonCorrections(p,mat,ethscaled);
-            res += x1 - x2*ethscaled/kinEnergy;
-          } 
+        res += (res0 - res1)*eth/escaled;
+      } 
+
+      // low energy correction for ions
+      /*
+      if(isIon) {
+        G4double length = 1.*nm;
+        const G4Region* r = 0;
+        const G4MaterialCutsCouple* couple = FindCouple(mat, r);
+        G4double eloss = res*length;
+        G4double niel  = 0.0;
+        currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length);
+        res = eloss/length; 
         
-          if(verbose > 1) {
-            G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV
-                   << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
-                   << G4endl;
-          }
+        if(verbose > 1) {
+          G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV
+                 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
+                 << G4endl;
         }
       }
+      */
+    }
       
-      if(verbose > 0) {
-        G4cout << "E(MeV)= " << kinEnergy/MeV
-               << " DEDX(MeV/mm)= " << res*mm/MeV
-               << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
-               << " cut(MeV)= " << cut/MeV
-               << "  " <<  p->GetParticleName()
-               << " in " <<  currentMaterialName
-               << " Zi^2= " << chargeSquare
-               << G4endl;
-      }
+    if(verbose > 0) {
+      G4cout << "E(MeV)= " << kinEnergy/MeV
+             << " DEDX(MeV/mm)= " << res*mm/MeV
+             << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
+             << " cut(MeV)= " << cut/MeV
+             << "  " <<  p->GetParticleName()
+             << " in " <<  currentMaterialName
+             << " Zi^2= " << chargeSquare
+             << G4endl;
     }
   }
@@ -745 +749 @@
 {
   if(p != currentParticle) {
+
+    // new particle
     currentParticle = p;
+    dynParticle.SetDefinition(const_cast<G4ParticleDefinition*>(p));
+    dynParticle.SetKineticEnergy(kinEnergy);
     baseParticle    = 0;
     currentParticleName = p->GetParticleName();
@@ -753 +761 @@
     currentProcess  = FindEnergyLossProcess(p);
     currentProcessName = "";
-    if(currentProcess) currentProcessName = currentProcess->GetProcessName();
-
-    if(p->GetParticleType() == "nucleus" 
-       && currentParticleName != "deuteron"  
-       && currentParticleName != "triton"
-       && currentParticleName != "alpha+"
-       && currentParticleName != "helium"
-       && currentParticleName != "hydrogen"
-      ) {
-      baseParticle = theGenericIon;
-      massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
-      isIon = true;
-      //      G4cout << p->GetParticleName()
-      // << " in " << currentMaterial->GetName()
-      //       << "  e= " << kinEnergy << G4endl;
-      chargeSquare =
-        ionEffCharge->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy);
-      //G4cout << "q2= " << chargeSquare << G4endl;
-    } else {
-      isIon = false;
-      if(currentProcess) {
-        baseParticle    = currentProcess->BaseParticle();
-
-        if(baseParticle) {
-          massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
-          G4double q = baseParticle->GetPDGCharge()/eplus;
-          chargeSquare /= (q*q);
-        }
-      }
-    }
-  }
+    isIon = false;
+
+    // ionisation process exist
+    if(currentProcess) {
+      currentProcessName = currentProcess->GetProcessName();
+      baseParticle = currentProcess->BaseParticle();
+
+      // base particle is used
+      if(baseParticle) {
+        massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
+        G4double q = p->GetPDGCharge()/baseParticle->GetPDGCharge();
+        chargeSquare = q*q;
+      } 
+
+      if(p->GetParticleType()   == "nucleus" 
+         && currentParticleName != "deuteron"  
+         && currentParticleName != "triton"
+         && currentParticleName != "alpha+"
+         && currentParticleName != "helium"
+         && currentParticleName != "hydrogen"
+         ) {
+        isIon = true;
+        massRatio = theGenericIon->GetPDGMass()/p->GetPDGMass();
+        baseParticle = theGenericIon;
+        //      G4cout << p->GetParticleName()
+        // << " in " << currentMaterial->GetName()
+        //       << "  e= " << kinEnergy << G4endl;
+      }
+    }
+  }
+
+  // Effective charge for ions
   if(isIon) {
     chargeSquare =
-     ionEffCharge->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy)
+      corr->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy)
       * corr->EffectiveChargeCorrection(p,currentMaterial,kinEnergy);
     if(currentProcess) {
@@ -810 +819 @@
     p = currentParticle;
   }
+  return p;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+const G4ParticleDefinition* G4EmCalculator::FindIon(G4int Z, G4int A)
+{
+  const G4ParticleDefinition* p = 
+    G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
   return p;
 }
@@ -897 +915 @@
     G4String partname =  p->GetParticleName();
     const G4ParticleDefinition* part = p;
-    if(isIon) part = theGenericIon;
-
+    if(isIon) { part = theGenericIon; }
+
+    // energy loss process
     G4LossTableManager* lManager = G4LossTableManager::Instance();
     const std::vector<G4VEnergyLossProcess*> vel = 
-      lManager->GetEnergyLossProcessVector();
+    lManager->GetEnergyLossProcessVector();
     G4int n = vel.size();
     for(G4int i=0; i<n; i++) {
       if((vel[i])->GetProcessName() == currentName && 
          (vel[i])->Particle() == part) 
-      {
-        currentLambda = (vel[i])->LambdaTable();
-        isApplicable    = true;
-        break;
-      }
-    }
+        {
+          currentLambda = (vel[i])->LambdaTable();
+          isApplicable    = true;
+          break;
+        }
+    }
+  
+    // discrete process
     if(!currentLambda) {
       const std::vector<G4VEmProcess*> vem = lManager->GetEmProcessVector();
@@ -925 +946 @@
       }
     }
+
+    // msc process
     if(!currentLambda) {
       const std::vector<G4VMultipleScattering*> vmsc = 
@@ -957 +980 @@
   const G4ParticleDefinition* part = p;
   G4double scaledEnergy = kinEnergy*massRatio;
-  if(isIon) part = theGenericIon; 
+  if(isIon) { part = theGenericIon; } 
 
   if(verbose > 1) {
@@ -967 +990 @@
   }
 
-  // Search for the process
+  // Search for energy loss process
   currentName = processName;
   currentModel = 0;
@@ -976 +999 @@
     lManager->GetEnergyLossProcessVector();
   G4int n = vel.size();
+  G4VEnergyLossProcess* elproc = 0;
   for(G4int i=0; i<n; i++) {
     //    G4cout << "i= " << i << " part= " 
     //  << (vel[i])->Particle()->GetParticleName()
     //     << "   proc= " << (vel[i])->GetProcessName()  << G4endl;
-    if((vel[i])->GetProcessName() == currentName && 
-       (vel[i])->Particle() == part) 
-    {
-      const G4ParticleDefinition* bp = (vel[i])->BaseParticle();
-      //      G4cout << "i= " << i << " bp= " << bp << G4endl;
-      if(!bp) {
-        currentModel = (vel[i])->SelectModelForMaterial(scaledEnergy, idx);
-        loweModel = (vel[i])->SelectModelForMaterial(keV, idx);
-        isApplicable    = true;
-        break;
+    if((vel[i])->GetProcessName() == currentName) {
+      if(baseParticle) {
+        if((vel[i])->Particle() == baseParticle) {
+          elproc = vel[i];
+          break;
+        }
       } else {
-        for(G4int j=0; j<n; j++) {
-          if((vel[j])->Particle() == bp) {
-            currentModel = (vel[j])->SelectModelForMaterial(scaledEnergy, idx);
-            loweModel = (vel[j])->SelectModelForMaterial(keV, idx);
-            isApplicable    = true;
-            break;
-          }
+        if((vel[i])->Particle() == part) {
+          elproc = vel[i];
+          break;
         }
       }
     }
   }
+  if(elproc) {
+    isApplicable = true;
+    currentModel = elproc->SelectModelForMaterial(scaledEnergy, idx);
+    G4double eth = currentModel->LowEnergyLimit();
+    if(eth > 0.0) { loweModel = elproc->SelectModelForMaterial(eth-keV, idx); }       
+  }
+
+  // Search for discrete process
   if(!currentModel) {
     const std::vector<G4VEmProcess*> vem = lManager->GetEmProcessVector();
@@ -1010 +1034 @@
       {
         currentModel = (vem[i])->SelectModelForMaterial(kinEnergy, idx);
-        loweModel = (vem[i])->SelectModelForMaterial(keV, idx);
+        G4double eth = currentModel->LowEnergyLimit();
+        if(eth > 0.0) { loweModel = (vem[i])->SelectModelForMaterial(eth-keV, idx); }       
         isApplicable    = true;
         break;
@@ -1016 +1041 @@
     }
   }
+
+  // Search for msc process
   if(!currentModel) {
     const std::vector<G4VMultipleScattering*> vmsc = 
@@ -1025 +1052 @@
       {
         currentModel = (vmsc[i])->SelectModelForMaterial(kinEnergy, idx);
-        loweModel = (vmsc[i])->SelectModelForMaterial(keV, idx);
+        G4double eth = currentModel->LowEnergyLimit();
+        if(eth > 0.0) { loweModel = (vmsc[i])->SelectModelForMaterial(eth-keV, idx); }       
         isApplicable    = true;
         break;
@@ -1043 +1071 @@
   G4String partname =  p->GetParticleName();
   const G4ParticleDefinition* part = p;
+  
   if(p->GetParticleType() == "nucleus" && 
      partname != "deuteron" && 
-     partname != "triton") part = G4GenericIon::GenericIon(); 
+     partname != "triton") { part = theGenericIon; } 
   
   G4LossTableManager* lManager = G4LossTableManager::Instance();
@@ -1052 +1081 @@
   G4int n = vel.size();
   for(G4int i=0; i<n; i++) {
-    if((vel[i])->Particle() == part) {
+    if( (vel[i])->Particle() == part ) {
       elp = vel[i];
       break;

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4EmConfigurator.cc,v 1.4 2009/02/26 11:33:33 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4EmConfigurator.cc,v 1.5 2009/07/20 17:06:48 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -223 +223 @@
           // set fluctuation model for ionisation processes
           if(fluc && ptype == eloss) {
-            G4VEnergyLossProcess* p = reinterpret_cast<G4VEnergyLossProcess*>(proc);
+            G4VEnergyLossProcess* p = static_cast<G4VEnergyLossProcess*>(proc);
             p->SetFluctModel(fluc);
          
@@ -262 +262 @@
           // added model
           if(ptype == eloss) {
-            G4VEnergyLossProcess* p = reinterpret_cast<G4VEnergyLossProcess*>(proc);
+            G4VEnergyLossProcess* p = static_cast<G4VEnergyLossProcess*>(proc);
             p->AddEmModel(index,mod,fluc,reg);
             //G4cout << "### Added eloss model order= " << index << " for " 
             //     << particleName << " and " << processName << "  " << mod << G4endl;
           } else if(ptype == discrete) {
-            G4VEmProcess* p = reinterpret_cast<G4VEmProcess*>(proc);
+            G4VEmProcess* p = static_cast<G4VEmProcess*>(proc);
             p->AddEmModel(index,mod,reg);
           } else if(ptype == msc) {
-            G4VMultipleScattering* p = reinterpret_cast<G4VMultipleScattering*>(proc);
+            G4VMultipleScattering* p = static_cast<G4VMultipleScattering*>(proc);
             p->AddEmModel(index,mod,reg);
           }

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4EmCorrections.cc,v 1.51 2008/12/18 13:01:44 gunter Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4EmCorrections.cc,v 1.54 2009/10/29 17:56:36 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -687 +687 @@
     massFactor = proton_mass_c2/p->GetPDGMass();
     idx = -1;
-    G4int dz = 1000;
 
     for(G4int i=0; i<nIons; i++) {
-      if(materialList[i] == mat) {
-        G4int delz = currentZ - Zion[i];
-        if(delz < 0) delz = -delz;
-        if(delz < dz) {
-          idx = i;
-          dz = delz;
-          if(0 == delz) break;
-        }
+      if(materialList[i] == mat && currentZ == Zion[i]) {
+        idx = i;
+        break;
       }
     }
-    //    G4cout << " idx= " << idx << " dz= " << dz << G4endl;
-    if(idx > 0) {
+    //    G4cout << " idx= " << idx << " dz= " << G4endl;
+    if(idx >= 0) {
       if(!ionList[idx]) BuildCorrectionVector(); 
       if(ionList[idx])  curVector = stopData[idx];
-    }
+    } else { return factor; }
   }
   if(curVector) {
-    G4bool b;
-    factor = curVector->GetValue(ekin*massFactor,b);
+    factor = curVector->Value(ekin*massFactor);
     if(verbose > 1) {
       G4cout << "E= " << ekin << " factor= " << factor << " massfactor= " 
@@ -771 +764 @@
            << " massRatio= " << massRatio << G4endl;
   }
-  G4bool b;
 
   G4PhysicsLogVector* vv = 
@@ -777 +769 @@
   vv->SetSpline(true);
   G4double e, eion, dedx, dedx1;
-  G4double eth0 = v->GetLowEdgeEnergy(0);
+  G4double eth0 = v->Energy(0);
   G4double escal = eth/massRatio;
   G4double qe = 
@@ -790 +782 @@
   //     << " dedxt1= " << dedx1t << G4endl;   
 
-  for(G4int i=0; i<nbinCorr; i++) {
-    e = vv->GetLowEdgeEnergy(i);
+  for(G4int i=0; i<=nbinCorr; i++) {
+    e = vv->Energy(i);
     escal = e/massRatio;
     eion  = escal/A;
     if(eion <= eth0) {
-      dedx = v->GetValue(eth0, b)*std::sqrt(eion/eth0);
+      dedx = v->Value(eth0)*std::sqrt(eion/eth0);
     } else {
-      dedx = v->GetValue(eion, b);
+      dedx = v->Value(eion);
     }
     qe = effCharge.EffectiveChargeSquareRatio(curParticle,curMaterial,escal); 

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4EmElementSelector.cc,v 1.10 2009/05/26 16:59:35 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4EmElementSelector.cc,v 1.11 2009/09/29 11:31:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -68 +68 @@
   if(nElmMinusOne > 0) {
     xSections.reserve(n);
-    for(G4int i=0; i<n; i++) {
+    for(G4int i=0; i<n; ++i) {
       G4PhysicsLogVector* v = new G4PhysicsLogVector(lowEnergy,highEnergy,nbins);
       //v->SetSpline(spline);
@@ -82 +82 @@
 {
   if(nElmMinusOne > 0) {
-    for(G4int i=0; i<=nElmMinusOne; i++) {
+    for(G4int i=0; i<=nElmMinusOne; ++i) {
       delete xSections[i];
     }
@@ -105 +105 @@
 
   // loop over bins
-  for(G4int j=0; j<nbins; j++) {
-    G4double e = (xSections[0])->GetLowEdgeEnergy(j);
+  for(G4int j=0; j<=nbins; ++j) {
+    G4double e = (xSections[0])->Energy(j);
     model->SetupForMaterial(part, material, e);
     cross = 0.0;
     //G4cout << "j= " << j << " e(MeV)= " << e/MeV << G4endl;
-    for (i=0; i<=nElmMinusOne; i++) {
+    for (i=0; i<=nElmMinusOne; ++i) {
       cross += theAtomNumDensityVector[i]*      
         model->ComputeCrossSectionPerAtom(part, (*theElementVector)[i], e, 
@@ -120 +120 @@
   // xSections start from null, so use probabilities from the next bin
   if(DBL_MIN >= (*xSections[nElmMinusOne])[0]) {
-    for (i=0; i<=nElmMinusOne; i++) {
+    for (i=0; i<=nElmMinusOne; ++i) {
       xSections[i]->PutValue(0, (*xSections[i])[1]);
     }
   }
   // xSections ends with null, so use probabilities from the previous bin
-  if(DBL_MIN >= (*xSections[nElmMinusOne])[nbins-1]) {
-    for (i=0; i<=nElmMinusOne; i++) {
-      xSections[i]->PutValue(nbins-1, (*xSections[i])[nbins-2]);
+  if(DBL_MIN >= (*xSections[nElmMinusOne])[nbins]) {
+    for (i=0; i<=nElmMinusOne; ++i) {
+      xSections[i]->PutValue(nbins, (*xSections[i])[nbins-1]);
     }
   }
   // perform normalization
-  for(G4int j=0; j<nbins; j++) {
+  for(G4int j=0; j<=nbins; ++j) {
     cross = (*xSections[nElmMinusOne])[j];
     // only for positive X-section 
     if(cross > DBL_MIN) {
-      for (i=0; i<nElmMinusOne; i++) {
-        xSections[i]->PutValue(j, (*xSections[i])[j]/cross);
+      for (i=0; i<nElmMinusOne; ++i) {
+        G4double x = (*xSections[i])[j]/cross;
+        xSections[i]->PutValue(j, x);
       }
     }
@@ -158 +159 @@
          << (*theElementVector)[nElmMinusOne]->GetName() 
          << G4endl;
+  G4cout << G4endl;
 }
 

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4EmModelManager.cc,v 1.49 2009/04/17 10:35:32 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4EmModelManager.cc,v 1.58 2009/10/29 18:07:08 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -62 +62 @@
 // 12-04-07 Add verbosity at destruction (V.Ivanchenko)
 // 08-04-08 Fixed and simplified initialisation of G4RegionModel (VI)
+// 03-08-09 Create internal vectors only it is needed (VI)
 //
 // Class Description:
@@ -85 +86 @@
   theListOfModelIndexes = new G4int [nModelsForRegion];
   lowKineticEnergy      = new G4double [nModelsForRegion+1];
-  for (G4int i=0; i<nModelsForRegion; i++) {
+  for (G4int i=0; i<nModelsForRegion; ++i) {
     theListOfModelIndexes[i] = indx[i];
     lowKineticEnergy[i] = lowE[i];
@@ -114 +115 @@
 #include "G4Positron.hh"
 #include "G4UnitsTable.hh"
+#include "G4DataVector.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -120 +122 @@
   nEmModels(0),
   nRegions(0),
-  nCouples(0),
-  idxOfRegionModels(0),
-  setOfRegionModels(0),
-  minSubRange(0.1),
   particle(0),
   verboseLevel(0)
 {
-  models.clear();
-  flucModels.clear();
-  regions.clear();
-  orderOfModels.clear();
-  maxCutInRange    = 12.*cm;
   maxSubCutInRange = 0.7*mm;
-  theGamma = G4Gamma::Gamma();
-  thePositron = G4Positron::Positron();
   models.reserve(4);
   flucModels.reserve(4);
@@ -140 +131 @@
   orderOfModels.reserve(4);
   isUsed.reserve(4);
+  severalModels = true;
+  currRegionModel = 0;
+  currModel = 0;
 }
 
@@ -146 +140 @@
 G4EmModelManager::~G4EmModelManager()
 {
+  verboseLevel = 0; // no verbosity at destruction
   Clear();
 }
@@ -156 +151 @@
     G4cout << "G4EmModelManager::Clear()" << G4endl;
   }
-
-  theCuts.clear();
-  theSubCuts.clear();
-  if(idxOfRegionModels) delete [] idxOfRegionModels;
-  if(setOfRegionModels && nRegions) {
-    for(G4int i=0; i<nRegions; i++) {
-      delete (setOfRegionModels[i]);
-    }
-    delete [] setOfRegionModels;
-  }
-  idxOfRegionModels = 0;
-  setOfRegionModels = 0;
+  size_t n = setOfRegionModels.size();
+  if(n > 0) {
+    for(size_t i=0; i<n; ++i) {
+      delete setOfRegionModels[i];
+      setOfRegionModels[i] = 0;
+    }
+  }
 }
 
@@ -195 +185 @@
 {
   if (nEmModels > 0) {
-    for(G4int i=0; i<nEmModels; i++) {
+    for(G4int i=0; i<nEmModels; ++i) {
       if(nam == models[i]->GetName()) {
         models[i]->SetLowEnergyLimit(emin);
@@ -226 +216 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-const G4DataVector* G4EmModelManager::Initialise(const G4ParticleDefinition* p,
-                                                 const G4ParticleDefinition* sp,
-                                                 G4double theMinSubRange,
-                                                 G4int val)
+const G4DataVector* 
+G4EmModelManager::Initialise(const G4ParticleDefinition* p,
+                             const G4ParticleDefinition* secondaryParticle,
+                             G4double minSubRange,
+                             G4int val)
 {
   verboseLevel = val;
@@ -238 +229 @@
   }
   // Are models defined?
-  if(!nEmModels) {
-    G4Exception("G4EmModelManager::Initialise without any model defined for "+partname);
-  }
+  if(nEmModels < 1) {
+    G4Exception("G4EmModelManager::Initialise: no model defined for " + partname);
+  }
+
   particle = p;
-  secondaryParticle = sp;
-  minSubRange = theMinSubRange;
-
-  Clear();
+  Clear(); // needed if run is not first
+
   G4RegionStore* regionStore = G4RegionStore::GetInstance();
   const G4Region* world = 
@@ -256 +246 @@
   G4bool isWorld = false;
 
-  for (G4int ii=0; ii<nEmModels; ii++) {
+  for (G4int ii=0; ii<nEmModels; ++ii) {
     const G4Region* r = regions[ii];
     if ( r == 0 || r == world) {
@@ -264 +254 @@
       G4bool newRegion = true;
       if (nRegions>1) {
-        for (G4int j=1; j<nRegions; j++) {
+        for (G4int j=1; j<nRegions; ++j) {
           if ( r == setr[j] ) newRegion = false;
         }
@@ -274 +264 @@
     }
   }
+  // Are models defined?
+  if(!isWorld) {
+    G4Exception("G4EmModelManager::Initialise: no models defined for " + 
+                partname + " in the World volume");
+  }
 
   G4ProductionCutsTable* theCoupleTable=
     G4ProductionCutsTable::GetProductionCutsTable();
-  G4int numOfCouples = theCoupleTable->GetTableSize();
-  if(nRegions > 1) idxOfRegionModels = new G4int[numOfCouples];
-  setOfRegionModels = new G4RegionModels*[nRegions];
+  size_t numOfCouples = theCoupleTable->GetTableSize();
+
+  // prepare vectors, shortcut for the case of only 1 model
+  if(nRegions > 1 && nEmModels > 1) {
+    if(numOfCouples > idxOfRegionModels.size()) idxOfRegionModels.resize(numOfCouples);
+  }
+  size_t nr = 1;
+  if(nEmModels > 1) nr = nRegions;
+  if(nr > setOfRegionModels.size()) setOfRegionModels.resize(nr);
 
   std::vector<G4int>    modelAtRegion(nEmModels);
@@ -287 +288 @@
 
   // Order models for regions
-  for (G4int reg=0; reg<nRegions; reg++) {
+  for (G4int reg=0; reg<nRegions; ++reg) {
     const G4Region* region = setr[reg];
     G4int n = 0;
 
-    if(isWorld || 0 < reg) {
-
-      for (G4int ii=0; ii<nEmModels; ii++) {
-
-        G4VEmModel* model = models[ii];
-        if ( region == regions[ii] ) {
-
-          G4double tmin = model->LowEnergyLimit();
-          G4double tmax = model->HighEnergyLimit();
-          G4int  ord    = orderOfModels[ii];
-          G4bool push   = true;
-          G4bool insert = false;
-          G4int  idx    = n;
-
-          if(1 < verboseLevel) {
-            G4cout << "Model #" << ii 
-                   << " <" << model->GetName() << "> for region <";
-            if (region) G4cout << region->GetName();
-            G4cout << ">  "
-                   << " tmin(MeV)= " << tmin/MeV
-                   << "; tmax(MeV)= " << tmax/MeV
-                   << "; order= " << ord
-                   << G4endl;
-          }
+    for (G4int ii=0; ii<nEmModels; ++ii) {
+
+      G4VEmModel* model = models[ii];
+      if ( region == regions[ii] ) {
+
+        G4double tmin = model->LowEnergyLimit();
+        G4double tmax = model->HighEnergyLimit();
+        G4int  ord    = orderOfModels[ii];
+        G4bool push   = true;
+        G4bool insert = false;
+        G4int  idx    = n;
+
+        if(1 < verboseLevel) {
+          G4cout << "Model #" << ii 
+                 << " <" << model->GetName() << "> for region <";
+          if (region) G4cout << region->GetName();
+          G4cout << ">  "
+                 << " tmin(MeV)= " << tmin/MeV
+                 << "; tmax(MeV)= " << tmax/MeV
+                 << "; order= " << ord
+                 << G4endl;
+        }
         
-          if(n > 0) {
-
-            // extend energy range to previous models
-            tmin = std::min(tmin, eHigh[n-1]);
-            tmax = std::max(tmax, eLow[0]);
-            //G4cout << "tmin= " << tmin << "  tmax= " 
-            //     << tmax << "  ord= " << ord <<G4endl;
-            // empty energy range
-            if( tmax - tmin <= eV) push = false;
-            // low-energy model
-            else if (tmax == eLow[0]) {
-              push = false;
-              insert = true;
-              idx = 0;
-              // resolve intersections
-            } else if(tmin < eHigh[n-1]) { 
-              // compare order
-              for(G4int k=0; k<n; k++) {
-                // new model has lower application 
-                if(ord >= modelOrd[k]) {
-                  if(tmin < eHigh[k]  && tmin >= eLow[k]) tmin = eHigh[k];
-                  if(tmax <= eHigh[k] && tmax >  eLow[k]) tmax = eLow[k];
-                  if(tmax > eHigh[k] && tmin < eLow[k]) {
-                    if(tmax - eHigh[k] > eLow[k] - tmin) tmin = eHigh[k];
-                    else tmax = eLow[k];
-                  }
-                  if( tmax - tmin <= eV) {
-                    push = false;
-                    break;
-                  }
+        if(n > 0) {
+
+          // extend energy range to previous models
+          tmin = std::min(tmin, eHigh[n-1]);
+          tmax = std::max(tmax, eLow[0]);
+          //G4cout << "tmin= " << tmin << "  tmax= " 
+          //       << tmax << "  ord= " << ord <<G4endl;
+          // empty energy range
+          if( tmax - tmin <= eV) push = false;
+          // low-energy model
+          else if (tmax == eLow[0]) {
+            push = false;
+            insert = true;
+            idx = 0;
+            // resolve intersections
+          } else if(tmin < eHigh[n-1]) { 
+            // compare order
+            for(G4int k=0; k<n; ++k) {
+              // new model has lower application 
+              if(ord >= modelOrd[k]) {
+                if(tmin < eHigh[k]  && tmin >= eLow[k]) tmin = eHigh[k];
+                if(tmax <= eHigh[k] && tmax >  eLow[k]) tmax = eLow[k];
+                if(tmax > eHigh[k] && tmin < eLow[k]) {
+                  if(tmax - eHigh[k] > eLow[k] - tmin) tmin = eHigh[k];
+                  else tmax = eLow[k];
+                }
+                if( tmax - tmin <= eV) {
+                  push = false;
+                  break;
                 }
               }
-              //G4cout << "tmin= " << tmin << "  tmax= " 
-              //     << tmax << "  push= " << push << " idx= " << idx <<G4endl;
-              if(push) {
-                if (tmax == eLow[0]) {
+            }
+            //G4cout << "tmin= " << tmin << "  tmax= " 
+            //     << tmax << "  push= " << push << " idx= " << idx <<G4endl;
+            if(push) {
+              if (tmax == eLow[0]) {
+                push = false;
+                insert = true;
+                idx = 0;
+                // continue resolve intersections
+              } else if(tmin < eHigh[n-1]) { 
+                // last energy interval
+                if(tmin > eLow[n-1] && tmax >= eHigh[n-1]) {
+                  eHigh[n-1] = tmin;
+                  // first energy interval
+                } else if(tmin <= eLow[0] && tmax < eHigh[0]) {
+                  eLow[0] = tmax;
                   push = false;
                   insert = true;
                   idx = 0;
-                  // continue resolve intersections
-                } else if(tmin < eHigh[n-1]) { 
-                  // last energy interval
-                  if(tmin > eLow[n-1] && tmax >= eHigh[n-1]) {
-                    eHigh[n-1] = tmin;
-                    // first energy interval
-                  } else if(tmin <= eLow[0] && tmax < eHigh[0]) {
-                    eLow[0] = tmax;
-                    push = false;
-                    insert = true;
-                    idx = 0;
-                  } else {
-                    // find energy interval to replace
-                    for(G4int k=0; k<n; k++) { 
-                      if(tmin <= eLow[k] && tmax >= eHigh[k]) {
-                        push = false;
-                        modelAtRegion[k] = ii;
-                        modelOrd[k] = ord;
-                        isUsed[ii] = 1;
-                      } 
-                    }
+                } else {
+                  // find energy interval to replace
+                  for(G4int k=0; k<n; ++k) { 
+                    if(tmin <= eLow[k] && tmax >= eHigh[k]) {
+                      push = false;
+                      modelAtRegion[k] = ii;
+                      modelOrd[k] = ord;
+                      isUsed[ii] = 1;
+                    } 
                   }
                 }
@@ -381 +379 @@
             }
           }
-          if(insert) {
-            for(G4int k=n-1; k>=idx; k--) {            
-              modelAtRegion[k+1] = modelAtRegion[k];
-              modelOrd[k+1] = modelOrd[k];
-              eLow[k+1]  = eLow[k];
-              eHigh[k+1] = eHigh[k];
-            }
-          }
-          //G4cout << "push= " << push << " insert= " << insert 
-          //<< " idx= " << idx <<G4endl;
-          if (push || insert) {
-            n++;
-            modelAtRegion[idx] = ii;
-            modelOrd[idx] = ord;
-            eLow[idx]  = tmin;
-            eHigh[idx] = tmax;
-            isUsed[ii] = 1;
+        }
+        if(insert) {
+          for(G4int k=n-1; k>=idx; --k) {    
+            modelAtRegion[k+1] = modelAtRegion[k];
+            modelOrd[k+1] = modelOrd[k];
+            eLow[k+1]  = eLow[k];
+            eHigh[k+1] = eHigh[k];
           }
         }
-      }
-    } else {
-      n = 1;
-      models.push_back(0);
-      modelAtRegion.push_back(nEmModels);
-      eLow.push_back(0.0);
-      eHigh.push_back(DBL_MAX);
+        //G4cout << "push= " << push << " insert= " << insert 
+        //<< " idx= " << idx <<G4endl;
+        if (push || insert) {
+          ++n;
+          modelAtRegion[idx] = ii;
+          modelOrd[idx] = ord;
+          eLow[idx]  = tmin;
+          eHigh[idx] = tmax;
+          isUsed[ii] = 1;
+        }
+      }
     }
     eLow[0] = 0.0;
@@ -415 +407 @@
       if (region) G4cout << region->GetName();
       G4cout << ">  Elow(MeV)= ";
-      for(G4int ii=0; ii<=n; ii++) {G4cout << eLow[ii]/MeV << " ";}
+      for(G4int ii=0; ii<=n; ++ii) {G4cout << eLow[ii]/MeV << " ";}
       G4cout << G4endl;
     }
     G4RegionModels* rm = new G4RegionModels(n, modelAtRegion, eLow, region);
     setOfRegionModels[reg] = rm;
+    if(1 == nEmModels) break;
   }
 
@@ -425 +418 @@
 
   // Access to materials and build cuts
-  for(G4int i=0; i<numOfCouples; i++) {
+  size_t idx = 1;
+  if(secondaryParticle) {
+    if( secondaryParticle == G4Gamma::Gamma() )           idx = 0;
+    else if( secondaryParticle == G4Positron::Positron()) idx = 2;
+  }
+
+  if(numOfCouples > theCuts.size()) {theCuts.resize(numOfCouples);}
+  if(minSubRange < 1.0 && numOfCouples > theSubCuts.size()) {
+    theSubCuts.resize(numOfCouples);
+  }
+  for(size_t i=0; i<numOfCouples; ++i) {
 
     const G4MaterialCutsCouple* couple = 
@@ -433 +436 @@
  
     G4int reg = 0;
-    if(nRegions > 1) {
+    if(nRegions > 1 && nEmModels > 1) {
       reg = nRegions;
-      do {reg--;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts()));
+      do {--reg;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts()));
       idxOfRegionModels[i] = reg;
     }
@@ -446 +449 @@
     }
 
-    G4double cut = DBL_MAX;
+    G4double cut = (*theCoupleTable->GetEnergyCutsVector(idx))[i]; 
     G4double subcut = DBL_MAX;
     if(secondaryParticle) {
-      size_t idx = 1;
-      if( secondaryParticle == theGamma ) idx = 0;
-      cut = (*theCoupleTable->GetEnergyCutsVector(idx))[i];
-      if( secondaryParticle == thePositron && cut < DBL_MAX )
-        cut += (*theCoupleTable->GetEnergyCutsVector(2))[i] + 
-               2.0*electron_mass_c2;
 
       // compute subcut 
-      if( cut < DBL_MAX ) subcut = minSubRange*cut;
-      if(pcuts->GetProductionCut(idx) < maxCutInRange) {
+      if( cut < DBL_MAX && minSubRange < 1.0) {
+        subcut = minSubRange*cut;
+        G4double rcut = std::min(minSubRange*pcuts->GetProductionCut(idx), 
+                                 maxSubCutInRange);
         G4double tcutmax = 
-          theCoupleTable->ConvertRangeToEnergy(secondaryParticle,
-                                               material,maxSubCutInRange);
+          theCoupleTable->ConvertRangeToEnergy(secondaryParticle,material,rcut);
         if(tcutmax < subcut) subcut = tcutmax;
       }
@@ -467 +465 @@
 
     G4int nm = setOfRegionModels[reg]->NumberOfModels();
-    for(G4int j=0; j<nm; j++) {
+    for(G4int j=0; j<nm; ++j) {
 
       G4VEmModel* model = models[setOfRegionModels[reg]->ModelIndex(j)];
@@ -484 +482 @@
       }
     }
-    theCuts.push_back(cut);
-    theSubCuts.push_back(subcut);
-  }
-
-  for(G4int jj=0; jj<nEmModels; jj++) {
+    theCuts[i] = cut;
+    if(minSubRange < 1.0) theSubCuts[i] = subcut;
+  }
+
+  // initialize models
+  G4int nn = 0;
+  severalModels = true;
+  for(G4int jj=0; jj<nEmModels; ++jj) {
     if(1 == isUsed[jj]) {
-      models[jj]->Initialise(particle, theCuts);
+      ++nn;
+      currModel = models[jj];
+      currModel->Initialise(particle, theCuts);
       if(flucModels[jj]) flucModels[jj]->InitialiseMe(particle);
     }
   }
+  if(1 == nn) severalModels = false;
 
   if(1 < verboseLevel) {
@@ -510 +514 @@
                                       G4EmTableType tType)
 {
-
-  G4double e;
-
   size_t i = couple->GetIndex();
   G4double cut = theCuts[i];
-  G4double subcut = 0.0;
+  G4double emin = 0.0;
 
   if(fTotal == tType) cut = DBL_MAX;
-  else if(fSubRestricted == tType) subcut = theSubCuts[i];
+  else if(fSubRestricted == tType) {
+    emin = cut;
+    if(theSubCuts.size() > 0) emin = theSubCuts[i];
+  }
 
   if(1 < verboseLevel) {
@@ -524 +528 @@
            << couple->GetMaterial()->GetName()
            << "  cut(MeV)= " << cut
-           << "  subcut(MeV)= " << subcut
+           << "  emin(MeV)= " << emin
            << "  Type " << tType
            << "  for " << particle->GetParticleName()
@@ -531 +535 @@
 
   G4int reg  = 0;
-  if(nRegions > 1) reg = idxOfRegionModels[i];
+  if(nRegions > 1 && nEmModels > 1) reg = idxOfRegionModels[i];
   const G4RegionModels* regModels = setOfRegionModels[reg];
   G4int nmod = regModels->NumberOfModels();
 
-  // vectors to provide continues dE/dx
-  G4DataVector factor(nmod);
-  G4DataVector eLow(nmod+1);
-  G4DataVector dedxLow(nmod);
-  G4DataVector dedxHigh(nmod);
-
-  if(1 < verboseLevel) {
-      G4cout << "There are " << nmod << " models for "
-             << couple->GetMaterial()->GetName()
-             << " at the region #" << reg
-             << G4endl;
-  }
-
-  // calculate factors to provide continuity of energy loss
-
-
-  factor[0] = 1.0;
-  G4int j;
-
-  G4int totBinsLoss = aVector->GetVectorLength();
-
-  dedxLow[0] = 0.0;
-  eLow[0]    = 0.0;
-
-  e = regModels->LowEdgeEnergy(1);
-  eLow[1]    = e;
-  G4VEmModel* model = models[regModels->ModelIndex(0)]; 
-  dedxHigh[0] = 0.0;
-  if(model && cut > subcut) {
-    dedxHigh[0] = model->ComputeDEDX(couple,particle,e,cut);
-    if(subcut > 0.0) {
-      dedxHigh[0] -= model->ComputeDEDX(couple,particle,e,subcut);
-    }
-  }
-  if(nmod > 1) {
-    for(j=1; j<nmod; j++) {
-
-      e = regModels->LowEdgeEnergy(j);
-      eLow[j]   = e;
-      G4int idx = regModels->ModelIndex(j); 
-
-      dedxLow[j] = models[idx]->ComputeDEDX(couple,particle,e,cut);
-      if(subcut > 0.0) {
-        dedxLow[j] -= models[idx]->ComputeDEDX(couple,particle,e,subcut);
-      }
-      if(subcut == cut) dedxLow[j] = 0.0;
-
-      e = regModels->LowEdgeEnergy(j+1);
-      eLow[j+1] = e;
-      dedxHigh[j] = models[idx]->ComputeDEDX(couple,particle,e,cut);
-      if(subcut > 0.0) {
-        dedxHigh[j] -= models[idx]->ComputeDEDX(couple,particle,e,subcut);
-      }
-      if(subcut == cut) dedxHigh[j] = 0.0;
-    }
-    if(1 < verboseLevel) {
-      G4cout << " model #0"  
-             << "  dedx(" << eLow[0] << ")=  " << dedxLow[0]
-             << "  dedx(" << eLow[1] << ")=  " << dedxHigh[0]
-             << G4endl;
-    }
-
-    for(j=1; j<nmod; j++) {
-      if(dedxLow[j] > 0.0) {
-        factor[j] = (dedxHigh[j-1]/dedxLow[j] - 1.0)*eLow[j];
-      } else  factor[j] = 0.0;
-      if(1 < verboseLevel) {
-        G4cout << " model #" << j 
-               << "  dedx(" << eLow[j] << ")=  " << dedxLow[j]
-               << "  dedx(" << eLow[j+1] << ")=  " << dedxHigh[j]
-               << "  factor= " << factor[j]/eLow[j]
-               << G4endl;
-      }
-    }
-
-    if(2 < verboseLevel) {
-      G4cout << "Loop over " << totBinsLoss << " bins start " << G4endl;
-    }
-  }
-
   // Calculate energy losses vector
-  for(j=0; j<totBinsLoss; j++) {
-
-    G4double e = aVector->GetLowEdgeEnergy(j);
-    G4double fac = 1.0;
-
-      // Choose a model of energy losses
+
+  //G4cout << "nmod= " << nmod << G4endl;
+  size_t totBinsLoss = aVector->GetVectorLength();
+  for(size_t j=0; j<totBinsLoss; ++j) {
+
+    G4double e = aVector->Energy(j);
+    G4double del = 0.0;
+
+    // Choose a model of energy losses
     G4int k = 0;
-    if (nmod > 1 && e > eLow[1]) {
-      do {
-        k++;
-        fac *= (1.0 + factor[k]/e);
-      } while (k+1 < nmod && e > eLow[k+1]);
-    }
-
-    model = models[regModels->ModelIndex(k)];
-    G4double dedx = 0.0;
-    G4double dedx0 = 0.0;
-    if(model && cut > subcut) {
-      dedx = model->ComputeDEDX(couple,particle,e,cut); 
-      dedx0 = dedx;
-      if(subcut > 0.0) dedx -= model->ComputeDEDX(couple,particle,e,subcut); 
-      dedx *= fac;
-    }
+    if (nmod > 1) {
+      k = nmod;
+      do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k));
+      //G4cout << "k= " << k << G4endl;
+      if(k > 0) {
+        G4double elow = regModels->LowEdgeEnergy(k);
+        G4double dedx1 = ComputeDEDX(models[regModels->ModelIndex(k-1)],
+                                     couple,elow,cut,emin);
+        G4double dedx2 = ComputeDEDX(models[regModels->ModelIndex(k)],
+                                     couple,elow,cut,emin);
+        del = (dedx1 - dedx2)*elow/e;
+        //G4cout << "elow= " << elow 
+        //       << " dedx1= " << dedx1 << " dedx2= " << dedx2 << G4endl;
+      }
+    }
+    G4double dedx = ComputeDEDX(models[regModels->ModelIndex(k)],
+                                couple,e,cut,emin) + del;
 
     if(dedx < 0.0) dedx = 0.0;
     if(2 < verboseLevel) {
-        G4cout << "Material= " << couple->GetMaterial()->GetName()
-               << "   E(MeV)= " << e/MeV
-               << "  dEdx(MeV/mm)= " << dedx*mm/MeV
-               << "  dEdx0(MeV/mm)= " << dedx0*mm/MeV
-               << "  fac= " << fac
-               << G4endl;
+      G4cout << "Material= " << couple->GetMaterial()->GetName()
+             << "   E(MeV)= " << e/MeV
+             << "  dEdx(MeV/mm)= " << dedx*mm/MeV
+             << "  del= " << del*mm/MeV<< " k= " << k 
+             << " modelIdx= " << regModels->ModelIndex(k)
+             << G4endl;
     }
     aVector->PutValue(j, dedx);
@@ -660 +588 @@
                                         G4EmTableType tType)
 {
-  G4double e;
-
   size_t i = couple->GetIndex();
   G4double cut  = theCuts[i];
@@ -667 +593 @@
   if (fSubRestricted == tType) {
     tmax = cut;
-    cut  = theSubCuts[i];
-  }
-
+    if(theSubCuts.size() > 0) cut  = theSubCuts[i];
+  }
+
+  G4int reg  = 0;
+  if(nRegions > 1 && nEmModels > 1) reg  = idxOfRegionModels[i];
+  const G4RegionModels* regModels = setOfRegionModels[reg];
+  G4int nmod = regModels->NumberOfModels();
   if(1 < verboseLevel) {
-    G4cout << "G4EmModelManager::FillLambdaVector() for particle "
+    G4cout << "G4EmModelManager::FillLambdaVector() for "
            << particle->GetParticleName()
            << " in " << couple->GetMaterial()->GetName()
-           << "  Ecut(MeV)= " << cut
-           << "  Emax(MeV)= " << tmax
-           << "  Type " << tType   
+           << " Ecut(MeV)= " << cut
+           << " Emax(MeV)= " << tmax
+           << " Type " << tType   
+           << " nmod= " << nmod
            << G4endl;
   }
 
-  G4int reg  = 0;
-  if(nRegions > 1) reg  = idxOfRegionModels[i];
-  const G4RegionModels* regModels = setOfRegionModels[reg];
-  G4int nmod = regModels->NumberOfModels();
-
-  // vectors to provide continues dE/dx
-  G4DataVector factor(nmod);
-  G4DataVector eLow(nmod+1);
-  G4DataVector sigmaLow(nmod);
-  G4DataVector sigmaHigh(nmod);
-
-  if(2 < verboseLevel) {
-      G4cout << "There are " << nmod << " models for "
-             << couple->GetMaterial()->GetName() << G4endl;
-  }
-
-  // calculate factors to provide continuity of energy loss
-  factor[0] = 1.0;
-  G4int j;
-  G4int totBinsLambda = aVector->GetVectorLength();
-
-  sigmaLow[0] = 0.0;
-  eLow[0]     = 0.0;
-
-  e = regModels->LowEdgeEnergy(1);
-  eLow[1]     = e;
-  G4VEmModel* model = models[regModels->ModelIndex(0)];
-  sigmaHigh[0] = 0.0;
-  if(model) sigmaHigh[0] = model->CrossSection(couple,particle,e,cut,tmax);
-
-  if(2 < verboseLevel) {
-      G4cout << "### For material " << couple->GetMaterial()->GetName()
-             << "  " << nmod
-             << " models"
-             << " Ecut(MeV)= " << cut/MeV
-             << " Emax(MeV)= " << e/MeV
-             << " nbins= "  << totBinsLambda
-             << G4endl;
-      G4cout << " model #0   eUp= " << e 
-             << "  sigmaUp= " << sigmaHigh[0] << G4endl;
-  }
-
-
-  if(nmod > 1) {
-
-    for(j=1; j<nmod; j++) {
-
-      e  = regModels->LowEdgeEnergy(j);
-      eLow[j]   = e;
-      G4int idx = regModels->ModelIndex(j); 
-
-      sigmaLow[j] = models[idx]->CrossSection(couple,particle,e,cut,tmax);
-      e  = regModels->LowEdgeEnergy(j+1);
-      eLow[j+1] = e;
-      sigmaHigh[j] = models[idx]->CrossSection(couple,particle,e,cut,tmax);
-    }
-    if(1 < verboseLevel) {
-      G4cout << " model #0"  
-             << "  sigma(" << eLow[0] << ")=  " << sigmaLow[0]
-             << "  sigma(" << eLow[1] << ")=  " << sigmaHigh[0]
-             << G4endl;
-    }
-    for(j=1; j<nmod; j++) {
-      if(sigmaLow[j] > 0.0) {
-        factor[j] = (sigmaHigh[j-1]/sigmaLow[j] - 1.0)*eLow[j];
-      } else  factor[j] = 0.0;
-      if(1 < verboseLevel) {
-        G4cout << " model #" << j 
-               << "  sigma(" << eLow[j] << ")=  " << sigmaLow[j]
-               << "  sigma(" << eLow[j+1] << ")=  " << sigmaHigh[j]
-               << "  factor= " << factor[j]/eLow[j]
-               << G4endl;
-      }
-    }
-  }
 
   // Calculate lambda vector
-  for(j=0; j<totBinsLambda; j++) {
-
-    e = aVector->GetLowEdgeEnergy(j);
-
-    // Choose a model of energy losses
+  size_t totBinsLambda = aVector->GetVectorLength();
+  for(size_t j=0; j<totBinsLambda; ++j) {
+
+    G4double e = aVector->Energy(j);
+
+    G4double del = 0.0;
+
+    // Choose a model 
     G4int k = 0;
-    G4double fac = 1.0;
-    if (nmod > 1 && e > eLow[1]) {
-      do {
-        k++;
-        fac *= (1.0 + factor[k]/e);
-      } while ( k+1 < nmod && e > eLow[k+1] );
-    }
-
-    model = models[regModels->ModelIndex(k)];
-    G4double cross = 0.0;
-    if(model) cross = model->CrossSection(couple,particle,e,cut,tmax)*fac;
+    G4VEmModel* mod = models[regModels->ModelIndex(0)]; 
+    if (nmod > 1) {
+      k = nmod;
+      do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k));
+      if(k > 0) {
+        G4double elow = regModels->LowEdgeEnergy(k);
+        G4VEmModel* m = models[regModels->ModelIndex(k-1)]; 
+        G4double xs1  = m->CrossSection(couple,particle,elow,cut,tmax);
+        mod = models[regModels->ModelIndex(k)]; 
+        G4double xs2 = mod->CrossSection(couple,particle,elow,cut,tmax);
+        del = (xs1 - xs2)*elow/e;
+      }
+    }
+    G4double cross = mod->CrossSection(couple,particle,e,cut,tmax) + del;
+    
     if(j==0 && startFromNull) cross = 0.0;
 
@@ -779 +642 @@
       G4cout << "FillLambdaVector: " << j << ".   e(MeV)= " << e/MeV
              << "  cross(1/mm)= " << cross*mm
-             << " fac= " << fac << " k= " << k 
-             << " model= " << regModels->ModelIndex(k)
+             << " del= " << del*mm << " k= " << k 
+             << " modelIdx= " << regModels->ModelIndex(k)
              << G4endl;
     }
@@ -794 +657 @@
 {
   if(verb == 0) return;
-  for(G4int i=0; i<nRegions; i++) {
+  for(G4int i=0; i<nRegions; ++i) {
     G4RegionModels* r = setOfRegionModels[i];
     const G4Region* reg = r->Region();
-    //    if(verb > -1 || nRegions > 1) {
-    // }
     G4int n = r->NumberOfModels();  
-    if(verb > -1 || n > 0) {
+    if(n > 0) {
       G4cout << "      ===== EM models for the G4Region  " << reg->GetName()
              << " ======" << G4endl;;
-      for(G4int j=0; j<n; j++) {
+      for(G4int j=0; j<n; ++j) {
         const G4VEmModel* m = models[r->ModelIndex(j)];
         G4cout << std::setw(20);
@@ -813 +674 @@
       }  
     }
-  }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+    if(1 == nEmModels) break; 
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

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

@@ -25 +25 @@
 //
 // $Id: G4EmMultiModel.cc,v 1.6 2007/05/22 17:31:58 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4EmProcessOptions.cc,v 1.26 2009/02/18 14:43:27 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4EmProcessOptions.cc,v 1.27 2009/10/29 19:25:28 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -516 +516 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+void G4EmProcessOptions::SetFactorForAngleLimit(G4double val)
+{
+  theManager->SetFactorForAngleLimit(val);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
 void G4EmProcessOptions::SetLPMFlag(G4bool val)
 {

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4EmSaturation.cc,v 1.9 2008/11/12 15:37:33 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4EmSaturation.cc,v 1.10 2009/09/25 09:16:40 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -46 +46 @@
 
 #include "G4EmSaturation.hh"
-#include "G4Gamma.hh"
-#include "G4Electron.hh"
-#include "G4Neutron.hh"
-#include "G4Proton.hh"
 #include "G4LossTableManager.hh"
 #include "G4NistManager.hh"
 #include "G4Material.hh"
 #include "G4MaterialCutsCouple.hh"
+#include "G4Electron.hh"
+#include "G4Proton.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -66 +64 @@
   curChargeSq = 1.0;
   nMaterials = 0;
-  Initialise();
+  electron = 0;
+  Initialise(); 
 }
 
@@ -90 +89 @@
   if(bfactor > 0.0) { 
 
-    // atomic relaxations
-    if(p == gamma) {
+    G4int pdgCode = p->GetPDGEncoding();
+    // atomic relaxations for gamma incident
+    if(22 == pdgCode) {
       evis /= (1.0 + bfactor*edep/manager->GetRange(electron,edep,couple));
 
@@ -101 +101 @@
       if(nloss < 0.0) nloss = 0.0;
       G4double eloss = edep - nloss;
-      if(p == neutron || eloss < 0.0 || length <= 0.0) {
+
+      // neutrons
+      if(2112 == pdgCode || eloss < 0.0 || length <= 0.0) {
         nloss = edep;
         eloss = 0.0;
@@ -111 +113 @@
       // non-ionizing energy loss
       if(nloss > 0.0) {
+        if(!proton) {proton = G4Proton::Proton();}
         G4double escaled = nloss*curRatio;
         G4double s = manager->GetRange(proton,escaled,couple)/curChargeSq; 
@@ -146 +149 @@
 G4double G4EmSaturation::FindBirksCoefficient(const G4Material* mat)
 {
+  // electron should exist in any case
+  if(!manager) {
+    manager = G4LossTableManager::Instance();
+    nist    = G4NistManager::Instance();
+    electron= G4Electron::Electron();
+    proton  = 0;
+  }
+
   if(mat == curMaterial) return curBirks;
 
@@ -161 +172 @@
       return curBirks;
     }
-  }
-
-  if(!manager) {
-    manager = G4LossTableManager::Instance();
-    nist    = G4NistManager::Instance();
-    gamma   = G4Gamma::Gamma();
-    electron= G4Electron::Electron();
-    proton  = G4Proton::Proton();
-    neutron = G4Neutron::Neutron();
   }
 

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

@@ -25 +25 @@
 //
 //
-// $Id: G4EnergyLossMessenger.cc,v 1.37 2009/02/18 14:43:27 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4EnergyLossMessenger.cc,v 1.38 2009/10/29 19:25:28 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -196 +196 @@
   dedxCmd->SetGuidance("Set number of bins for DEDX tables");
   dedxCmd->SetParameterName("binsDEDX",true);
-  dedxCmd->SetDefaultValue(120);
+  dedxCmd->SetDefaultValue(77);
   dedxCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
 
@@ -202 +202 @@
   lamCmd->SetGuidance("Set number of bins for Lambda tables");
   lamCmd->SetParameterName("binsL",true);
-  lamCmd->SetDefaultValue(120);
+  lamCmd->SetDefaultValue(77);
   lamCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
 
@@ -242 +242 @@
   frCmd->SetParameterName("Fr",true);
   frCmd->SetRange("Fr>0");
-  frCmd->SetDefaultValue(0.02);
+  frCmd->SetDefaultValue(0.04);
   frCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
 
@@ -252 +252 @@
   fgCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
 
+  mscfCmd = new G4UIcmdWithADouble("/process/msc/FactorForAngleLimit",this);
+  mscfCmd->SetGuidance("Set factor for computation of a limit for -t (invariant trasfer)");
+  mscfCmd->SetParameterName("Fact",true);
+  mscfCmd->SetRange("Fact>0");
+  mscfCmd->SetDefaultValue(1.);
+  mscfCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
+
   skinCmd = new G4UIcmdWithADouble("/process/msc/Skin",this);
   skinCmd->SetGuidance("Set skin parameter for msc processes");
@@ -258 +265 @@
 
   angCmd = new G4UIcmdWithADoubleAndUnit("/process/msc/ThetaLimit",this);
-  angCmd->SetGuidance("Set the limit on the polar angle");
+  angCmd->SetGuidance("Set the limit on the polar angle for msc and single scattering");
   angCmd->SetParameterName("theta",true);
   angCmd->SetUnitCategory("Angle");
@@ -297 +304 @@
   delete skinCmd;
   delete angCmd;
+  delete mscfCmd;
 }
 
@@ -432 +440 @@
     G4UImanager::GetUIpointer()->ApplyCommand("/run/physicsModified");
   }
+  if (command == mscfCmd) {
+    opt->SetFactorForAngleLimit(mscfCmd->GetNewDoubleValue(newValue));
+    G4UImanager::GetUIpointer()->ApplyCommand("/run/physicsModified");
+  }
   if (command == angCmd) { 
     opt->SetPolarAngleLimit(angCmd->GetNewDoubleValue(newValue));

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

@@ -26 +26 @@
 //
 // $Id: G4EnergyLossTables.cc,v 1.34 2008/07/08 10:57:22 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4LossTableBuilder.cc,v 1.28 2009/02/18 16:24:47 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4LossTableBuilder.cc,v 1.32 2009/08/11 17:24:53 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -47 +47 @@
 // 16-01-07 Fill new (not old) DEDX table (V.Ivanchenko)
 // 12-02-07 Use G4LPhysicsFreeVector for the inverse range table (V.Ivanchenko)
+// 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
 //
 // Class Description:
@@ -85 +86 @@
   if(0 >= n_vectors) return;
 
-  G4bool b;
-
-  G4PhysicsVector* pv = (*(list[0]))[0];
-  size_t nbins = pv->GetVectorLength();
-  G4double elow = pv->GetLowEdgeEnergy(0);
-  G4double ehigh = pv->GetLowEdgeEnergy(nbins);
+  G4PhysicsLogVector* pv0 = static_cast<G4PhysicsLogVector*>((*(list[0]))[0]);
+  size_t npoints = pv0->GetVectorLength();
   for (size_t i=0; i<n_vectors; i++) {
 
-    pv = new G4PhysicsLogVector(elow, ehigh, nbins);
+    G4PhysicsLogVector* pv = new G4PhysicsLogVector(*pv0);
+    //    pv = new G4PhysicsLogVector(elow, ehigh, npoints-1);
     pv->SetSpline(splineFlag);
-    for (size_t j=0; j<nbins; j++) {
+    for (size_t j=0; j<npoints; j++) {
       G4double dedx = 0.0;
-      G4double energy = pv->GetLowEdgeEnergy(j);
-
       for (size_t k=0; k<n_processes; k++) {
-        dedx += ((*(list[k]))[i])->GetValue(energy, b);
+        G4PhysicsVector* pv1   = (*(list[k]))[i];
+        dedx += (*pv1)[j];
       }
       pv->PutValue(j, dedx);
-      G4PhysicsTableHelper::SetPhysicsVector(dedxTable, i, pv);
     }
+    if(splineFlag) pv->FillSecondDerivatives();
+    G4PhysicsTableHelper::SetPhysicsVector(dedxTable, i, pv);
   }
 }
@@ -118 +116 @@
   if(!n_vectors) return;
 
-  G4bool b;
   size_t n = 100;
   G4double del = 1.0/(G4double)n;
@@ -125 +122 @@
 
     if (rangeTable->GetFlag(i) || !isIonisation) {
-      G4PhysicsVector* pv = (*dedxTable)[i];
-      size_t nbins = pv->GetVectorLength();
-      size_t bin0  = 0;
-      G4double elow = pv->GetLowEdgeEnergy(0);
-      G4double ehigh = pv->GetLowEdgeEnergy(nbins);
-      G4double dedx1  = pv->GetValue(elow, b);
+      G4PhysicsLogVector* pv = 
+        static_cast<G4PhysicsLogVector*>((*dedxTable)[i]);
+      size_t npoints = pv->GetVectorLength();
+      size_t bin0    = 0;
+      G4double elow  = pv->Energy(0);
+      G4double ehigh = pv->Energy(npoints-1);
+      G4double dedx1 = pv->Value(elow);
 
       //G4cout << "nbins= " << nbins << " dedx1= " << dedx1 << G4endl;
@@ -136 +134 @@
       // protection for specific cases dedx=0
       if(dedx1 == 0.0) {
-        for (size_t k=1; k<nbins-1; k++) {
+        for (size_t k=1; k<npoints; k++) {
           bin0++;
-          elow  = pv->GetLowEdgeEnergy(k);
-          dedx1 = pv->GetValue(elow, b);
+          elow  = pv->Energy(k);
+          dedx1 = (*pv)[k];
           if(dedx1 > 0.0) break;
         }
-        nbins -= bin0;
-      }
-
-      //G4cout << "nbins= " << nbins << " elow= " << elow << " ehigh= " << ehigh << G4endl;
+        npoints -= bin0;
+      }
+      // G4cout<<"New Range vector" << G4endl;
+      // G4cout<<"nbins= "<<npoints-1<<" elow= "<<elow<<" ehigh= "<<ehigh<<G4endl;
       // initialisation of a new vector
-      G4PhysicsLogVector* v = new G4PhysicsLogVector(elow, ehigh, nbins);
+      if(npoints < 2) npoints = 2;
+      G4PhysicsLogVector* v;
+      if(0 == bin0) { v = new G4PhysicsLogVector(*pv); }
+      else          { v = new G4PhysicsLogVector(elow, ehigh, npoints-1); }
       // dedx is exect zero
-      if(nbins <= 2) {
+      if(2 == npoints) {
         v->PutValue(0,1000.);
         v->PutValue(1,2000.);
@@ -162 +163 @@
       G4double energy1 = elow;
 
-      for (size_t j=1; j<nbins; j++) {
-
-        G4double energy2 = pv->GetLowEdgeEnergy(j+bin0);
+      for (size_t j=1; j<npoints; j++) {
+
+        G4double energy2 = pv->Energy(j+bin0);
         G4double de      = (energy2 - energy1) * del;
         G4double energy  = energy2 + de*0.5;
+        G4double sum = 0.0;
         for (size_t k=0; k<n; k++) {
           energy -= de;
-          dedx1 = pv->GetValue(energy, b);
-          if(dedx1 > 0.0) range += de/dedx1;
+          dedx1 = pv->Value(energy);
+          if(dedx1 > 0.0) sum += de/dedx1;
         }
-
+        range += sum;
         //      G4cout << "Range i= " <<i << " j= " << j << G4endl;
         v->PutValue(j,range);
         energy1 = energy2;
       }
+      if(splineFlag) v->FillSecondDerivatives();
       G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
     }
@@ -191 +194 @@
   size_t n_vectors = rangeTable->length();
   if(!n_vectors) return;
-  G4bool b;
 
   for (size_t i=0; i<n_vectors; i++) {
@@ -197 +199 @@
     if (invRangeTable->GetFlag(i) || !isIonisation) {
       G4PhysicsVector* pv = (*rangeTable)[i];
-      size_t nbins   = pv->GetVectorLength();
-      G4double elow  = pv->GetLowEdgeEnergy(0);
-      G4double ehigh = pv->GetLowEdgeEnergy(nbins-1);
-      G4double rlow  = pv->GetValue(elow, b);
-      G4double rhigh = pv->GetValue(ehigh, b);
+      size_t npoints = pv->GetVectorLength();
+      G4double rlow  = (*pv)[0];
+      G4double rhigh = (*pv)[npoints-1];
       
-      G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(nbins,rlow,rhigh);
+      G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(npoints,rlow,rhigh);
       v->SetSpline(splineFlag);
 
-      for (size_t j=0; j<nbins; j++) {
-        G4double e  = pv->GetLowEdgeEnergy(j);
-        G4double r  = pv->GetValue(e, b);
+      for (size_t j=0; j<npoints; j++) {
+        G4double e  = pv->Energy(j);
+        G4double r  = (*pv)[j];
         v->PutValues(j,r,e);
       }
-      v->PutValues(nbins,rhigh+rlow,ehigh);
+      if(splineFlag) v->FillSecondDerivatives();
 
       G4PhysicsTableHelper::SetPhysicsVector(invRangeTable, i, v);

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4LossTableManager.cc,v 1.96 2009/04/09 16:10:57 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4LossTableManager.cc,v 1.97 2009/10/29 19:25:28 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -176 +176 @@
   splineFlag = true;
   bremsTh = DBL_MAX;
+  factorForAngleLimit = 1.0;
   verbose = 1;
   tableBuilder->SetSplineFlag(splineFlag);
@@ -920 +921 @@
 }
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4LossTableManager::SetFactorForAngleLimit(G4double val) 
+{
+  if(val > 0.0) factorForAngleLimit = val;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4LossTableManager::FactorForAngleLimit() const
+{
+  return factorForAngleLimit;
+}
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 

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

@@ -25 +25 @@
 //
 // $Id: G4VEmFluctuationModel.cc,v 1.4 2009/02/19 11:25:50 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4VEmModel.cc,v 1.27 2009/05/26 15:00:49 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4VEmModel.cc,v 1.30 2009/09/23 14:42:47 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -62 +62 @@
 G4VEmModel::G4VEmModel(const G4String& nam):
   fluc(0), name(nam), lowLimit(0.1*keV), highLimit(100.0*TeV), 
+  eMinActive(0.0),eMaxActive(DBL_MAX),
   polarAngleLimit(0.0),secondaryThreshold(DBL_MAX),theLPMflag(false),
   pParticleChange(0),nuclearStopping(false),
@@ -94 +95 @@
   } else {
     p = new G4ParticleChangeForLoss();
+    pParticleChange = p;
   }
   return p;
@@ -107 +109 @@
   } else {
     p = new G4ParticleChangeForGamma();
+    pParticleChange = p;
   }
   return p;
@@ -132 +135 @@
   // initialise vector
   for(G4int i=0; i<numOfCouples; i++) {
-    const G4MaterialCutsCouple* couple =
-      theCoupleTable->GetMaterialCutsCouple(i);
-    const G4Material* material = couple->GetMaterial();
-    G4int idx = couple->GetIndex();
+    currentCouple = theCoupleTable->GetMaterialCutsCouple(i);
+    const G4Material* material = currentCouple->GetMaterial();
+    G4int idx = currentCouple->GetIndex();
 
     // selector already exist check if should be deleted
@@ -202 +204 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
+void G4VEmModel::DefineForRegion(const G4Region*) 
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
 G4double G4VEmModel::MinEnergyCut(const G4ParticleDefinition*,
                                   const G4MaterialCutsCouple*)
@@ -249 +256 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-void G4VEmModel::DefineForRegion(const G4Region*) 
-{}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
 void G4VEmModel::SetupForMaterial(const G4ParticleDefinition*,
                                   const G4Material*, G4double)

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4VEmProcess.cc,v 1.66 2009/04/17 10:35:32 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4VEmProcess.cc,v 1.79 2009/11/10 20:30:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -52 +52 @@
 // 12-04-07 remove double call to Clear model manager (V.Ivanchenko)
 // 27-10-07 Virtual functions moved to source (V.Ivanchenko)
+// 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
 //
 // Class Description:
@@ -172 +173 @@
 {
   G4int n = emModels.size();
-  if(index >= n) for(G4int i=n; i<index+1; i++) {emModels.push_back(0);}
+  if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} }
   emModels[index] = p;
 }
@@ -181 +182 @@
 {
   G4VEmModel* p = 0;
-  if(index >= 0 && index <  G4int(emModels.size())) p = emModels[index];
+  if(index >= 0 && index <  G4int(emModels.size())) { p = emModels[index]; }
   return p;
 }
@@ -218 +219 @@
     Clear();
     InitialiseProcess(particle);
+
+    // initialisation of models
+    G4int nmod = modelManager->NumberOfModels();
+    for(G4int i=0; i<nmod; ++i) {
+      G4VEmModel* mod = modelManager->GetModel(i);
+      mod->SetPolarAngleLimit(polarAngleLimit);
+      if(mod->HighEnergyLimit() > maxKinEnergy) {
+        mod->SetHighEnergyLimit(maxKinEnergy);
+      }
+    }
+
     theCuts = modelManager->Initialise(particle,secondaryParticle,2.,verboseLevel);
     const G4ProductionCutsTable* theCoupleTable=
@@ -224 +236 @@
     theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
     theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
+
+    // prepare tables
     if(buildLambdaTable){
       theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
@@ -237 +251 @@
     idxDERegions = new G4bool[numOfCouples];
 
-    for (size_t j=0; j<numOfCouples; j++) {
+    for (size_t j=0; j<numOfCouples; ++j) {
 
       const G4MaterialCutsCouple* couple =
@@ -243 +257 @@
       const G4ProductionCuts* pcuts = couple->GetProductionCuts();
       G4bool reg = false;
-      for(G4int i=0; i<nDERegions; i++) {
+      for(G4int i=0; i<nDERegions; ++i) {
         if(deRegions[i]) {
           if(pcuts == deRegions[i]->GetProductionCuts()) reg = true;
@@ -253 +267 @@
   if (1 < verboseLevel && nDERegions>0) {
     G4cout << " Deexcitation is activated for regions: " << G4endl;
-    for (G4int i=0; i<nDERegions; i++) {
+    for (G4int i=0; i<nDERegions; ++i) {
       const G4Region* r = deRegions[i];
       G4cout << "           " << r->GetName() << G4endl;
@@ -264 +278 @@
 void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
 {
+  G4String partname = part.GetParticleName();
   if(1 < verboseLevel) {
     G4cout << "G4VEmProcess::BuildPhysicsTable() for "
            << GetProcessName()
-           << " and particle " << part.GetParticleName()
+           << " and particle " << partname
            << " buildLambdaTable= " << buildLambdaTable
            << G4endl;
@@ -276 +291 @@
     FindLambdaMax();
   }
-  if(0 < verboseLevel) PrintInfoDefinition();
+
+  // reduce printout for nuclear stopping
+  G4bool gproc = true;
+  if(GetProcessName() == "nuclearStopping" && 
+     partname != "GenericIon" && partname != "alpha") { gproc = false; } 
+
+  if(gproc && 0 < verboseLevel) { PrintInfoDefinition(); }
 
   if(1 < verboseLevel) {
     G4cout << "G4VEmProcess::BuildPhysicsTable() done for "
            << GetProcessName()
-           << " and particle " << part.GetParticleName()
+           << " and particle " << partname
            << G4endl;
   }
@@ -301 +322 @@
         G4ProductionCutsTable::GetProductionCutsTable();
   size_t numOfCouples = theCoupleTable->GetTableSize();
-  for(size_t i=0; i<numOfCouples; i++) {
+
+  G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
+
+  G4PhysicsLogVector* aVector = 0;
+  G4PhysicsLogVector* bVector = 0;
+
+  for(size_t i=0; i<numOfCouples; ++i) {
 
     if (theLambdaTable->GetFlag(i)) {
 
       // create physics vector and fill it
-      const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
-      G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
+      const G4MaterialCutsCouple* couple = 
+        theCoupleTable->GetMaterialCutsCouple(i);
+      if(!bVector) {
+        aVector = 
+          static_cast<G4PhysicsLogVector*>(LambdaPhysicsVector(couple));
+        bVector = aVector;
+      } else {
+        aVector = new G4PhysicsLogVector(*bVector);
+      }
+        //      G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
+      aVector->SetSpline(splineFlag);
       modelManager->FillLambdaVector(aVector, couple, startFromNull);
+      if(splineFlag) aVector->FillSecondDerivatives();
       G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
     }
@@ -364 +401 @@
   if(previousStepSize <= DBL_MIN) theNumberOfInteractionLengthLeft = -1.0;
   InitialiseStep(track);
+  if(!currentModel->IsActive(preStepKinEnergy)) return x;
 
   if(preStepKinEnergy < mfpKinEnergy) {
@@ -444 +482 @@
   }
 
-  SelectModel(finalT);
+  SelectModel(finalT, currentCoupleIndex);
+  if(!currentModel->IsActive(finalT)) return &fParticleChange;
   if(useDeexcitation) {
-    currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);
+    currentModel->SetDeexcitationFlag(idxDERegions[currentCoupleIndex]);
   }
   /*  
@@ -464 +503 @@
                                   currentCouple, 
                                   track.GetDynamicParticle(),
-                                  (*theCuts)[currentMaterialIndex]);
+                                  (*theCuts)[currentCoupleIndex]);
 
   // save secondaries
@@ -473 +512 @@
     G4double edep = fParticleChange.GetLocalEnergyDeposit();
      
-    for (G4int i=0; i<num; i++) {
+    for (G4int i=0; i<num; ++i) {
       G4DynamicParticle* dp = secParticles[i];
       const G4ParticleDefinition* p = dp->GetDefinition();
@@ -480 +519 @@
       if(applyCuts) {
         if (p == theGamma) {
-          if (e < (*theCutsGamma)[currentMaterialIndex]) good = false;
+          if (e < (*theCutsGamma)[currentCoupleIndex]) good = false;
 
         } else if (p == theElectron) {
-          if (e < (*theCutsElectron)[currentMaterialIndex]) good = false;
+          if (e < (*theCutsElectron)[currentCoupleIndex]) good = false;
 
         } else if (p == thePositron) {
-          if (electron_mass_c2 < (*theCutsGamma)[currentMaterialIndex] &&
-              e < (*theCutsPositron)[currentMaterialIndex]) {
+          if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] &&
+              e < (*theCutsPositron)[currentCoupleIndex]) {
             good = false;
             e += 2.0*electron_mass_c2;
@@ -539 +578 @@
 G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
                                           const G4String& directory,
-                                                G4bool ascii)
+                                          G4bool ascii)
 {
   if(1 < verboseLevel) {
@@ -558 +597 @@
   if ( yes ) {
     if (0 < verboseLevel) {
-      G4cout << "Lambda table for " << particleName << " is Retrieved from <"
+      G4cout << "Lambda table for " << particleName 
+             << " is Retrieved from <"
              << filename << ">"
              << G4endl;
@@ -564 +604 @@
     if((G4LossTableManager::Instance())->SplineFlag()) {
       size_t n = theLambdaTable->length();
-      for(size_t i=0; i<n; i++) {(* theLambdaTable)[i]->SetSpline(true);}
+      for(size_t i=0; i<n; ++i) {
+        if((* theLambdaTable)[i]) {
+          (* theLambdaTable)[i]->SetSpline(true);
+        }
+      }
     }
   } else {
@@ -587 +631 @@
   // the region is in the list
   if (nDERegions) {
-    for (G4int i=0; i<nDERegions; i++) {
+    for (G4int i=0; i<nDERegions; ++i) {
       if (reg == deRegions[i]) {
         if(!val) deRegions[i] = 0;
@@ -613 +657 @@
   DefineMaterial(couple);
   G4double cross = 0.0;
-  G4bool b;
   if(theLambdaTable) {
-    cross = (((*theLambdaTable)[currentMaterialIndex])->
-                           GetValue(kineticEnergy, b));
+    cross = (((*theLambdaTable)[currentCoupleIndex])->Value(kineticEnergy));
   } else {
-    SelectModel(kineticEnergy);
+    SelectModel(kineticEnergy, currentCoupleIndex);
     cross = currentModel->CrossSectionPerVolume(currentMaterial,
                                                 particle,kineticEnergy);
@@ -650 +692 @@
   theEnergyOfCrossSectionMax = new G4double [n];
   theCrossSectionMax = new G4double [n];
-  G4bool b;
-
-  for (size_t i=0; i<n; i++) {
+
+  for (size_t i=0; i<n; ++i) {
     pv = (*theLambdaTable)[i];
     emax = DBL_MAX;
@@ -658 +699 @@
     if(pv) {
       size_t nb = pv->GetVectorLength();
-      emax = pv->GetLowEdgeEnergy(nb);
+      emax = DBL_MAX;
       smax = 0.0;
-      for (size_t j=0; j<nb; j++) {
-        e = pv->GetLowEdgeEnergy(j);
-        s = pv->GetValue(e,b);
-        if(s > smax) {
-          smax = s;
-          emax = e;
+      if(nb > 0) {
+        for (size_t j=0; j<nb; ++j) {
+          e = pv->Energy(j);
+          s = (*pv)(j);
+          if(s > smax) {
+            smax = s;
+            emax = e;
+          }
         }
       }

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

@@ -26 +26 @@
 //
 // $Id: G4VEnergyLoss.cc,v 1.46 2006/06/29 19:55:21 gunter Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4VEnergyLossProcess.cc,v 1.149 2009/04/17 10:35:32 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4VEnergyLossProcess.cc,v 1.158 2009/10/29 18:07:08 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -107 +107 @@
 // 25-04-07 move initialisation of safety helper to BuildPhysicsTable (VI)
 // 27-10-07 Virtual functions moved to source (V.Ivanchenko)
+// 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
 //
 // Class Description:
@@ -230 +231 @@
   vstrag->SetSpline(true);
   G4double s[nrbins] = {-0.2, -0.85, -1.3, -1.578, -1.76, -1.85, -1.9};
-  for(G4int i=0; i<nrbins; i++) {vstrag->PutValue(i, s[i]);}
+  for(G4int i=0; i<nrbins; ++i) {vstrag->PutValue(i, s[i]);}
 }
 
@@ -351 +352 @@
 {
   G4int n = emModels.size();
-  if(index >= n) for(G4int i=n; i<index+1; i++) {emModels.push_back(0);}
+  if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} }
   emModels[index] = p;
 }
@@ -360 +361 @@
 {
   G4VEmModel* p = 0;
-  if(index >= 0 && index <  G4int(emModels.size())) p = emModels[index];
+  if(index >= 0 && index <  G4int(emModels.size())) { p = emModels[index]; }
   return p;
 }
@@ -466 +467 @@
     G4double q = initialCharge/baseParticle->GetPDGCharge();
     chargeSqRatio = q*q;
-    if(chargeSqRatio > 0.0) reduceFactor = 1.0/(chargeSqRatio*massRatio);
+    if(chargeSqRatio > 0.0) { reduceFactor = 1.0/(chargeSqRatio*massRatio); }
+  }
+
+  // initialisation of models
+  G4int nmod = modelManager->NumberOfModels();
+  for(G4int i=0; i<nmod; ++i) {
+    G4VEmModel* mod = modelManager->GetModel(i);
+    if(mod->HighEnergyLimit() > maxKinEnergy) {
+      mod->SetHighEnergyLimit(maxKinEnergy);
+    }
   }
 
@@ -483 +493 @@
     if(nDERegions>0) idxDERegions = new G4bool[numOfCouples];
   
-    for (size_t j=0; j<numOfCouples; j++) {
+    for (size_t j=0; j<numOfCouples; ++j) {
 
       const G4MaterialCutsCouple* couple = 
@@ -491 +501 @@
       if(nSCoffRegions>0) {
         G4bool reg = false;
-        for(G4int i=0; i<nSCoffRegions; i++) {
+        for(G4int i=0; i<nSCoffRegions; ++i) {
           if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = true;
         }
@@ -498 +508 @@
       if(nDERegions>0) {
         G4bool reg = false;
-        for(G4int i=0; i<nDERegions; i++) {
+        for(G4int i=0; i<nDERegions; ++i) {
           if( pcuts == deRegions[i]->GetProductionCuts()) reg = true;
         }
@@ -517 +527 @@
     if (nSCoffRegions) {
       G4cout << " SubCutoff Regime is ON for regions: " << G4endl;
-      for (G4int i=0; i<nSCoffRegions; i++) {
+      for (G4int i=0; i<nSCoffRegions; ++i) {
         const G4Region* r = scoffRegions[i];
         G4cout << "           " << r->GetName() << G4endl;
@@ -524 +534 @@
     if (nDERegions) {
       G4cout << " Deexcitation is ON for regions: " << G4endl;
-      for (G4int i=0; i<nDERegions; i++) {
+      for (G4int i=0; i<nDERegions; ++i) {
         const G4Region* r = deRegions[i];
         G4cout << "           " << r->GetName() << G4endl;
@@ -556 +566 @@
     }
   }
+
+  // Added tracking cut to avoid tracking artifacts
+  if(isIonisation) fParticleChange.SetLowEnergyLimit(lowestKinEnergy);
 
   if(1 < verboseLevel) {
@@ -613 +626 @@
   if(!table) return table;
 
-  for(size_t i=0; i<numOfCouples; i++) {
+  G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
+  G4PhysicsLogVector* aVector = 0;
+  G4PhysicsLogVector* bVector = 0;
+
+  for(size_t i=0; i<numOfCouples; ++i) {
 
     if(1 < verboseLevel) {
@@ -624 +641 @@
       const G4MaterialCutsCouple* couple = 
         theCoupleTable->GetMaterialCutsCouple(i);
-      G4PhysicsVector* aVector = new G4PhysicsLogVector(emin, emax, bin);
-      aVector->SetSpline((G4LossTableManager::Instance())->SplineFlag());
+      if(!bVector) {
+        aVector = new G4PhysicsLogVector(emin, emax, bin);
+        bVector = aVector;
+      } else {
+        aVector = new G4PhysicsLogVector(*bVector);
+      }
+      //      G4PhysicsVector* aVector = new G4PhysicsLogVector(emin, emax, bin);
+      aVector->SetSpline(splineFlag);
 
       modelManager->FillDEDXVector(aVector, couple, tType);
+      if(splineFlag) aVector->FillSecondDerivatives();
 
       // Insert vector for this material into the table
@@ -676 +700 @@
   size_t numOfCouples = theCoupleTable->GetTableSize();
 
-  for(size_t i=0; i<numOfCouples; i++) {
+  G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
+
+  for(size_t i=0; i<numOfCouples; ++i) {
 
     if (table->GetFlag(i)) {
@@ -686 +712 @@
       if(fSubRestricted == tType) cut = (*theSubCuts)[i]; 
       G4PhysicsVector* aVector = LambdaPhysicsVector(couple, cut);
+      aVector->SetSpline(splineFlag);
+
       modelManager->FillLambdaVector(aVector, couple, true, tType);
+      if(splineFlag) aVector->FillSecondDerivatives();
 
       // Insert vector for this material into the table
@@ -787 +816 @@
   // the region is in the list
   if (nSCoffRegions) {
-    for (G4int i=0; i<nSCoffRegions; i++) {
+    for (G4int i=0; i<nSCoffRegions; ++i) {
       if (reg == scoffRegions[i]) {
         if(!val) deRegions[i] = 0;
@@ -799 +828 @@
     useSubCutoff = true;
     scoffRegions.push_back(reg);
-    nSCoffRegions++;
+    ++nSCoffRegions;
   } else {
     useSubCutoff = false;
@@ -815 +844 @@
   // the region is in the list
   if (nDERegions) {
-    for (G4int i=0; i<nDERegions; i++) {
+    for (G4int i=0; i<nDERegions; ++i) {
       if (reg == deRegions[i]) {
         if(!val) deRegions[i] = 0;
@@ -827 +856 @@
     useDeexcitation = true;
     deRegions.push_back(reg);
-    nDERegions++;
+    ++nDERegions;
   } else {
     useDeexcitation = false;
@@ -883 +912 @@
   preStepScaledEnergy = preStepKinEnergy*massRatio;
   SelectModel(preStepScaledEnergy);
+  if(!currentModel->IsActive(preStepScaledEnergy)) return x;
 
   if(isIon) {
@@ -951 +981 @@
   fParticleChange.InitializeForAlongStep(track);
   // The process has range table - calculate energy loss
-  if(!isIonisation) return &fParticleChange;
+  if(!isIonisation || !currentModel->IsActive(preStepScaledEnergy)) {
+    return &fParticleChange;
+  }
 
   // Get the actual (true) Step length
@@ -1076 +1108 @@
         /*
         if(nProcesses > 0) {
-          for(G4int i=0; i<nProcesses; i++) {
+          for(G4int i=0; i<nProcesses; ++i) {
             (scProcesses[i])->SampleSubCutSecondaries(
                 scTracks, step, (scProcesses[i])->
@@ -1088 +1120 @@
           G4ThreeVector mom = dynParticle->GetMomentum();
           fParticleChange.SetNumberOfSecondaries(n);
-          for(G4int i=0; i<n; i++) {
+          for(G4int i=0; i<n; ++i) {
             G4Track* t = scTracks[i];
             G4double e = t->GetKineticEnergy();
@@ -1182 +1214 @@
   const G4DynamicParticle* dp = track->GetDynamicParticle();
   G4double e = dp->GetKineticEnergy()*massRatio;
-  G4bool b;
-  G4double cross = chargeSqRatio*(((*theSubLambdaTable)[idx])->GetValue(e,b));
+  G4double cross = chargeSqRatio*(((*theSubLambdaTable)[idx])->Value(e));
   G4double length = step.GetStepLength();
 
@@ -1220 +1251 @@
     G4ThreeVector r = prepoint + fragment*dr;
     std::vector<G4DynamicParticle*>::iterator it;
-    for(it=secParticles.begin(); it!=secParticles.end(); it++) {
+    for(it=secParticles.begin(); it!=secParticles.end(); ++it) {
 
       G4bool addSec = true;
@@ -1226 +1257 @@
       // do not track very low-energy delta-electrons
       if(theSecondaryRangeTable && (*it)->GetDefinition() == theElectron) {
-        G4bool b;
         G4double ekin = (*it)->GetKineticEnergy();
-        G4double rg = ((*theSecondaryRangeTable)[idx]->GetValue(ekin, b));
+        G4double rg = ((*theSecondaryRangeTable)[idx]->Value(ekin));
         //          if(rg < currentMinSafety) {
         if(rg < safetyHelper->ComputeSafety(r)) {
@@ -1265 +1295 @@
 
   G4double postStepScaledEnergy = finalT*massRatio;
+
+  if(!currentModel->IsActive(postStepScaledEnergy)) return &fParticleChange;
   /*
   if(-1 < verboseLevel) {
@@ -1283 +1315 @@
              << " < " << lx << " (postLambda) "
              << G4endl;
-      nWarnings++;
+      ++nWarnings;
     }
     */
@@ -1308 +1340 @@
   if(num > 0) {
     fParticleChange.SetNumberOfSecondaries(num);
-    for (G4int i=0; i<num; i++) {
+    for (G4int i=0; i<num; ++i) {
       fParticleChange.AddSecondary(secParticles[i]);
     }
@@ -1390 +1422 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
 
-G4bool G4VEnergyLossProcess::RetrievePhysicsTable(
-       const G4ParticleDefinition* part, const G4String& directory,
-       G4bool ascii)
+G4bool 
+G4VEnergyLossProcess::RetrievePhysicsTable(const G4ParticleDefinition* part, 
+                                           const G4String& directory,
+                                           G4bool ascii)
 {
   G4bool res = true;
@@ -1405 +1438 @@
   if(particle == part) {
 
-    //    G4bool yes = true;
     if ( !baseParticle ) {
 
@@ -1412 +1444 @@
         {fpi = false;}
 
-      if(!RetrieveTable(part,theIonisationTable,ascii,directory,"Ionisation",false)) 
+      // ionisation table keeps individual dEdx and not sum of sub-processes
+      if(!RetrieveTable(part,theDEDXTable,ascii,directory,"Ionisation",false)) 
         {fpi = false;}
 
@@ -1440 +1473 @@
 
       if(!fpi) yes = false;
-      if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes)) 
+      if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes))
         {res = false;}
     }
@@ -1465 +1498 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
 
-G4bool G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part, 
-                                           G4PhysicsTable* aTable, G4bool ascii,
-                                           const G4String& directory,
-                                           const G4String& tname,
-                                           G4bool mandatory)
+G4bool 
+G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part, 
+                                    G4PhysicsTable* aTable, 
+                                    G4bool ascii,
+                                    const G4String& directory,
+                                    const G4String& tname,
+                                    G4bool mandatory)
 {
   G4bool res = true;
@@ -1475 +1510 @@
   G4bool yes = aTable->ExistPhysicsTable(filename);
   if(yes) {
+    if(!aTable) aTable = G4PhysicsTableHelper::PreparePhysicsTable(0);
     yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii);
+
     if((G4LossTableManager::Instance())->SplineFlag()) {
       size_t n = aTable->length();
-      for(size_t i=0; i<n; i++) {(*aTable)[i]->SetSpline(true);}
+      for(size_t i=0; i<n; ++i) {
+        if((*aTable)[i]) {
+          (*aTable)[i]->SetSpline(true);
+        }
+      }
     }
   }
@@ -1525 +1566 @@
   DefineMaterial(couple);
   G4double cross = 0.0;
-  G4bool b;
   if(theLambdaTable) {
-    cross = 
-      ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b);
+    cross = ((*theLambdaTable)[currentMaterialIndex])->Value(kineticEnergy);
   } else {
     SelectModel(kineticEnergy);
-    cross = 
-      currentModel->CrossSectionPerVolume(currentMaterial,
-                                          particle, kineticEnergy,
-                                          (*theCuts)[currentMaterialIndex]);
+    cross = currentModel->CrossSectionPerVolume(currentMaterial,
+                                                particle, kineticEnergy,
+                                                (*theCuts)[currentMaterialIndex]);
   }
 
@@ -1587 +1625 @@
                  const G4MaterialCutsCouple* couple, G4double cut)
 {
-  //  G4double cut  = (*theCuts)[couple->GetIndex()];
-  //  G4int nbins = nLambdaBins;
   G4double tmin = 
     std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
@@ -1595 +1631 @@
   G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
   v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
-
   return v;
 }
@@ -1607 +1642 @@
   if(p->GetProcessName() != "eBrem") add = false;
   if(add && nProcesses > 0) {
-    for(G4int i=0; i<nProcesses; i++) {
+    for(G4int i=0; i<nProcesses; ++i) {
       if(p == scProcesses[i]) {
         add = false;
@@ -1616 +1651 @@
   if(add) {
     scProcesses.push_back(p);
-    nProcesses++;
+    ++nProcesses;
     if (1 < verboseLevel) { 
       G4cout << "### The process " << p->GetProcessName() 
@@ -1636 +1671 @@
       G4PhysicsVector* pv = (*p)[0];
       G4double emax = maxKinEnergyCSDA;
-      G4bool b;
       theDEDXAtMaxEnergy = new G4double [n];
 
-      for (size_t i=0; i<n; i++) {
+      for (size_t i=0; i<n; ++i) {
         pv = (*p)[i];
-        G4double dedx = pv->GetValue(emax, b);
+        G4double dedx = pv->Value(emax);
         theDEDXAtMaxEnergy[i] = dedx;
         //G4cout << "i= " << i << " emax(MeV)= " << emax/MeV<< " dedx= " 
@@ -1671 +1705 @@
     G4PhysicsVector* pv = (*p)[0];
     G4double emax = maxKinEnergyCSDA;
-    G4bool b;
     theRangeAtMaxEnergy = new G4double [n];
 
-    for (size_t i=0; i<n; i++) {
+    for (size_t i=0; i<n; ++i) {
       pv = (*p)[i];
-      G4double r2 = pv->GetValue(emax, b);
+      G4double r2 = pv->Value(emax);
       theRangeAtMaxEnergy[i] = r2;
       //G4cout << "i= " << i << " e2(MeV)= " << emax/MeV << " r2= " 
@@ -1744 +1777 @@
     theEnergyOfCrossSectionMax = new G4double [n];
     theCrossSectionMax = new G4double [n];
-    G4bool b;
-
-    for (size_t i=0; i<n; i++) {
+
+    for (size_t i=0; i<n; ++i) {
       pv = (*p)[i];
       emax = DBL_MAX;
@@ -1752 +1784 @@
       if(pv) {
         size_t nb = pv->GetVectorLength();
-        emax = pv->GetLowEdgeEnergy(nb);
-        for (size_t j=0; j<nb; j++) {
-          e = pv->GetLowEdgeEnergy(j);
-          s = pv->GetValue(e,b);
-          if(s > smax) {
-            smax = s;
-            emax = e;
+        if(nb > 0) {
+          for (size_t j=0; j<nb; ++j) {
+            e = pv->Energy(j);
+            s = (*pv)(j);
+            if(s > smax) {
+              smax = s;
+              emax = e;
+            }
           }
         }

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4VMscModel.cc,v 1.12 2009/05/10 19:36:19 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4VMscModel.cc,v 1.13 2009/07/20 17:32:47 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -60 +60 @@
   facrange(0.04),
   facgeom(2.5),
-  facsafety(0.25),
+  facsafety(0.3),
   skin(3.0),
   dtrl(0.05),

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4VMultipleScattering.cc,v 1.66 2009/05/27 11:55:02 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4VMultipleScattering.cc,v 1.77 2009/10/29 18:07:08 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -56 +56 @@
 // 27-10-07 Virtual functions moved to source (V.Ivanchenko)
 // 11-03-08 Set skin value does not effect step limit type (V.Ivanchenko)
+// 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
 //
 // Class Description:
@@ -98 +99 @@
   facgeom(2.5),
   latDisplasment(true),
+  isIon(false),
   currentParticle(0),
   currentCouple(0)
@@ -144 +146 @@
   if(p) p->SetParticleChange(pParticleChange);
 }
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4VEmModel* G4VMultipleScattering::GetModelByIndex(G4int idx, G4bool ver)
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4VMultipleScattering::SetModel(G4VMscModel* p, G4int index)
+{
+  G4int n = mscModels.size();
+  if(index >= n) { for(G4int i=n; i<=index; ++i) {mscModels.push_back(0);} }
+  mscModels[index] = p;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4VMscModel* G4VMultipleScattering::Model(G4int index)
+{
+  G4VMscModel* p = 0;
+  if(index >= 0 && index <  G4int(mscModels.size())) { p = mscModels[index]; }
+  return p;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4VEmModel* 
+G4VMultipleScattering::GetModelByIndex(G4int idx, G4bool ver) const
 {
   return modelManager->GetModel(idx, ver);
@@ -170 +190 @@
     size_t numOfCouples = theCoupleTable->GetTableSize();
 
-    for (size_t i=0; i<numOfCouples; i++) {
+    G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
+
+    G4PhysicsLogVector* aVector = 0;
+    G4PhysicsLogVector* bVector = 0;
+
+    for (size_t i=0; i<numOfCouples; ++i) {
 
       if (theLambdaTable->GetFlag(i)) {
@@ -176 +201 @@
         const G4MaterialCutsCouple* couple = 
           theCoupleTable->GetMaterialCutsCouple(i);
-        G4PhysicsVector* aVector = PhysicsVector(couple);
+        if(!bVector) {
+          aVector = static_cast<G4PhysicsLogVector*>(PhysicsVector(couple));
+          bVector = aVector;
+        } else {
+          aVector = new G4PhysicsLogVector(*bVector);
+        }       
+        //G4PhysicsVector* aVector = PhysicsVector(couple);
+        aVector->SetSpline(splineFlag);
         modelManager->FillLambdaVector(aVector, couple, false);
+        if(splineFlag) aVector->FillSecondDerivatives();
         G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
       }
@@ -212 +245 @@
        part.GetParticleSubType() == "generic") {
       firstParticle = G4GenericIon::GenericIon();
+      isIon = true;
     } else {
       firstParticle = &part;
+      if(part.GetParticleType() == "nucleus" || 
+         part.GetPDGMass() > GeV) {isIon = true;} 
     } 
-
+    // limitations for ions
+    if(isIon) {
+      SetStepLimitType(fMinimal);
+      SetLateralDisplasmentFlag(false);
+      SetBuildLambdaTable(false);
+    }
     currentParticle = &part;
   }
@@ -232 +273 @@
 
     InitialiseProcess(firstParticle);
-    if(buildLambdaTable)
+
+    // initialisation of models
+    G4int nmod = modelManager->NumberOfModels();
+    for(G4int i=0; i<nmod; ++i) {
+      G4VMscModel* msc = static_cast<G4VMscModel*>(modelManager->GetModel(i));
+      if(isIon) {
+        msc->SetStepLimitType(fMinimal);
+        msc->SetLateralDisplasmentFlag(false);
+        msc->SetRangeFactor(0.2);
+      } else {
+        msc->SetStepLimitType(StepLimitType());
+        msc->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
+        msc->SetSkin(Skin());
+        msc->SetRangeFactor(RangeFactor());
+        msc->SetGeomFactor(GeomFactor());
+      }
+      msc->SetPolarAngleLimit(polarAngleLimit);
+      if(msc->HighEnergyLimit() > maxKinEnergy) {
+        msc->SetHighEnergyLimit(maxKinEnergy);
+      }
+    }
+
+    modelManager->Initialise(firstParticle, G4Electron::Electron(), 
+                             10.0, verboseLevel);
+
+    // prepare tables
+    if(buildLambdaTable) {
       theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
-    const G4DataVector* theCuts = 
-      modelManager->Initialise(firstParticle, 
-                               G4Electron::Electron(), 
-                               10.0, verboseLevel);
-
-    if(2 < verboseLevel) G4cout << theCuts << G4endl;
-
+    }
   }
 }
@@ -277 +338 @@
                              G4double,
                              G4double currentMinimalStep,
-                             G4double& currentSafety,
+                             G4double&,
                              G4GPILSelection* selection)
 {
   // get Step limit proposed by the process
-  valueGPILSelectionMSC = NotCandidateForSelection;
-  G4double steplength = GetMscContinuousStepLimit(track,
-                                                  track.GetKineticEnergy(),
-                                                  currentMinimalStep,
-                                                  currentSafety);
-  // G4cout << "StepLimit= " << steplength << G4endl;
-  // set return value for G4GPILSelection
-  *selection = valueGPILSelectionMSC;
-  return  steplength;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4double G4VMultipleScattering::PostStepGetPhysicalInteractionLength(
-              const G4Track&, G4double, G4ForceCondition* condition)
-{
-  *condition = Forced;
-  return DBL_MAX;
+  *selection = NotCandidateForSelection;
+  G4double x = currentMinimalStep;
+  DefineMaterial(track.GetMaterialCutsCouple());
+  G4double ekin = track.GetKineticEnergy();
+  if(isIon) { ekin *= proton_mass_c2/track.GetDefinition()->GetPDGMass(); }
+  currentModel = static_cast<G4VMscModel*>(SelectModel(ekin));
+  if(x > 0.0 && ekin > 0.0 && currentModel->IsActive(ekin)) {
+    G4double tPathLength = 
+      currentModel->ComputeTruePathLengthLimit(track, theLambdaTable, x);
+    if (tPathLength < x) *selection = CandidateForSelection;
+    x = currentModel->ComputeGeomPathLength(tPathLength);
+    //  G4cout << "tPathLength= " << tPathLength
+    //         << " stepLimit= " << x
+    //        << " currentMinimalStep= " << currentMinimalStep<< G4endl;
+  }
+  return x;
 }
 
@@ -309 +368 @@
                                        G4double& currentSafety)
 {
-  return GetMscContinuousStepLimit(track,previousStepSize,currentMinimalStep,
-                                   currentSafety);
+  G4GPILSelection* selection = 0;
+  return AlongStepGetPhysicalInteractionLength(track,previousStepSize,currentMinimalStep,
+                                               currentSafety, selection);
 }
 
@@ -320 +380 @@
   *condition = Forced;
   return DBL_MAX;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4VParticleChange* G4VMultipleScattering::AlongStepDoIt(const G4Track&,
-                                                        const G4Step& step)
-{
-  fParticleChange.ProposeTrueStepLength(
-    currentModel->ComputeTrueStepLength(step.GetStepLength()));
-  return &fParticleChange;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4VParticleChange* G4VMultipleScattering::PostStepDoIt(const G4Track& track,
-                                                       const G4Step& step)
-{
-  fParticleChange.Initialize(track);
-  currentModel->SampleScattering(track.GetDynamicParticle(),
-                                 step.GetPostStepPoint()->GetSafety());
-  return &fParticleChange;
 }
 
@@ -384 +423 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-G4bool G4VMultipleScattering::RetrievePhysicsTable(const G4ParticleDefinition* part,
-                                                   const G4String& directory,
-                                                         G4bool ascii)
+G4bool 
+G4VMultipleScattering::RetrievePhysicsTable(const G4ParticleDefinition* part,
+                                            const G4String& directory,
+                                            G4bool ascii)
 {
   if(0 < verboseLevel) {
-//    G4cout << "========================================================" << G4endl;
     G4cout << "G4VMultipleScattering::RetrievePhysicsTable() for "
            << part->GetParticleName() << " and process "
@@ -401 +440 @@
 
   G4String filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
-  yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,filename,ascii);
+  yes = 
+    G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,filename,ascii);
   if ( yes ) {
     if (0 < verboseLevel) {
-        G4cout << "Lambda table for " << part->GetParticleName() << " is retrieved from <"
+        G4cout << "Lambda table for " << part->GetParticleName() 
+               << " is retrieved from <"
                << filename << ">"
                << G4endl;
@@ -410 +451 @@
     if((G4LossTableManager::Instance())->SplineFlag()) {
       size_t n = theLambdaTable->length();
-      for(size_t i=0; i<n; i++) {(* theLambdaTable)[i]->SetSpline(true);}
+      for(size_t i=0; i<n; ++i) {
+        if((* theLambdaTable)[i]) {
+          (* theLambdaTable)[i]->SetSpline(true);
+        }
+      }
     }
   } else {
     if (1 < verboseLevel) {
-        G4cout << "Lambda table for " << part->GetParticleName() << " in file <"
+        G4cout << "Lambda table for " << part->GetParticleName() 
+               << " in file <"
                << filename << "> is not exist"
                << G4endl;

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4ionEffectiveCharge.cc,v 1.24 2008/12/18 13:01:46 gunter Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4ionEffectiveCharge.cc,v 1.25 2009/10/29 16:57:39 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -105 +105 @@
   // Fast ions or hadrons
   G4double reducedEnergy = kineticEnergy * proton_mass_c2/mass ;
+
+  //G4cout << "e= " << reducedEnergy << " Zi= " << Zi << "  " << material->GetName() << G4endl;
+
   if( reducedEnergy > Zi*energyHighLimit || Zi < 1.5 || !material) return charge;
 
   G4double z    = material->GetIonisation()->GetZeffective();
-  //  reducedEnergy = std::max(reducedEnergy,energyLowLimit);
+  reducedEnergy = std::max(reducedEnergy,energyLowLimit);
 
   // Helium ion case