Changeset 961 for trunk/source/processes/electromagnetic/muons/src/G4EnergyLossForExtrapolator.cc

Timestamp:

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

Author:

garnier

Message:

update processes

File:

• trunk/source/processes/electromagnetic/muons/src/G4EnergyLossForExtrapolator.cc (modified) (11 diffs)

trunk/source/processes/electromagnetic/muons/src/G4EnergyLossForExtrapolator.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4EnergyLossForExtrapolator.cc,v 1.13 2007/07/28 13:44:25 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4EnergyLossForExtrapolator.cc,v 1.18 2008/11/13 14:14:07 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 //---------------------------------------------------------------------------
@@ -67 +67 @@
 #include "G4MuBremsstrahlungModel.hh"
 #include "G4ProductionCuts.hh"
+#include "G4LossTableManager.hh"
+#include "G4WentzelVIModel.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -88 +90 @@
   delete invRangePositron;
   delete invRangeProton;
+  delete mscElectron;
   delete cuts;
 }
@@ -100 +103 @@
   if(!isInitialised) Initialisation();
   G4double kinEnergyFinal = kinEnergy;
-  if(mat && part) {
-    G4double step = ComputeTrueStep(mat,part,kinEnergy,stepLength);
+  if(SetupKinematics(part, mat, kinEnergy)) {
+    G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
     G4double r  = ComputeRange(kinEnergy,part);
     if(r <= step) {
@@ -125 +128 @@
   G4double kinEnergyFinal = kinEnergy;
 
-  if(mat && part) {
-    G4double step = ComputeTrueStep(mat,part,kinEnergy,stepLength);
+  if(SetupKinematics(part, mat, kinEnergy)) {
+    G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
     G4double r  = ComputeRange(kinEnergy,part);
 
@@ -141 +144 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-G4double G4EnergyLossForExtrapolator::ComputeTrueStep(const G4Material* mat, 
-                                                      const G4ParticleDefinition* part,
-                                                      G4double kinEnergy, G4double stepLength)
-{
+G4double G4EnergyLossForExtrapolator::TrueStepLength(G4double kinEnergy, 
+                                                     G4double stepLength,
+                                                     const G4Material* mat, 
+                                                     const G4ParticleDefinition* part)
+{
+  G4double res = stepLength;
+  if(!isInitialised) Initialisation();
+  if(SetupKinematics(part, mat, kinEnergy)) {
+    if(part == electron || part == positron) {
+      G4double x = stepLength*ComputeValue(kinEnergy, mscElectron);
+      if(x < 0.2) res *= (1.0 + 0.5*x + x*x/3.0);
+      else if(x < 0.9999) res = -std::log(1.0 - x)*stepLength/x;
+      else res = ComputeRange(kinEnergy,part);
+    
+    } else {
+      res = ComputeTrueStep(mat,part,kinEnergy,stepLength);
+    }
+  }
+  return res;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4bool G4EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part, 
+                                                    const G4Material* mat, 
+                                                    G4double kinEnergy)
+{
+  if(!part || !mat || kinEnergy < keV) return false;
   if(!isInitialised) Initialisation();
   G4bool flag = false;
@@ -182 +209 @@
     if(tmax > maxEnergyTransfer) tmax = maxEnergyTransfer;
   }
-  G4double theta = ComputeScatteringAngle(stepLength);
-  return stepLength*std::sqrt(1.0 + 0.625*theta*theta);
+  return true;
 } 
 
@@ -231 +257 @@
   invRangeMuon     = PrepareTable();
   invRangeProton   = PrepareTable();
+  mscElectron      = PrepareTable();
 
   G4LossTableBuilder builder; 
@@ -262 +289 @@
   builder.BuildInverseRangeTable(rangeProton, invRangeProton);  
 
+  ComputeTrasportXS(electron, mscElectron);
 }
 
@@ -273 +301 @@
 
     G4PhysicsVector* v = new G4PhysicsLogVector(emin, emax, nbins);
+    v->SetSpline(G4LossTableManager::Instance()->SplineFlag());
     table->push_back(v);
   }
@@ -488 +517 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+void G4EnergyLossForExtrapolator::ComputeTrasportXS(const G4ParticleDefinition* part, 
+                                                    G4PhysicsTable* table)
+{
+  G4DataVector v;
+  G4WentzelVIModel* msc = new G4WentzelVIModel();
+  msc->SetPolarAngleLimit(CLHEP::pi);
+  msc->Initialise(part, v);
+
+  mass    = part->GetPDGMass();
+  charge2 = 1.0;
+  currentParticle = part;
+
+  const G4MaterialTable* mtable = G4Material::GetMaterialTable();
+
+  if(0<verbose) {
+    G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName() 
+           << G4endl;
+  }
+ 
+  for(G4int i=0; i<nmat; i++) {  
+
+    const G4Material* mat = (*mtable)[i];
+    if(1<verbose)
+      G4cout << "i= " << i << "  mat= " << mat->GetName() << G4endl;
+    G4PhysicsVector* aVector = (*table)[i];
+    for(G4int j=0; j<nbins; j++) {
+        
+       G4double e = aVector->GetLowEdgeEnergy(j);
+       G4double xs = msc->CrossSectionPerVolume(mat,part,e);
+       aVector->PutValue(j,xs);
+       if(1<verbose) {
+         G4cout << "j= " << j << "  e(MeV)= " << e/MeV  
+                << " xs(1/mm)= " << xs*mm << G4endl;
+       }
+    }
+  }
+  delete msc;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+