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

Timestamp:

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

Author:

garnier

Message:

update processes

File:

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

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4BraggModel.cc,v 1.16 2007/07/28 13:30:53 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4BraggModel.cc,v 1.21 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -50 +50 @@
 // 15-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
 // 25-04-06 Add stopping data from PSTAR (V.Ivanchenko)
+// 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio, 
+//          CorrectionsAlongStep needed for ions(V.Ivanchenko)
 
 // Class Description:
@@ -67 +69 @@
 #include "G4Electron.hh"
 #include "G4ParticleChangeForLoss.hh"
+#include "G4LossTableManager.hh"
+#include "G4EmCorrections.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -74 +78 @@
 G4BraggModel::G4BraggModel(const G4ParticleDefinition* p, const G4String& nam)
   : G4VEmModel(nam),
-  particle(0),
-  protonMassAMU(1.007276),
-  iMolecula(0),
-  isIon(false)
+    particle(0),
+    protonMassAMU(1.007276),
+    iMolecula(0),
+    isIon(false),
+    isInitialised(false)
 {
   if(p) SetParticle(p);
+  SetHighEnergyLimit(2.0*MeV);
+
   lowestKinEnergy  = 1.0*keV;
   theZieglerFactor = eV*cm2*1.0e-15;
@@ -104 +111 @@
 {
   if(p != particle) SetParticle(p);
-  G4String pname = particle->GetParticleName();
-  if(particle->GetParticleType() == "nucleus" && 
-     pname != "deuteron" && pname != "triton") isIon = true;
-
-  if(pParticleChange)
-    fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>
-                                                              (pParticleChange);
-  else
-    fParticleChange = new G4ParticleChangeForLoss();
+
+  // always false before the run
+  SetDeexcitationFlag(false);
+
+  if(!isInitialised) {
+    isInitialised = true;
+
+    G4String pname = particle->GetParticleName();
+    if(particle->GetParticleType() == "nucleus" && 
+       pname != "deuteron" && pname != "triton") isIon = true;
+
+    corr = G4LossTableManager::Instance()->EmCorrections();
+
+    if(pParticleChange) {
+      fParticleChange = 
+        reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
+    } else {
+      fParticleChange = new G4ParticleChangeForLoss();
+    }
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4BraggModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
+                                            const G4Material* mat,
+                                            G4double kineticEnergy)
+{
+  // this method is called only for ions
+  G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
+  GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
+  return q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4BraggModel::GetParticleCharge(const G4ParticleDefinition* p,
+                                         const G4Material* mat,
+                                         G4double kineticEnergy)
+{
+  // this method is called only for ions
+  return corr->GetParticleCharge(p,mat,kineticEnergy);
 }
 
@@ -123 +163 @@
                                                  G4double maxKinEnergy)
 {
-
   G4double cross     = 0.0;
   G4double tmax      = MaxSecondaryEnergy(p, kineticEnergy);
-  G4double maxEnergy = min(tmax,maxKinEnergy);
+  G4double maxEnergy = std::min(tmax,maxKinEnergy);
   if(cutEnergy < tmax) {
 
@@ -203 +242 @@
 
   return dedx;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4BraggModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
+                                        const G4DynamicParticle* dp,
+                                        G4double& eloss,
+                                        G4double&,
+                                        G4double length)
+{
+  if(nuclearStopping) {
+
+    G4double preKinEnergy = dp->GetKineticEnergy();
+    G4double e = preKinEnergy - eloss*0.5;
+    if(e < 0.0) e = preKinEnergy*0.5;
+    G4double nloss = length*corr->NuclearDEDX(dp->GetDefinition(),
+                                              couple->GetMaterial(),
+                                              e,false);
+
+    // too big energy loss
+    if(eloss + nloss > preKinEnergy) {
+      nloss *= (preKinEnergy/(eloss + nloss));
+      eloss = preKinEnergy;
+    } else {
+      eloss += nloss;
+    }
+    /*
+    G4cout << "G4ionIonisation::CorrectionsAlongStep: e= " << preKinEnergy
+           << " de= " << eloss << " NIEL= " << nloss 
+           << " dynQ= " << dp->GetCharge()/eplus << G4endl;
+    */
+    fParticleChange->ProposeNonIonizingEnergyDeposit(nloss);
+  }
 }
 
@@ -268 +340 @@
 
   vdp->push_back(delta);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4BraggModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
+                                          G4double kinEnergy)
+{
+  if(pd != particle) SetParticle(pd);
+  G4double tau  = kinEnergy/mass;
+  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
+                  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
+  return tmax;
 }