Changeset 961 for trunk/source/processes/electromagnetic/standard/src/G4hIonisation.cc

Timestamp:

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

Author:

garnier

Message:

update processes

File:

• trunk/source/processes/electromagnetic/standard/src/G4hIonisation.cc (modified) (8 diffs)

trunk/source/processes/electromagnetic/standard/src/G4hIonisation.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4hIonisation.cc,v 1.70 2008/01/14 11:59:45 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4hIonisation.cc,v 1.82 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -78 +78 @@
 //          positive from pi+ and p (VI)
 // 14-01-07 use SetEmModel() and SetFluctModel() from G4VEnergyLossProcess (mma)
+// 12-09-08 Removed CorrectionsAlongStep (VI)
 //
 // -------------------------------------------------------------------
@@ -90 +91 @@
 #include "G4BraggModel.hh"
 #include "G4BetheBlochModel.hh"
+#include "G4IonFluctuations.hh"
 #include "G4UniversalFluctuation.hh"
 #include "G4BohrFluctuations.hh"
@@ -95 +97 @@
 #include "G4PionPlus.hh"
 #include "G4PionMinus.hh"
-#include "G4LossTableManager.hh"
+#include "G4KaonPlus.hh"
+#include "G4KaonMinus.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -103 +106 @@
 G4hIonisation::G4hIonisation(const G4String& name)
   : G4VEnergyLossProcess(name),
