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

Timestamp:

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

Author:

garnier

Message:

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

File:

• trunk/source/processes/electromagnetic/utils/src/G4VEnergyLossProcess.cc (modified) (41 diffs)

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4VEnergyLossProcess.cc,v 1.158 2009/10/29 18:07:08 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4VEnergyLossProcess.cc,v 1.167 2010/05/26 10:41:34 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -336 +336 @@
 {
   modelManager->AddEmModel(order, p, fluc, region);
-  if(p) p->SetParticleChange(pParticleChange, fluc);
+  if(p) { p->SetParticleChange(pParticleChange, fluc); }
 }
 
@@ -406 +406 @@
     particle = ∂
     if(part.GetParticleType() == "nucleus") {
-      if(!theGenericIon) theGenericIon = G4GenericIon::GenericIon();
+      if(!theGenericIon) { theGenericIon = G4GenericIon::GenericIon(); }
       if(particle == theGenericIon) { isIon = true; }
       else if(part.GetPDGCharge() > eplus) {
@@ -426 +426 @@
       lManager->RegisterExtraParticle(&part, this);
     }
+    if(1 < verboseLevel) {
+      G4cout << "### G4VEnergyLossProcess::PreparePhysicsTable() end for "
+             << part.GetParticleName() << G4endl;
+    }
     return;
   }
 
   Clean();
-  lManager->EmConfigurator()->AddModels();
+  lManager->PreparePhysicsTable(&part, this);
 
   // Base particle and set of models can be defined here
@@ -502 +506 @@
         G4bool reg = false;
         for(G4int i=0; i<nSCoffRegions; ++i) {
-          if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = true;
+          if( pcuts == scoffRegions[i]->GetProductionCuts()) { reg = true; }
         }
         idxSCoffRegions[j] = reg;
@@ -509 +513 @@
         G4bool reg = false;
         for(G4int i=0; i<nDERegions; ++i) {
-          if( pcuts == deRegions[i]->GetProductionCuts()) reg = true;
+          if( pcuts == deRegions[i]->GetProductionCuts()) { reg = true; }
         }
         idxDERegions[j] = reg;
@@ -515 +519 @@
     }
   }
-
-  lManager->EnergyLossProcessIsInitialised(particle, this);
 
   if (1 < verboseLevel) {
@@ -568 +570 @@
 
   // Added tracking cut to avoid tracking artifacts
-  if(isIonisation) fParticleChange.SetLowEnergyLimit(lowestKinEnergy);
+  if(isIonisation) { fParticleChange.SetLowEnergyLimit(lowestKinEnergy); }
 
   if(1 < verboseLevel) {
@@ -574 +576 @@
            << GetProcessName()
            << " and particle " << part.GetParticleName();
-    if(isIonisation) G4cout << "  isIonisation  flag = 1";
+    if(isIonisation) { G4cout << "  isIonisation  flag = 1"; }
     G4cout << G4endl;
   }
@@ -590 +592 @@
   }
   G4PhysicsTable* table = 0;
-  G4double emin = minKinEnergy;
   G4double emax = maxKinEnergy;
   G4int bin = nBins;
@@ -619 +620 @@
     G4cout << numOfCouples << " materials"
            << " minKinEnergy= " << minKinEnergy
-           << " maxKinEnergy= " << maxKinEnergy
+           << " maxKinEnergy= " << emax
+           << " nbin= " << bin
            << " EmTableType= " << tType
            << " table= " << table
