Changeset 961 for trunk/source/processes/electromagnetic/utils/src/G4VEmProcess.cc

Timestamp:

Apr 6, 2009, 12:21:12 PM (15 years ago)

Author:

garnier

Message:

update processes

File:

• trunk/source/processes/electromagnetic/utils/src/G4VEmProcess.cc (modified) (12 diffs)

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4VEmProcess.cc,v 1.48 2007/10/29 08:38:58 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4VEmProcess.cc,v 1.62 2009/02/19 09:57:36 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -82 +82 @@
 
 G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type):
-                      G4VDiscreteProcess(name, type),
-  selectedModel(0),                   
+  G4VDiscreteProcess(name, type),
+  secondaryParticle(0),
+  buildLambdaTable(true),
   theLambdaTable(0),
   theEnergyOfCrossSectionMax(0),
   theCrossSectionMax(0),
-  particle(0),
-  secondaryParticle(0),
-  nLambdaBins(90),
-  lambdaFactor(0.8),
-  currentCouple(0),
   integral(false),
-  buildLambdaTable(true),
   applyCuts(false),
   startFromNull(true),
-  nRegions(0)
+  useDeexcitation(false),
+  nDERegions(0),
+  idxDERegions(0),
+  currentModel(0),
+  particle(0),
+  currentCouple(0)
 {
   SetVerboseLevel(1);
+
+  // Size of tables assuming spline
   minKinEnergy = 0.1*keV;
-  maxKinEnergy = 100.0*GeV;
+  maxKinEnergy = 100.0*TeV;
+  nLambdaBins  = 84;
+
+  // default lambda factor
+  lambdaFactor  = 0.8;
+
+  // default limit on polar angle
+  polarAngleLimit = 0.0;
+
+  // particle types
   theGamma     = G4Gamma::Gamma();
   theElectron  = G4Electron::Electron();
@@ -126 +137 @@
   delete modelManager;
   (G4LossTableManager::Instance())->DeRegister(this);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4VEmProcess::Clear()
+{
+  delete [] theEnergyOfCrossSectionMax;
+  delete [] theCrossSectionMax;
+  delete [] idxDERegions;
+  theEnergyOfCrossSectionMax = 0;
+  theCrossSectionMax = 0;
+  idxDERegions = 0;
+  currentCouple = 0;
+  preStepLambda = 0.0;
+  mfpKinEnergy  = DBL_MAX;
+  deRegions.clear();
+  nDERegions = 0;
 }
 
@@ -153 +181 @@
       theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
   }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-void G4VEmProcess::Clear()
-{
-  if(theEnergyOfCrossSectionMax) delete [] theEnergyOfCrossSectionMax;
-  if(theCrossSectionMax) delete [] theCrossSectionMax;
-  theEnergyOfCrossSectionMax = 0;
-  theCrossSectionMax = 0;
-  currentCouple = 0;
-  preStepLambda = 0.0;
-  mfpKinEnergy  = DBL_MAX;
+  // Sub Cutoff and Deexcitation
+  if (nDERegions>0) {
+
+    const G4ProductionCutsTable* theCoupleTable=
+          G4ProductionCutsTable::GetProductionCutsTable();
+    size_t numOfCouples = theCoupleTable->GetTableSize();
+
+    idxDERegions = new G4bool[numOfCouples];
+
+    for (size_t j=0; j<numOfCouples; j++) {
+
+      const G4MaterialCutsCouple* couple =
+        theCoupleTable->GetMaterialCutsCouple(j);
+      const G4ProductionCuts* pcuts = couple->GetProductionCuts();
+      G4bool reg = false;
+      for(G4int i=0; i<nDERegions; i++) {
+        if(deRegions[i]) {
+          if(pcuts == deRegions[i]->GetProductionCuts()) reg = true;
+        }
+      }
+      idxDERegions[j] = reg;
+    }
+  }
+  if (1 < verboseLevel && nDERegions>0) {
+    G4cout << " Deexcitation is activated for regions: " << G4endl;
+    for (G4int i=0; i<nDERegions; i++) {
+      const G4Region* r = deRegions[i];
+      G4cout << "           " << r->GetName() << G4endl;
+    }
+  }
 }
 
@@ -228 +273 @@
       G4cout << *theLambdaTable << G4endl;
     }
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4VEmProcess::PrintInfoDefinition()
+{
+  if(verboseLevel > 0) {
+    G4cout << G4endl << GetProcessName() << ":   for  "
+           << particle->GetParticleName();
+    if(integral) G4cout << ", integral: 1 ";
+    if(applyCuts) G4cout << ", applyCuts: 1 ";
+    G4cout << "    SubType= " << GetProcessSubType() << G4endl;
+    if(buildLambdaTable) {
+      G4cout << "      Lambda tables from "
+             << G4BestUnit(minKinEnergy,"Energy") 
+             << " to "
+             << G4BestUnit(maxKinEnergy,"Energy")
+             << " in " << nLambdaBins << " bins, spline: " 
+             << (G4LossTableManager::Instance())->SplineFlag()
+             << G4endl;
+    }
+    PrintInfo();
+    modelManager->DumpModelList(verboseLevel);
+  }
+
+  if(verboseLevel > 2 && buildLambdaTable) {
+    G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
+    if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
   }
 }