-    theParticle(0),
-    theBaseParticle(0),
-    isInitialised(false)
-{
-  SetStepFunction(0.2, 1*mm);
-  SetIntegral(true);
-  SetVerboseLevel(1);
+    isInitialised(false),
+    nuclearStopping(true)
+{
+  //  SetStepFunction(0.2, 1.0*mm);
+  //SetIntegral(true);
+  //SetVerboseLevel(1);
+  SetProcessSubType(fIonisation);
   mass = 0.0;
   ratio = 0.0;
-  corr = G4LossTableManager::Instance()->EmCorrections();  
-  nuclearStopping = true;
 }
 
@@ -120 +121 @@
 G4hIonisation::~G4hIonisation()
 {}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4bool G4hIonisation::IsApplicable(const G4ParticleDefinition& p)
+{
+  return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV &&
+         !p.IsShortLived());
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4hIonisation::MinPrimaryEnergy(const G4ParticleDefinition*,
+                                         const G4Material*,
+                                         G4double cut)
+{
+  G4double x = 0.5*cut/electron_mass_c2;
+  G4double g = x*ratio + std::sqrt((1. + x)*(1. + x*ratio*ratio));
+  return mass*(g - 1.0);
+}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....  
@@ -127 +147 @@
                     const G4ParticleDefinition* bpart)
 {
-  if(isInitialised) return;
-
-  theParticle = part;
-
-  G4String pname = part->GetParticleName();
-
-  // standard base particles
-  if(part == bpart || pname == "proton" ||
-     pname == "anti_proton" || pname == "pi+" || pname == "pi-" ) 
-    theBaseParticle = 0;
-
-  // select base particle 
-  else if(bpart == 0) {
-
-    if(pname == "kaon+")      theBaseParticle = G4PionPlus::PionPlus();
-    else if(pname == "kaon-") theBaseParticle = G4PionMinus::PionMinus();
-    else if(part->GetPDGCharge() > 0.0) theBaseParticle = G4Proton::Proton();
-    else theBaseParticle = G4AntiProton::AntiProton();
-
-  } else theBaseParticle = bpart;
-
-  SetBaseParticle(theBaseParticle);
-  SetSecondaryParticle(G4Electron::Electron());
-  mass  = theParticle->GetPDGMass();
-  ratio = electron_mass_c2/mass;
-  massratio = 1.0;
-  if(theBaseParticle) massratio = theBaseParticle->GetPDGMass()/mass; 
-
-  if (!EmModel(1)) SetEmModel(new G4BraggModel(),1);
-  EmModel(1)->SetLowEnergyLimit(100*eV);
-  eth = 2.0*MeV*mass/proton_mass_c2;
-  ethnuc = eth*50.0;
-  EmModel(1)->SetHighEnergyLimit(eth);
-  if (!FluctModel()) SetFluctModel(new G4UniversalFluctuation());
-  AddEmModel(1, EmModel(1), FluctModel());
-
-  if (!EmModel(2)) SetEmModel(new G4BetheBlochModel(),2);  
-  EmModel(2)->SetLowEnergyLimit(eth);
-  EmModel(2)->SetHighEnergyLimit(100*TeV);
-  AddEmModel(2, EmModel(2), FluctModel());  
-
-  isInitialised = true;
+  if(!isInitialised) {
+
+    const G4ParticleDefinition* theBaseParticle = 0;
+    G4String pname = part->GetParticleName();
+
+    // standard base particles
+    if(part == bpart || pname == "proton" ||
+       pname == "anti_proton" || 
+       pname == "pi+" || pname == "pi-" || 
+       pname == "kaon+" || pname == "kaon-") 
+      { 
+        theBaseParticle = 0;
+      }
+    // select base particle 
+    else if(bpart == 0) {
+
+      if(part->GetPDGSpin() == 0.0)
+        if(part->GetPDGCharge() > 0.0 ) {
+          theBaseParticle = G4KaonPlus::KaonPlus();
+        } else {
+          theBaseParticle = G4KaonMinus::KaonMinus();
+        }
+      else if(part->GetPDGCharge() > 0.0) {
+        theBaseParticle = G4Proton::Proton();
+      } else {
+        theBaseParticle = G4AntiProton::AntiProton();
+      }
+      // base particle defined by interface
+    } else { 
+      theBaseParticle = bpart;
+    }
+    SetBaseParticle(theBaseParticle);
+    SetSecondaryParticle(G4Electron::Electron());
+
+    mass  = part->GetPDGMass();
+    ratio = electron_mass_c2/mass;
+
+    if(mass < 900.*MeV) nuclearStopping = false;
+
+    if (!EmModel(1)) SetEmModel(new G4BraggModel(),1);
+    EmModel(1)->SetLowEnergyLimit(MinKinEnergy());
+
+    // model limit defined for protons
+    eth = (EmModel(1)->HighEnergyLimit())*mass/proton_mass_c2;
+    EmModel(1)->SetHighEnergyLimit(eth);
+    AddEmModel(1, EmModel(1), new G4IonFluctuations());
+
+    if (!FluctModel()) SetFluctModel(new G4UniversalFluctuation());
+
+    if (!EmModel(2)) SetEmModel(new G4BetheBlochModel(),2);  
+    EmModel(2)->SetLowEnergyLimit(eth);
+    EmModel(2)->SetHighEnergyLimit(MaxKinEnergy());
+    AddEmModel(2, EmModel(2), FluctModel());  
+
+    isInitialised = true;
+  }
+  EmModel(1)->ActivateNuclearStopping(nuclearStopping);
+  EmModel(2)->ActivateNuclearStopping(nuclearStopping);
 }
 
@@ -175 +211 @@
 void G4hIonisation::PrintInfo()
 {
-  if(EmModel(1) && EmModel(2))
-    G4cout << "      Scaling relation is used from proton dE/dx and range."
-           << "\n      Delta cross sections and sampling from " 
-           << EmModel(2)->GetName() << " model for scaled energy > "
-           << eth/MeV << " MeV"
-           << "\n      Parametrisation from "
-           << EmModel(1)->GetName() << " for protons below."
-           << " NuclearStopping= " << nuclearStopping
+  if(EmModel(1) && EmModel(2)) {
+    G4cout << "      NuclearStopping= " << nuclearStopping
            << G4endl;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-void G4hIonisation::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
-                                         const G4DynamicParticle* dp,
-                                         G4double& eloss,
-                                         G4double& s)
-{
-  G4double ekin = dp->GetKineticEnergy();
-  if(nuclearStopping && ekin < ethnuc) {
-    G4double nloss = s*corr->NuclearDEDX(theParticle,couple->GetMaterial(),
-                                         ekin - eloss*0.5);
-    eloss += nloss;
-    //  G4cout << "G4ionIonisation::CorrectionsAlongStep: e= " << preKinEnergy
-    //     << " de= " << eloss << " NIEL= " << nloss << G4endl;
-    fParticleChange.ProposeNonIonizingEnergyDeposit(nloss);
   }
 }