[Gate-users] Bugs in DoseActor

[Javier Ors Anía]

Hi, thanks for sending the patch, I already had seen it in the user
list. But I need to tell you that there still rests a pair of
bugs/problems in DoseActor:

-Problem #1: The first problem is with this line in GateDoseActor.cc:

    mDoseImage.SetScaleFactor(1e12/mDoseImage.GetVoxelVolume());

This line is there to convert from energy/density, which is that the
Actor is collecting during the simulation, to energy/mass in Gray units,
which is what we want. The problem is that this scaling is applied at a
final stage just before saving the result in this lines of the file
GateImageWithStatistics.cc:

  if (!mIsValuesMustBeScaled) mValueImage.Write(mFilename);
  else {
    GateImage::iterator po = mScaledValueImage.begin();
    GateImage::iterator pi = mValueImage.begin();
    GateImage::const_iterator pe = mValueImage.end(); 
    while (pi != pe) {
      *po = (*pi)*mScaleFactor;
      ++pi;
      ++po;
    }
    mScaledValueImage.Write(mFilename);
    SetScaleFactor(factor);
  }

  if (mIsSquaredImageEnabled) {	
    UpdateSquaredImage();
    mSquaredImage.Write(mSquaredFilename);
  }
  if (mIsUncertaintyImageEnabled) {
    if (!mIsSquaredImageEnabled) UpdateSquaredImage();
    UpdateUncertaintyImage(numberOfEvents);
    mUncertaintyImage.Write(mUncertaintyFilename);
  }  
}

The first problem, as you can see, is that this correction is only
applied to the Dose image, by creating mScaledValueImage and saving it
instead of mValueImage. BUT, if you have chosen to save also the Squared
Dose image in the output, the mentioned correction is NOT being applied
to that image, so you are going to get the Squared Dose values in
different units than the Dose ones. If you try to externally reconstruct
the Uncertainty image, you are going to get wrong results unless you
also apply the corresponding correction to squared values (or reverse
the correction on dose values).

This first problem is only affecting to someone who wants to reconstruct
the uncertainty image externally (for example because is sending
different jobs to a grid and need to join the results, as was my case).
But now comes...


-Problem #2. This is an even more important problem than previous one,
but related to, because it is affecting the simulation output itself. 

As you may have noticed, in the Radiotherapy/Novice1 example, the
Uncertainty output has some strange issues, like relative uncertainties
higher than 1 in some pixels. The relative uncertainties are calculated
in this lines that you just changed in you last patch, maybe trying to
solve this issue?...

void GateImageWithStatistic::UpdateUncertaintyImage(int numberOfEvents)
{
  
  GateImage::iterator po = mUncertaintyImage.begin();
  GateImage::iterator pi = mValueImage.begin();
  GateImage::iterator pii = mSquaredImage.begin();
  GateImage::const_iterator pe = mValueImage.end();
  int N = numberOfEvents;
  while (pi != pe) {
    double squared = (*pii);
    double mean = (*pi);	
...	
    if (mean != 0.0 && N != 1 && squared != 0.0){
      *po = sqrt( (1.0/(N-1))*(squared/N - pow(mean/N, 2)))/(mean/N);
    }
    else *po = 1;

...
	
    ++po;
    ++pi;
    ++pii;
  }
}

This strange things in the uncertainty results are not because of the
first problem, as we are taking the dose values from mValueImage, which
is never scaled, so both squared and "mean" (in my opinion should not be
called mean as long as it is the total dose value, not the total
dose/N), are in the same default internal GEANT4 units, and the relative
uncertainty should be ok. BUT...

The problem is that if you calculate the uncertainty with this internal
GEANT4 units (without the 1e12 scale factor), the order of magnitude for
the squared dose image in some voxels can be as low as e-37, just check
it in the Squared-Dose output. 
This is not a problem during the simulation, as internally GEANT4 works
with double precision. The problem is that when we pass this values to
the GateImageWithStatistics class, this values are converted to simple
precision (float), and e-37 is on the edge of the simple precision low
limit, so we are calculating the uncertainty image with dose squared
values that have lost precision, and this is the reason because you see
values higher than 1 and maybe other artifacts that should not happen in
the uncertainty image.

The way I've followed to solve both problems at the same time is, first,
in GateDoseActor.cc comment the line:

//    mDoseImage.SetScaleFactor(1e12/mDoseImage.GetVoxelVolume())

to avoid the scaling at the final stage, and then I've moved the scaling
to the function UserSteppingActionInVoxel, to the line where the dose is
calculated:

    dose = edep/density;       

by changing it to:

dose = edep/density*1e12/mDoseImage.GetVoxelVolume(); 

In this manner the dose values are passed to the GateImageWithStatistics
class in Gray units, the squared dose values are also stored in this
units (solving problem #1), and they never reach the low limit for
single precision numbers (solving problem #2). 

P.S: I'm also sending this mail to the list so that other affected users
can check it.