Changeset 1055 for trunk/source/processes/electromagnetic/lowenergy/src/G4LivermoreComptonModel.cc

Timestamp:

May 28, 2009, 4:26:57 PM (15 years ago)

Author:

garnier

Message:

maj sur la beta de geant 4.9.3

File:

• trunk/source/processes/electromagnetic/lowenergy/src/G4LivermoreComptonModel.cc (modified) (15 diffs)

trunk/source/processes/electromagnetic/lowenergy/src/G4LivermoreComptonModel.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4LivermoreComptonModel.cc,v 1.1 2008/10/30 14:17:46 sincerti Exp $
-// GEANT4 tag $Name: geant4-09-02 $
-//
+// $Id: G4LivermoreComptonModel.cc,v 1.6 2009/04/18 18:29:34 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+//
+//
+// Author: Sebastien Inserti
+//         30 October 2008
+//
+// History:
+// --------
+// 18 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
+//                  - apply internal high-energy limit only in constructor 
+//                  - do not apply low-energy limit (default is 0)
+//                  - remove GetMeanFreePath method and table
+//                  - added protection against numerical problem in energy sampling 
+//                  - use G4ElementSelector
 
 #include "G4LivermoreComptonModel.hh"
@@ -37 +49 @@
 
 G4LivermoreComptonModel::G4LivermoreComptonModel(const G4ParticleDefinition*,
-                                             const G4String& nam)
-:G4VEmModel(nam),isInitialised(false)
+                                                 const G4String& nam)
+  :G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),
+   scatterFunctionData(0),crossSectionHandler(0)
 {
-  lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ?
+  lowEnergyLimit = 250 * eV; 
   highEnergyLimit = 100 * GeV;
-  SetLowEnergyLimit(lowEnergyLimit);
+  //  SetLowEnergyLimit(lowEnergyLimit);
   SetHighEnergyLimit(highEnergyLimit);
 
-  verboseLevel= 0;
+  verboseLevel=0 ;
   // Verbosity scale:
   // 0 = nothing 
@@ -52 +65 @@
   // 3 = calculation of cross sections, file openings, sampling of atoms
   // 4 = entering in methods
-  
-  G4cout << "Livermore Compton model is constructed " << G4endl
-         << "Energy range: "
-         << lowEnergyLimit / keV << " keV - "
-         << highEnergyLimit / GeV << " GeV"
-         << G4endl;
- 
+
+  if(  verboseLevel>0 ) { 
+    G4cout << "Livermore Compton model is constructed " << G4endl
+           << "Energy range: "
+           << lowEnergyLimit / eV << " eV - "
+           << highEnergyLimit / GeV << " GeV"
+           << G4endl;
+  }
 }
 
@@ -65 +79 @@
 G4LivermoreComptonModel::~G4LivermoreComptonModel()
 {  
-  delete meanFreePathTable;
-  delete crossSectionHandler;
-  delete scatterFunctionData;
+  if (crossSectionHandler) delete crossSectionHandler;
+  if (scatterFunctionData) delete scatterFunctionData;
 }
 
@@ -73 +86 @@
 
 void G4LivermoreComptonModel::Initialise(const G4ParticleDefinition* particle,
-                                       const G4DataVector& cuts)
+                                         const G4DataVector& cuts)
 {
   if (verboseLevel > 3)
     G4cout << "Calling G4LivermoreComptonModel::Initialise()" << G4endl;
 
-  InitialiseElementSelectors(particle,cuts);
-
-  // Energy limits
-  
-  if (LowEnergyLimit() < lowEnergyLimit)
-    {
-      G4cout << "G4LivermoreComptonModel: low energy limit increased from " << 
-        LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
-      SetLowEnergyLimit(lowEnergyLimit);
-    }
-
-  if (HighEnergyLimit() > highEnergyLimit)
-    {
-      G4cout << "G4LivermoreComptonModel: high energy limit decreased from " << 
-        HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
-      SetHighEnergyLimit(highEnergyLimit);
-    }
-
+  if (crossSectionHandler)
+  {
+    crossSectionHandler->Clear();
+    delete crossSectionHandler;
+  }
+  
   // Reading of data files - all materials are read
   
@@ -113 +114 @@
   shellData.LoadData(file);
 
-  meanFreePathTable = 0;
-  meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
-  
   if (verboseLevel > 2) 
     G4cout << "Loaded cross section files for Livermore Compton model" << G4endl;
 
-  G4cout << "Livermore Compton model is initialized " << G4endl
-         << "Energy range: "
-         << LowEnergyLimit() / keV << " keV - "
-         << HighEnergyLimit() / GeV << " GeV"
-         << G4endl;
-
-  //
-  
+  InitialiseElementSelectors(particle,cuts);
+
+  if(  verboseLevel>0 ) { 
+    G4cout << "Livermore Compton model is initialized " << G4endl
+           << "Energy range: "
+           << LowEnergyLimit() / eV << " eV - "
+           << HighEnergyLimit() / GeV << " GeV"
+           << G4endl;
+  }
+  //  
   if(isInitialised) return;
-
-  if(pParticleChange)
-    fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
-  else
-    fParticleChange = new G4ParticleChangeForGamma();
-    
+  fParticleChange = GetParticleChangeForGamma();
   isInitialised = true;
 }
