Update:

- Error when xraylib is not available ; - Delta map computation/Fresnel PCI can be disabled now.
OpenGATE · Sep 19, 2016 · 6dedb3c · 6dedb3c
1 parent 5012193
commit 6dedb3c
Show file tree

Hide file tree

Showing 4 changed files with 27 additions and 25 deletions.
diff --git a/source/digits_hits/include/GateFixedForcedDetectionActor.hh b/source/digits_hits/include/GateFixedForcedDetectionActor.hh
@@ -381,7 +381,6 @@ protected:
   double mInteractionSquaredIntegralOverDetector;
   G4String mSourceType;
   bool mGeneratePhotons;
-  bool mActivateFresnelDiffraction;
 
   /* Account for primary fluence weighting */
   InputImageType::Pointer PrimaryFluenceWeighting(const InputImageType::Pointer input);

diff --git a/source/digits_hits/include/GateFixedForcedDetectionFunctors.hh b/source/digits_hits/include/GateFixedForcedDetectionFunctors.hh
@@ -318,6 +318,8 @@ namespace GateFixedForcedDetectionFunctor
               {
               delta += (1 - Refractive_Index_Re(mat->GetElementVector()->at(i)->GetSymbol(), energyList[energy]/(keV), 1.0)) * mat->GetFractionVector()[i];
               }
+#else
+            G4Exception( "GateFixedForcedDetectionFunctors::CreateMaterialDeltaMap", "CreateMaterialDeltaMap", FatalException, "Xraylib is not available\n");
 #endif
             delta *= Density;
             it.Set(delta);
@@ -437,11 +439,6 @@ namespace GateFixedForcedDetectionFunctor
         m_generatePhotons = boolean;
         }
 
-      void setActivateFresnelDiffraction(const bool & boolean)
-        {
-        m_ActivateDiffraction = boolean;
-        }
-
     protected:
       inline void Accumulate(const rtk::ThreadIdType threadId,
                              float & output,
@@ -498,7 +495,6 @@ namespace GateFixedForcedDetectionFunctor
       unsigned int m_EnergyResolvedSliceSize;
       std::vector<std::vector<newPhoton> > m_PhotonList;
       bool m_generatePhotons;
-      bool m_ActivateDiffraction;
       };
 
     /* Most of the computation for the primary is done in this functor. After a ray
@@ -526,7 +522,6 @@ namespace GateFixedForcedDetectionFunctor
                              const VectorType & farthestPoint)
         {
         double *p = m_MaterialMu->GetPixelContainer()->GetBufferPointer();
-        double *q = m_MaterialDelta->GetPixelContainer()->GetBufferPointer();
         /* Multiply interpolation weights by step norm in MM to convert voxel
          intersection length to MM. */
         const double stepInMMNorm = stepInMM.GetNorm();
@@ -572,21 +567,25 @@ namespace GateFixedForcedDetectionFunctor
             vcl_exp(-rayIntegral) * (*m_EnergyWeightList)[i],
             m_EnergyList[i]);
 
-            double rayIntegral = 0.;
-            for (unsigned int j = 0; j < m_InterpolationWeights[threadId].size(); j++)
+            if (m_MaterialDelta != nullptr)
               {
-              rayIntegral += m_InterpolationWeights[threadId][j] * *q++;
+              double *q = m_MaterialDelta->GetPixelContainer()->GetBufferPointer();
+              double rayIntegral = 0.;
+              for (unsigned int j = 0; j < m_InterpolationWeights[threadId].size(); j++)
+                {
+                rayIntegral += m_InterpolationWeights[threadId][j] * *q++;
+                }
+
+              /* Matter wave : for photon, the formula is lambda = hc/E
+               * https://fr.wikipedia.org/wiki/Hypoth%C3%A8se_de_De_Broglie
+               */
+              double wavelength = h_Planck/(eV*s) * c_light/(m/s) / (m_EnergyList[i]/eV);
+              wavelength = wavelength * (m/mm); // To use the same unit as Geant4
+              rayIntegral *= (-2*itk::Math::pi/wavelength);
+              this->AccumulatePhase(output,
+              rayIntegral * (*m_EnergyWeightList)[i],
+              m_EnergyList[i]);
               }
-
-            /* Matter wave : for photon, the formula is lambda = hc/E
-             * https://fr.wikipedia.org/wiki/Hypoth%C3%A8se_de_De_Broglie
-             */
-            double wavelength = h_Planck/(eV*s) * c_light/(m/s) / (m_EnergyList[i]/eV);
-            wavelength = wavelength * (m/mm); // To use the same unit as Geant4
-            rayIntegral *= (-2*itk::Math::pi/wavelength);
-            this->AccumulatePhase(output,
-            rayIntegral * (*m_EnergyWeightList)[i],
-            m_EnergyList[i]);
             }
           }
 

diff --git a/source/digits_hits/src/GateFixedForcedDetectionActor.cc b/source/digits_hits/src/GateFixedForcedDetectionActor.cc
@@ -51,8 +51,7 @@ GateFixedForcedDetectionActor::GateFixedForcedDetectionActor(G4String name, G4in
     mInputRTKGeometryFilename(""),
     mEnergyResolvedBinSize(0),
     mSourceType("plane"),
-    mGeneratePhotons(false),
-    mActivateFresnelDiffraction(true)
+    mGeneratePhotons(false)
   {
   GateDebugMessageInc("Actor",4,"GateFixedForcedDetectionActor() -- begin"<<G4endl);
   pActorMessenger = new GateFixedForcedDetectionActorMessenger(this);
@@ -392,8 +391,11 @@ void GateFixedForcedDetectionActor::PreparePrimaryProjector(GeometryType::Pointe
   mPrimaryProjector->GetProjectedValueAccumulation().CreateMaterialMuMap(mEMCalculator,
                                                                          energyList,
                                                                          gate_image_volume);
-  mPrimaryProjector->GetProjectedValueAccumulation().CreateMaterialDeltaMap(energyList,
-                                                                            gate_image_volume);
+  if (mMaterialDeltaFilename != "" || mFresnelFilename != "")
+    {
+    mPrimaryProjector->GetProjectedValueAccumulation().CreateMaterialDeltaMap(energyList,
+                                                                              gate_image_volume);
+    }
   mPrimaryProjector->GetProjectedValueAccumulation().Init(mPrimaryProjector->GetNumberOfThreads());
   mPrimaryProjector->GetProjectedValueAccumulation().SetNumberOfPrimaries(mNoisePrimary);
   mPrimaryProjector->GetProjectedValueAccumulation().SetResponseDetector(&mEnergyResponseDetector);

diff --git a/source/physics/include/G4XrayBoundaryProcess.hh b/source/physics/include/G4XrayBoundaryProcess.hh
@@ -136,6 +136,8 @@ G4double G4XrayBoundaryProcess::GetRindex(G4Material *Material, G4double Energy)
 #ifdef GATE_USE_XRAYLIB
     for (unsigned int i = 0; i < Material->GetElementVector()->size(); ++i)
         delta += (1 - Refractive_Index_Re(Material->GetElementVector()->at(i)->GetSymbol(), Energy/(keV), 1.0)) * Material->GetFractionVector()[i];
+#else
+    G4Exception( "G4XrayBoundaryProcess::GetRindex", "GetRindex", FatalException, "Xraylib is not available\n");
 #endif
 
     return 1 - delta * Density;