@@ -642 +644 @@
         theCoupleTable->GetMaterialCutsCouple(i);
       if(!bVector) {
-        aVector = new G4PhysicsLogVector(emin, emax, bin);
+        aVector = new G4PhysicsLogVector(minKinEnergy, 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();
+      if(splineFlag) { aVector->FillSecondDerivatives(); }
 
       // Insert vector for this material into the table
@@ -701 +702 @@
 
   G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
+  G4PhysicsLogVector* aVector = 0;
+  G4PhysicsLogVector* bVector = 0;
 
   for(size_t i=0; i<numOfCouples; ++i) {
@@ -709 +712 @@
       const G4MaterialCutsCouple* couple = 
         theCoupleTable->GetMaterialCutsCouple(i);
-      G4double cut = (*theCuts)[i];
-      if(fSubRestricted == tType) cut = (*theSubCuts)[i]; 
-      G4PhysicsVector* aVector = LambdaPhysicsVector(couple, cut);
+      if(!bVector) {
+        aVector = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nBins);
+        bVector = aVector;
+      } else {
+        aVector = new G4PhysicsLogVector(*bVector);
+      }
       aVector->SetSpline(splineFlag);
 
       modelManager->FillLambdaVector(aVector, couple, true, tType);
-      if(splineFlag) aVector->FillSecondDerivatives();
+      if(splineFlag) { aVector->FillSecondDerivatives(); }
 
       // Insert vector for this material into the table
@@ -881 +887 @@
     if(x > finalRange && y < currentMinStep) { 
       x = y + finalRange*(1.0 - dRoverRange)*(2.0 - finalRange/fRange);
-    } else if (rndmStepFlag) {x = SampleRange();}
+    } else if (rndmStepFlag) { x = SampleRange(); }
     //G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
     //  <<" range= "<<fRange <<" cMinSt="<<currentMinStep
@@ -912 +918 @@
   preStepScaledEnergy = preStepKinEnergy*massRatio;
   SelectModel(preStepScaledEnergy);
-  if(!currentModel->IsActive(preStepScaledEnergy)) return x;
-
-  if(isIon) {
-    chargeSqRatio = 
-      currentModel->GetChargeSquareRatio(currPart,currentMaterial,preStepKinEnergy);
+  if(!currentModel->IsActive(preStepScaledEnergy)) { return x; }
+
+  if(isIon) { 
+    chargeSqRatio = currentModel->ChargeSquareRatio(track);
     reduceFactor  = 1.0/(chargeSqRatio*massRatio);
   }
   //G4cout << "q2= " << chargeSqRatio << " massRatio= " << massRatio << G4endl; 
   // initialisation for sampling of the interaction length 
-  if(previousStepSize <= DBL_MIN) theNumberOfInteractionLengthLeft = -1.0;
-  if(theNumberOfInteractionLengthLeft < 0.0) mfpKinEnergy = DBL_MAX;
+  if(previousStepSize <= DBL_MIN) { theNumberOfInteractionLengthLeft = -1.0; }
+  if(theNumberOfInteractionLengthLeft < 0.0) { mfpKinEnergy = DBL_MAX; }
 
   // compute mean free path
   if(preStepScaledEnergy < mfpKinEnergy) {
-    if (integral) ComputeLambdaForScaledEnergy(preStepScaledEnergy);
-    else  preStepLambda = GetLambdaForScaledEnergy(preStepScaledEnergy);
-    if(preStepLambda <= DBL_MIN) mfpKinEnergy = 0.0;
+    if (integral) { ComputeLambdaForScaledEnergy(preStepScaledEnergy); }
+    else  { preStepLambda = GetLambdaForScaledEnergy(preStepScaledEnergy); }
+    if(preStepLambda <= DBL_MIN) { mfpKinEnergy = 0.0; }
   }
 
@@ -940 +945 @@
       // subtract NumberOfInteractionLengthLeft
       SubtractNumberOfInteractionLengthLeft(previousStepSize);
-      if(theNumberOfInteractionLengthLeft < 0.)
+      if(theNumberOfInteractionLengthLeft < 0.) {
         theNumberOfInteractionLengthLeft = perMillion;
+      }
     }
 
@@ -966 +972 @@
       // subtract NumberOfInteractionLengthLeft
       SubtractNumberOfInteractionLengthLeft(previousStepSize);
-      if(theNumberOfInteractionLengthLeft < 0.)
+      if(theNumberOfInteractionLengthLeft < 0.) {
         theNumberOfInteractionLengthLeft = perMillion;
+      }
     }
     currentInteractionLength = DBL_MAX;
@@ -987 +994 @@
   // Get the actual (true) Step length
   G4double length = step.GetStepLength();
-  if(length <= DBL_MIN) return &fParticleChange;
+  if(length <= DBL_MIN) { return &fParticleChange; }
   G4double eloss  = 0.0;