@@ -148 +143 @@
     G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreComptonModel" << G4endl;
 
-  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
+  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;
+    
+  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);  
   return cs;
 }
@@ -155 +152 @@
 
 void G4LivermoreComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
-                                              const G4MaterialCutsCouple* couple,
-                                              const G4DynamicParticle* aDynamicGamma,
-                                              G4double,
-                                              G4double)
+                                                const G4MaterialCutsCouple* couple,
+                                                const G4DynamicParticle* aDynamicGamma,
+                                                G4double, G4double)
 {
+
   // The scattered gamma energy is sampled according to Klein - Nishina formula.
   // then accepted or rejected depending on the Scattering Function multiplied
@@ -174 +171 @@
   // (Nucl Phys 20(1960),15).
 
-  if (verboseLevel > 3)
-    G4cout << "Calling SampleSecondaries() of G4LivermoreComptonModel" << G4endl;
-
   G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
-  
-  if (photonEnergy0 <= lowEnergyLimit)
-  {
+
+  if (verboseLevel > 3) {
+    G4cout << "G4LivermoreComptonModel::SampleSecondaries() E(MeV)= " 
+           << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName() 
+           << G4endl;
+  }
+  
+  // low-energy gamma is absorpted by this process
+  if (photonEnergy0 <= lowEnergyLimit) 
+    {
       fParticleChange->ProposeTrackStatus(fStopAndKill);
       fParticleChange->SetProposedKineticEnergy(0.);
       fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
-      // SI - IS THE FOLLOWING RETURN NECESSARY ?
       return ;
-  }
+    }
 
   G4double e0m = photonEnergy0 / electron_mass_c2 ;
@@ -192 +192 @@
 
   // Select randomly one element in the current material
-  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
+  //  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
+  const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
+  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
+  G4int Z = (G4int)elm->GetZ();
 
   G4double epsilon0 = 1. / (1. + 2. * e0m);
@@ -212 +215 @@
       if ( alpha1/(alpha1+alpha2) > G4UniformRand())
       {
-        epsilon = std::exp(-alpha1 * G4UniformRand());  // std::pow(epsilon0,G4UniformRand())
+        // std::pow(epsilon0,G4UniformRand())
+        epsilon = std::exp(-alpha1 * G4UniformRand());  
         epsilonSq = epsilon * epsilon;
       }
@@ -238 +242 @@
   // Doppler broadening -  Method based on:
   // Y. Namito, S. Ban and H. Hirayama, 
-  // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code" 
-  // NIM A 349, pp. 489-494, 1994
+  // "Implementation of the Doppler Broadening of a Compton-Scattered Photon 
+  // into the EGS4 Code", NIM A 349, pp. 489-494, 1994
   
   // Maximum number of sampling iterations
@@ -257 +261 @@
       eMax = photonEnergy0 - bindingE;
      
-      // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
+      // Randomly sample bound electron momentum 
+      // (memento: the data set is in Atomic Units)
       G4double pSample = profileData.RandomSelectMomentum(Z,shell);
       // Rescale from atomic units
@@ -299 +304 @@
 
   if (photonEnergy1 > 0.)
-  {
-    fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
-  }
+    {
+      fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
+    }
   else
-  {
-    fParticleChange->SetProposedKineticEnergy(0.) ;
-    fParticleChange->ProposeTrackStatus(fStopAndKill);
-  }
+    {
+      photonEnergy1 = 0.;
+      fParticleChange->SetProposedKineticEnergy(0.) ;
+      fParticleChange->ProposeTrackStatus(fStopAndKill);   
+    }
 
   // Kinematics of the scattered electron
   G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE;
+
+  // protection against negative final energy: no e- is created
+  if(eKineticEnergy < 0.0) {
+    fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0 - photonEnergy1);
+    return;
+  }
   G4double eTotalEnergy = eKineticEnergy + electron_mass_c2;
 
@@ -329 +341 @@
   eDirection.rotateUz(photonDirection0);
 
-// SI - The range test has been removed wrt original G4LowEnergyCompton class
+  // SI - The range test has been removed wrt original G4LowEnergyCompton class
 
   fParticleChange->ProposeLocalEnergyDeposit(bindingE);
   
-  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),eDirection,eKineticEnergy) ;
+  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),
+                                                 eDirection,eKineticEnergy) ;
   fvect->push_back(dp);
 }
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4double G4LivermoreComptonModel::GetMeanFreePath(const G4Track& track,
-                                             G4double, // previousStepSize
-                                             G4ForceCondition*)
-{
-  const G4DynamicParticle* photon = track.GetDynamicParticle();
-  G4double energy = photon->GetKineticEnergy();
-  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
-  size_t materialIndex = couple->GetIndex();
-
-  G4double meanFreePath;
-  if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
-  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
-  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
-  return meanFreePath;
-}
-