Changeset 1055 for trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedComptonModel.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/G4LivermorePolarizedComptonModel.cc (modified) (12 diffs)

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

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4LivermorePolarizedComptonModel.cc,v 1.1 2008/10/30 14:16:35 sincerti Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4LivermorePolarizedComptonModel.cc,v 1.6 2009/05/03 08:29:55 sincerti Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
+// History:
+// --------
+// 02 May 2009   S Incerti as V. Ivanchenko proposed in G4LivermoreComptonModel.cc
+//
+// 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 "G4LivermorePolarizedComptonModel.hh"
@@ -38 +48 @@
 G4LivermorePolarizedComptonModel::G4LivermorePolarizedComptonModel(const G4ParticleDefinition*,
                                              const G4String& nam)
-:G4VEmModel(nam),isInitialised(false)
-{
-  lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ?
+:G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),scatterFunctionData(0),crossSectionHandler(0)
+{
+  lowEnergyLimit = 250 * eV; 
   highEnergyLimit = 100 * GeV;
-  SetLowEnergyLimit(lowEnergyLimit);
+  //SetLowEnergyLimit(lowEnergyLimit);
   SetHighEnergyLimit(highEnergyLimit);
   
@@ -53 +63 @@
   // 4 = entering in methods
 
+  if( verboseLevel>0 ) { 
   G4cout << "Livermore Polarized Compton is constructed " << G4endl
          << "Energy range: "
-         << lowEnergyLimit / keV << " keV - "
+         << lowEnergyLimit / eV << " eV - "
          << highEnergyLimit / GeV << " GeV"
          << G4endl;
-
+  }
 }
 
@@ -65 +76 @@
 G4LivermorePolarizedComptonModel::~G4LivermorePolarizedComptonModel()
 {  
-  delete meanFreePathTable;
-  delete crossSectionHandler;
-  delete scatterFunctionData;
+  if (meanFreePathTable)   delete meanFreePathTable;
+  if (crossSectionHandler) delete crossSectionHandler;
+  if (scatterFunctionData) delete scatterFunctionData;
 }
 
@@ -78 +89 @@
     G4cout << "Calling G4LivermorePolarizedComptonModel::Initialise()" << G4endl;
 
-  InitialiseElementSelectors(particle,cuts);
-
-  // Energy limits
-  
-  if (LowEnergyLimit() < lowEnergyLimit)
-    {
-      G4cout << "G4LivermorePolarizedComptonModel: low energy limit increased from " << 
-        LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
-      SetLowEnergyLimit(lowEnergyLimit);
-    }
-
-  if (HighEnergyLimit() > highEnergyLimit)
-    {
-      G4cout << "G4LivermorePolarizedComptonModel: 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
@@ -116 +115 @@
   shellData.LoadData(file);
 
-  //
   if (verboseLevel > 2) 
     G4cout << "Loaded cross section files for Livermore Polarized Compton model" << G4endl;
 
-  G4cout << "Livermore Polarized Compton model is initialized " << G4endl
+  InitialiseElementSelectors(particle,cuts);
+
+  if(  verboseLevel>0 ) { 
+    G4cout << "Livermore Polarized Compton model is initialized " << G4endl
          << "Energy range: "
-         << LowEnergyLimit() / keV << " keV - "
+         << 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;
 }
@@ -149 +146 @@
     G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" << G4endl;
 
+  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;
+
   G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
   return cs;
@@ -203 +202 @@
       fParticleChange->SetProposedKineticEnergy(0.);
       fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0);
-      // SI - IS THE FOLLOWING RETURN NECESSARY ?
       return;
     }
@@ -210 +208 @@
 
   // Select randomly one element in the current material
-
-  G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
+  //G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
+  const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
+  const G4Element* elm = SelectRandomAtom(couple,particle,gammaEnergy0);
+  G4int Z = (G4int)elm->GetZ();
 
   // Sample the energy and the polarization of the scattered photon
@@ -429 +429 @@
   else
     {
+      gammaEnergy1 = 0.;
       fParticleChange->SetProposedKineticEnergy(0.) ;
       fParticleChange->ProposeTrackStatus(fStopAndKill);
@@ -439 +440 @@
   G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
 
+  // SI -protection against negative final energy: no e- is created
+  // like in G4LivermoreComptonModel.cc
+  if(ElecKineEnergy < 0.0) {
+    fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0 - gammaEnergy1);
+    return;
+  }
+ 
   // SI - Removed range test
   
@@ -664 +672 @@
 }
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4double G4LivermorePolarizedComptonModel::GetMeanFreePath(const G4Track& track,
-                                                      G4double,
-                                                      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;
-}
-
-
+