-  G4double esecdep = 0.0;
  
   /*  
@@ -1071 +1077 @@
 
       // Check boundary
-      if(prePoint->GetStepStatus() == fGeomBoundary) yes = true;
+      if(prePoint->GetStepStatus() == fGeomBoundary) { yes = true; }
 
       // Check PrePoint
@@ -1083 +1089 @@
         }
 
-        if(preSafety < rcut) yes = true;
+        if(preSafety < rcut) { yes = true; }
 
         // Check PostPoint
@@ -1091 +1097 @@
             postSafety = 
               safetyHelper->ComputeSafety(step.GetPostStepPoint()->GetPosition());
-            if(postSafety < rcut) yes = true;
+            if(postSafety < rcut) { yes = true; }
           }
         }
@@ -1103 +1109 @@
         scTracks.clear();
         SampleSubCutSecondaries(scTracks, step, 
-                                currentModel,currentMaterialIndex, 
-                                esecdep);
+                                currentModel,currentMaterialIndex);
         // add bremsstrahlung sampling
         /*
@@ -1112 +1117 @@
                 scTracks, step, (scProcesses[i])->
                 SelectModelForMaterial(preStepKinEnergy, currentMaterialIndex),
-                currentMaterialIndex,esecdep);
+                currentMaterialIndex);
           }
         } 
@@ -1123 +1128 @@
             G4Track* t = scTracks[i];
             G4double e = t->GetKineticEnergy();
-            if (t->GetDefinition() == thePositron) e += 2.0*electron_mass_c2;
+            if (t->GetDefinition() == thePositron) { e += 2.0*electron_mass_c2; }
             esec += e;
             pParticleChange->AddSecondary(t);
@@ -1133 +1138 @@
 
   // Corrections, which cannot be tabulated
-  currentModel->CorrectionsAlongStep(currentCouple, dynParticle, 
-                                     eloss, esecdep, length);
+  if(isIon) {
+    G4double eadd = 0.0;
+    currentModel->CorrectionsAlongStep(currentCouple, dynParticle, 
+                                       eloss, eadd, length);
+  }
 
   // Sample fluctuations
@@ -1140 +1148 @@
     G4VEmFluctuationModel* fluc = currentModel->GetModelOfFluctuations();
     if(fluc && 
-      (eloss + esec + esecdep + lowestKinEnergy) < preStepKinEnergy) {
+      (eloss + esec + lowestKinEnergy) < preStepKinEnergy) {
 
       G4double tmax = 
@@ -1158 +1166 @@
     }
   }
-  // add low-energy subcutoff particles
-  eloss += esecdep;
-  if(eloss < 0.0) eloss = 0.0;
 
   // deexcitation
-  else if (useDeexcitation) {
-    if(idxDERegions[currentMaterialIndex]) {
+  if (useDeexcitation) {
+    if(idxDERegions[currentMaterialIndex] && eloss > 0.0) {
       currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
-      if(eloss < 0.0) eloss = 0.0;
+      if(eloss < 0.0) { eloss = 0.0; }
     }
   }
@@ -1203 +1208 @@
        const G4Step& step, 
        G4VEmModel* model,
-       G4int idx,
-       G4double& /*extraEdep*/) 
+       G4int idx) 
 {
   // Fast check weather subcutoff can work
   G4double subcut = (*theSubCuts)[idx];
   G4double cut = (*theCuts)[idx];
-  if(cut <= subcut) return;
+  if(cut <= subcut) { return; }
 
   const G4Track* track = step.GetTrack();
@@ -1218 +1222 @@
 
   // negligible probability to get any interaction
-  if(length*cross < perMillion) return;
+  if(length*cross < perMillion) { return; }
   /*      
   if(-1 < verboseLevel) 
@@ -1269 +1273 @@
       if(addSec) {
         G4Track* t = new G4Track((*it), pretime + fragment*dt, r);
-        //G4Track* t = new G4Track((*it), pretime, r);
         t->SetTouchableHandle(track->GetTouchableHandle());
         tracks.push_back(t);
@@ -1296 +1299 @@
   G4double postStepScaledEnergy = finalT*massRatio;
 
-  if(!currentModel->IsActive(postStepScaledEnergy)) return &fParticleChange;
+  if(!currentModel->IsActive(postStepScaledEnergy)) { return &fParticleChange; }
   /*
   if(-1 < verboseLevel) {
@@ -1506 +1509 @@
                                     G4bool mandatory)
 {
-  G4bool res = true;
+  G4bool isRetrieved = false;
   G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
-  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) {
-        if((*aTable)[i]) {
-          (*aTable)[i]->SetSpline(true);
+  if(aTable) {
+    if(aTable->ExistPhysicsTable(filename)) {
+      if(G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii)) {
+        isRetrieved = true;
+        if((G4LossTableManager::Instance())->SplineFlag()) {
+          size_t n = aTable->length();
+          for(size_t i=0; i<n; ++i) {
+            if((*aTable)[i]) { (*aTable)[i]->SetSpline(true); }
+          }
         }
-      }
-    }
-  }
-  if(yes) {
-    if (0 < verboseLevel) {
-      G4cout << tname << " table for " << part->GetParticleName() 
-             << " is Retrieved from <" << filename << ">"
-             << G4endl;
-    }
-  } else {
-    if(mandatory) res = false;
-    if(mandatory || 1 < verboseLevel) {
+        if (0 < verboseLevel) {
+          G4cout << tname << " table for " << part->GetParticleName() 
+                 << " is Retrieved from <" << filename << ">"
+                 << G4endl;
+        }
+      }
+    }
+  }
+  if(mandatory && !isRetrieved) {
+    if(0 < verboseLevel) {
       G4cout << tname << " table for " << part->GetParticleName() 
              << " from file <"
@@ -1536 +1536 @@
              << G4endl;
     }
-  }
-  return res;
+    return false;
+  }
+  return true;
 }
 
@@ -1574 +1575 @@
                                                 (*theCuts)[currentMaterialIndex]);
   }
-
+  if(cross < 0.0) { cross = 0.0; }
   return cross;
 }
@@ -1585 +1586 @@
   preStepLambda = GetLambdaForScaledEnergy(track.GetKineticEnergy()*massRatio);
   G4double x = DBL_MAX;
-  if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
+  if(DBL_MIN < preStepLambda) { x = 1.0/preStepLambda; }
   return x;
 }
@@ -1622 +1623 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(
-                 const G4MaterialCutsCouple* couple, G4double cut)
-{
+G4PhysicsVector* 
+G4VEnergyLossProcess::LambdaPhysicsVector(const G4MaterialCutsCouple* /*couple*/, 
+                                          G4double /*cut*/)
+{
+  /*
   G4double tmin = 
     std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
              minKinEnergy);
-  if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy;
+  if(tmin >= maxKinEnergy) { tmin = 0.5*maxKinEnergy; }
   G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
+  */
+  G4PhysicsVector* v = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nBins);
   v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
   return v;
@@ -1640 +1645 @@
 {
   G4bool add = true;
-  if(p->GetProcessName() != "eBrem") add = false;
+  if(p->GetProcessName() != "eBrem") { add = false; }
   if(add && nProcesses > 0) {
     for(G4int i=0; i<nProcesses; ++i) {
@@ -1699 +1704 @@
 void G4VEnergyLossProcess::SetCSDARangeTable(G4PhysicsTable* p)
 {
-  if(theCSDARangeTable != p) theCSDARangeTable = p;
+  if(theCSDARangeTable != p) { theCSDARangeTable = p; }
 
   if(p) {
@@ -1768 +1773 @@
            << " and process " << GetProcessName() << G4endl;
   }
-  if(theLambdaTable != p) theLambdaTable = p;
+  if(theLambdaTable != p) { theLambdaTable = p; }
   tablesAreBuilt = true;
 
@@ -1822 +1827 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+const G4Element* G4VEnergyLossProcess::GetCurrentElement() const
+{
+  const G4Element* elm = 0;
+  if(currentModel) { elm = currentModel->GetCurrentElement(); }
+  return elm;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+