Skip to content

Commit

Permalink
Merge pull request #646 from SimonRit/ffd_fixes
Browse files Browse the repository at this point in the history 
Fixed forced detection fixes
  • Loading branch information
@tbaudier
tbaudier authored Feb 5, 2024
2 parents f3232ef + ddc1885 commit 22c8dac
Show file tree
Hide file tree
Showing 2 changed files with 25 additions and 22 deletions.
36 changes: 18 additions & 18 deletions source/digits_hits/include/GateFixedForcedDetectionFunctors.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -543,7 +543,7 @@ namespace GateFixedForcedDetectionFunctor
const double & itkNotUsed(rayCastValue),
const VectorType & stepInMM,
const VectorType & itkNotUsed(source),
const VectorType & sourceToPixel,
const VectorType & pixelToSource,
const VectorType & nearestPoint,
const VectorType & farthestPoint)
{
Expand All @@ -559,7 +559,7 @@ namespace GateFixedForcedDetectionFunctor
/* The last material is the world material. One must fill the weight with
the length from source to nearest point and farthest point to pixel
point. */
VectorType worldVector = sourceToPixel + nearestPoint - farthestPoint;
VectorType worldVector = farthestPoint - nearestPoint - pixelToSource;
for (int i = 0; i < 3; i++)
{
worldVector[i] *= m_VolumeSpacing[i];
Expand Down Expand Up @@ -651,14 +651,14 @@ namespace GateFixedForcedDetectionFunctor
const double &itkNotUsed(rayCastValue),
const VectorType &stepInMM,
const VectorType &itkNotUsed(source),
const VectorType &sourceToPixel,
const VectorType &pixelToSource,
const VectorType &nearestPoint,
const VectorType &farthestPoint)
{
/* Compute ray length in world material
This is used to compute the length in world as well as the direction
of the ray in mm. */
VectorType worldVector = sourceToPixel + nearestPoint - farthestPoint;
VectorType worldVector = farthestPoint - nearestPoint - pixelToSource;
for (int i = 0; i < 3; i++)
{
worldVector[i] *= m_VolumeSpacing[i];
Expand Down Expand Up @@ -710,7 +710,7 @@ namespace GateFixedForcedDetectionFunctor

/* Final computation */
// double weight = std::exp(-rayIntegral) * DCScompton * GetSolidAngle(sourceToPixel) * (*m_ResponseDetector)(energy);
double weight = std::exp(-rayIntegral) * DCScompton * GetSolidAngle(sourceToPixel);
double weight = std::exp(-rayIntegral) * DCScompton * GetSolidAngle(pixelToSource);
if (m_generatePhotons)
{
//Accumulate(threadId, output, weight, m_Energy);
Expand All @@ -730,7 +730,7 @@ namespace GateFixedForcedDetectionFunctor
for (int i = 0; i < 3; i++)
{
photonDirection[i] = worldVector[i] / worldVectorNorm;
photonPosition[i] = sourceToPixel[i] * m_VolumeSpacing[i];
photonPosition[i] = -pixelToSource[i] * m_VolumeSpacing[i];
}
SavePhotonsparameters(threadId, photonPosition, photonDirection, weight, energy);
}
Expand Down Expand Up @@ -844,12 +844,12 @@ namespace GateFixedForcedDetectionFunctor
const double & itkNotUsed(rayCastValue),
const VectorType & stepInMM,
const VectorType & itkNotUsed(source),
const VectorType & sourceToPixel,
const VectorType & pixelToSource,
const VectorType & nearestPoint,
const VectorType & farthestPoint)
{
/* Compute ray length in world material. This is used to compute the length in world as well as the direction of the ray in mm. */
VectorType worldVector = sourceToPixel + nearestPoint - farthestPoint;
VectorType worldVector = farthestPoint - nearestPoint - pixelToSource;
for (int i = 0; i < 3; i++)
{
worldVector[i] *= m_VolumeSpacing[i];
Expand Down Expand Up @@ -886,7 +886,7 @@ namespace GateFixedForcedDetectionFunctor
}

/* Final computation */
double weight = std::exp(-rayIntegral) * DCSrayleigh * GetSolidAngle(sourceToPixel);
double weight = std::exp(-rayIntegral) * DCSrayleigh * GetSolidAngle(pixelToSource);
if (m_generatePhotons)
{
VectorType photonDirection;
Expand All @@ -905,7 +905,7 @@ namespace GateFixedForcedDetectionFunctor
for (int i = 0; i < 3; i++)
{
photonDirection[i] = worldVector[i] / worldVectorNorm;
photonPosition[i] = sourceToPixel[i] * m_VolumeSpacing[i];
photonPosition[i] = -pixelToSource[i] * m_VolumeSpacing[i];
}
SavePhotonsparameters(threadId, photonPosition, photonDirection, weight, m_Energy);
}
Expand Down Expand Up @@ -989,14 +989,14 @@ namespace GateFixedForcedDetectionFunctor
const double &itkNotUsed(rayCastValue),
const VectorType & stepInMM,
const VectorType & itkNotUsed(source),
const VectorType & sourceToPixel,
const VectorType & pixelToSource,
const VectorType & nearestPoint,
const VectorType & farthestPoint)
{
/* Compute ray length in world material
This is used to compute the length in world as well as the direction
of the ray in mm. */
VectorType worldVector = sourceToPixel + nearestPoint - farthestPoint;
VectorType worldVector = farthestPoint - nearestPoint - pixelToSource;
for (int i = 0; i < 3; i++)
{
worldVector[i] *= m_VolumeSpacing[i];
Expand All @@ -1023,7 +1023,7 @@ namespace GateFixedForcedDetectionFunctor
}

/* Final computation */
double weight = m_Weight * std::exp(-rayIntegral)*GetSolidAngle(sourceToPixel)/(4*itk::Math::pi);
double weight = m_Weight * std::exp(-rayIntegral)*GetSolidAngle(pixelToSource)/(4*itk::Math::pi);
if (m_generatePhotons)
{
VectorType photonDirection;
Expand All @@ -1042,7 +1042,7 @@ namespace GateFixedForcedDetectionFunctor
for (int i = 0; i < 3; i++)
{
photonDirection[i] = worldVector[i] / worldVectorNorm;
photonPosition[i] = sourceToPixel[i] * m_VolumeSpacing[i];
photonPosition[i] = -pixelToSource[i] * m_VolumeSpacing[i];
}
SavePhotonsparameters(threadId, photonPosition, photonDirection, weight, m_Energy);
}
Expand Down Expand Up @@ -1095,14 +1095,14 @@ namespace GateFixedForcedDetectionFunctor
const double & itkNotUsed(rayCastValue),
const VectorType & stepInMM,
const VectorType & itkNotUsed(source),
const VectorType & sourceToPixel,
const VectorType & pixelToSource,
const VectorType & nearestPoint,
const VectorType & farthestPoint)
{
/* Compute ray length in world material
This is used to compute the length in world as well as the direction
of the ray in mm. */
VectorType worldVector = sourceToPixel + nearestPoint - farthestPoint;
VectorType worldVector = farthestPoint - nearestPoint - pixelToSource;
for (int i = 0; i < 3; i++)
{
worldVector[i] *= m_VolumeSpacing[i];
Expand Down Expand Up @@ -1136,7 +1136,7 @@ namespace GateFixedForcedDetectionFunctor
}

/* Final computation */
double weight = m_Weight * std::exp(-rayIntegral)*GetSolidAngle(sourceToPixel)/(4*itk::Math::pi);
double weight = m_Weight * std::exp(-rayIntegral)*GetSolidAngle(pixelToSource)/(4*itk::Math::pi);
if (m_generatePhotons)
{
VectorType photonDirection;
Expand All @@ -1155,7 +1155,7 @@ namespace GateFixedForcedDetectionFunctor
for (int i = 0; i < 3; i++)
{
photonDirection[i] = worldVector[i] / worldVectorNorm;
photonPosition[i] = sourceToPixel[i] * m_VolumeSpacing[i];
photonPosition[i] = -pixelToSource[i] * m_VolumeSpacing[i];
}
SavePhotonsparameters(threadId, photonPosition, photonDirection, weight, m_Energy);
}
Expand Down
11 changes: 7 additions & 4 deletions source/digits_hits/src/GateFixedForcedDetectionActor.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -176,12 +176,15 @@ void GateFixedForcedDetectionActor::GetEnergyList(std::vector<double> & energyLi
}
auto energyHistogram = mSource->GetEneDist()->GetUserDefinedEnergyHisto();
double weightSum = 0.;
double energy = 0;
for (unsigned int i = 0; i < energyHistogram.GetVectorLength(); i++)
double leftEdge = energyHistogram.Energy(0);
for (unsigned int i = 1; i < energyHistogram.GetVectorLength(); i++)
{
energy = energyHistogram.Energy(i);
// See https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/GettingStarted/generalParticleSource.html?highlight=gps
double rightEdge = energyHistogram.Energy(i);
double energy = 0.5*(rightEdge+leftEdge);
leftEdge = rightEdge;
energyList.push_back(energy);
energyWeightList.push_back(energyHistogram.Value(energy));
energyWeightList.push_back(energyHistogram.Value(rightEdge));
weightSum += energyWeightList.back();
/* noise is desactivated */
if (mNoisePrimary == 0)
Expand Down

0 comments on commit 22c8dac

Please sign in to comment.