@@ -295 +369 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
-                                       G4double,
-                                       G4ForceCondition* condition)
-{
-  *condition = NotForced;
-  return G4VEmProcess::MeanFreePath(track);
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
                                               const G4Step&)
@@ -333 +397 @@
   }
 
-  G4VEmModel* currentModel = SelectModel(finalT);
-
+  SelectModel(finalT);
+  if(useDeexcitation) {
+    currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);
+  }
   /*  
   if(0 < verboseLevel) {
@@ -373 +439 @@
 
         } else if (p == thePositron) {
-          if (e < (*theCutsPositron)[currentMaterialIndex]) {
+          if (electron_mass_c2 < (*theCutsGamma)[currentMaterialIndex] &&
+              e < (*theCutsPositron)[currentMaterialIndex]) {
             good = false;
             e += 2.0*electron_mass_c2;
@@ -394 +461 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-void G4VEmProcess::PrintInfoDefinition()
-{
-  if(verboseLevel > 0) {
-    G4cout << G4endl << GetProcessName() << ": " ;
-    PrintInfo();
-    if(integral) {
-      G4cout << "      Integral mode is used  "<< G4endl;
-    }
-  }
-
-  if (!buildLambdaTable)  return;
-  
-  if(verboseLevel > 0) {
-    G4cout << "      tables are built for  "
-           << particle->GetParticleName()
-           << G4endl
-           << "      Lambda tables from "
-           << G4BestUnit(minKinEnergy,"Energy") 
-           << " to "
-           << G4BestUnit(maxKinEnergy,"Energy")
-           << " in " << nLambdaBins << " bins."
-           << G4endl;
-  }
-
-  if(verboseLevel > 1) {
-    G4cout << "Tables are built for " << particle->GetParticleName()
-           << G4endl;
-
-  if(verboseLevel > 2) {
-    G4cout << "LambdaTable address= " << theLambdaTable << G4endl;
-    if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
-    }
-  }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4double G4VEmProcess::MicroscopicCrossSection(G4double kineticEnergy,
-                                         const G4MaterialCutsCouple* couple)
-{
-  // Cross section per atom is calculated
-  DefineMaterial(couple);
-  G4double cross = 0.0;
-  G4bool b;
-  if(theLambdaTable) {
-    cross = (((*theLambdaTable)[currentMaterialIndex])->
-                           GetValue(kineticEnergy, b));
-
-    cross /= currentMaterial->GetTotNbOfAtomsPerVolume();
-  } else {
-    G4VEmModel* model = SelectModel(kineticEnergy);
-    cross = 
-      model->CrossSectionPerVolume(currentMaterial,particle,kineticEnergy);
-  }
-
-  return cross;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
                                        const G4String& directory,
@@ -508 +515 @@
              << G4endl;
     }
+    if((G4LossTableManager::Instance())->SplineFlag()) {
+      size_t n = theLambdaTable->length();
+      for(size_t i=0; i<n; i++) {(* theLambdaTable)[i]->SetSpline(true);}
+    }
   } else {
     if (1 < verboseLevel) {
@@ -517 +528 @@
 
   return yes;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4VEmProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
+{
+  G4RegionStore* regionStore = G4RegionStore::GetInstance();
+  const G4Region* reg = r;
+  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
+
+  // the region is in the list
+  if (nDERegions) {
+    for (G4int i=0; i<nDERegions; i++) {
+      if (reg == deRegions[i]) {
+        if(!val) deRegions[i] = 0;
+        return;
+      }
+    }
+  }
+
+  // new region 
+  if(val) {
+    useDeexcitation = true;
+    deRegions.push_back(reg);
+    nDERegions++;
+  } else {
+    useDeexcitation = false;
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
+                                             const G4MaterialCutsCouple* couple)
+{
+  // Cross section per atom is calculated
+  DefineMaterial(couple);
+  G4double cross = 0.0;
+  G4bool b;
+  if(theLambdaTable) {
+    cross = (((*theLambdaTable)[currentMaterialIndex])->
+                           GetValue(kineticEnergy, b));
+  } else {
+    SelectModel(kineticEnergy);
+    cross = currentModel->CrossSectionPerVolume(currentMaterial,
+                                                particle,kineticEnergy);
+  }
+
+  return cross;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
+                                       G4double,
+                                       G4ForceCondition* condition)
+{
+  *condition = NotForced;
+  return G4VEmProcess::MeanFreePath(track);
 }
 
@@ -568 +638 @@
   G4PhysicsVector* v = 
     new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins);
+  v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
   return v;
 }