Issues with Optical Properties in G4 v11.1.0

Hi there,

I am simulating the emission of scintillator light from protons in a segmented calorimeter made of polystyrene scintillator. The simulation is relatively simple, starting from example B1. I’ve been able to construct and run the simulation successfully with Geant4 version 11.0.3 on Intel MacOS 13.2, but when upgrading to G4 v11.1, I now encounter some issues with scintillator material properties.

Specifically, 11.1 seems to now require that energy vectors for refractive index, absorption spectrum and wavelength shifting spectrum, etc. be in increasing order, causing the error:

-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : mat215
issued by : G4MaterialPropertiesTable::AddProperty()
Energies in material property table must be in increasing order. Key: RINDEX Energy: 2.80507e-06
*** Fatal Exception *** core dump ***
**** Track information is not available at this moment
**** Step information is not available at this moment

-------- EEEE -------- G4Exception-END --------- EEEE -------

This made me realise that the files I use to import material data for refractive index, absorption length and wavelength shifting are in increasing wavelength order, which of course will give energy in decreasing order. Upon flipping the order of these vectors to remove this error but only running on G4 11.0.3, my simulation results change (blue curve in the below plots):

120MeV_100000.pdf (126.6 KB)

Regrettably, these data files have been passed down over the years from many different people working on simulation of scintillators, so I can’t say exactly where these definitions come from. I’m limited to 2 links on this post, so I cannot attach these files (for whatever use that might be) but I have added an excerpt from my code showing the material property table for the scintillator.

My question is why the results would be different if the order of these vectors is changed? If G4 wants the vectors in increasing energy order, how was version 11.0.3 interpreting this data to give “correct” behaviour?

This might also be a good opportunity to seek advice from the community on what’s actually necessary for good simulation of polystyrene scintillation, looking at the total light output for over millions of protons between energies 70 – 250 MeV? Our application is dosimetry in proton therapy.

It’s worth mentioning that a student working in our group with this simulation, @mayamaci, made a post on the forum about difficulties running this simulation on G4 11.0.3. but with Apple Silicon, again giving problems with scintillation. I am also investigating this, but first looking to fix issues with 11.1 on Intel MacOS described here. Perhaps information I’ve provided here can assist.

Many thanks for your help – please let me know if you need any further information and I’ll do my best to provide. I am somewhat new to this so apologies if I’ve missed something obvious!

Here is an excerpt from my detector construction showing the definition of the scintillator:

// Other definitions
G4double lightYield = 1;
G4double scintRes = 1;
G4double scintFastconst = 0.9ns;
G4double scintSlowconst = 2.1

// Define the energy spectrum of the scintillation photons
std::vector<G4double> photonEnergy = {
    2.034*eV, 2.068*eV, 2.103*eV, 2.139*eV,
    2.177*eV, 2.216*eV, 2.256*eV, 2.298*eV,
    2.341*eV, 2.386*eV, 2.433*eV, 2.481*eV,
    2.532*eV, 2.585*eV, 2.640*eV, 2.697*eV,
    2.757*eV, 2.820*eV, 2.885*eV, 2.954*eV,
    3.026*eV, 3.102*eV, 3.181*eV, 3.265*eV,
    3.353*eV, 3.446*eV, 3.545*eV, 3.649*eV,
    3.760*eV, 3.877*eV, 4.002*eV, 4.136*eV

std::vector<G4double> scintFast = {
    1., 1., 1., 1., 1.,
    1., 1., 1., 1., 1.,
    1., 1., 1., 1., 1.,
    1., 1., 1., 1., 1.,
    1., 1., 1., 1., 1.,
    1., 1., 1., 1., 1., 1., 1.

std::vector<G4double> scintSlow = {
    0.01, 1.00, 2.00, 3.00, 4.00,
    5.00, 6.00, 7.00, 8.00, 9.00,
    8.00, 7.00, 6.00, 4.00, 3.00,
    2.00, 1.00, 0.01, 1.00, 2.00,
    3.00, 4.00, 5.00, 6.00, 7.00,
    8.00, 9.00, 8.00, 7.00, 6.00, 5.00, 4.00

std::vector<G4double> reflectivity = {
    1., 1., 1., 1., 1.,
    1., 1., 1., 1., 1.,
    1., 1., 1., 1., 1.,
    1., 1., 1., 1., 1.,
    1., 1., 1., 1., 1.,
    1., 1., 1., 1., 1., 1., 1.

G4MaterialPropertiesTable* scintMPT = new G4MaterialPropertiesTable();

scintMPT->AddProperty("SCINTILLATIONCOMPONENT1", photonEnergy, scintFast);
scintMPT->AddProperty("SCINTILLATIONCOMPONENT2", photonEnergy, scintSlow);
scintMPT->AddProperty("RINDEX", ref_index_Energy, ref_index_value);
scintMPT->AddProperty("ABSLENGTH", absorbEnergy, Absorb);  //the bulk absorption spectrum
scintMPT->AddProperty("WLSCOMPONENT", wlsEnergy, wlsEmit);
scintMPT->AddProperty("REFLECTIVITY", photonEnergy, reflectivity); // doesn't seem to change results

scintMPT->AddConstProperty("SCINTILLATIONYIELD", lightYield/MeV);
scintMPT->AddConstProperty("RESOLUTIONSCALE", scintRes);
scintMPT->AddConstProperty("SCINTILLATIONTIMECONSTANT1", scintFastconst);
scintMPT->AddConstProperty("SCINTILLATIONTIMECONSTANT2", scintSlowconst);
scintMPT->AddConstProperty("SCINTILLATIONYIELD1", 1.);


// Set Birks Constant and create an instance of Birks Law
Sheet_mat->GetIonisation()->SetBirksConstant(0.09*mm/MeV); </code>


A Geant4 material property vector is actually a G4PhysicsVector, and G4PhysicsVector has always required that data are in order of increasing energy. Unfortunately this requirement never made it to the material property vector documentation. G4PhysicsVector interpolates (or fits a spline) to the discrete data, to find the property for all values of energy (in the range).

The history is that a material property table in order of decreasing energy worked well enough that noone noticed the error. It may have been that with closely spaced energy values, the interpolation was very close to the correct result. A recent change to G4PhysicsVector caused errors with MPT with decreasing energy, so we added the error message you are getting.

I don’t know what is causing the difference in results you observe. How to you light vs depth curve? Maybe with track->GetVertexPosition() ? How do the results depend on RINDEX, WLS, absorption? I guess, not strongly, so what happens if you remove all of that from the simulation? Just look at scintillation. Also, you’re sure you’re analyzing only scintillation photons, not Cerenkov?

I looked at the other post you mention, and as I wrote, there’s not enough information to comment. I would suggest starting with one of the optical examples (e.g. OpNovice) rather than example B1. Example B1 doesn’t have optical physics and you may be modifying it incorrectly.

Hi @dsawkey,

Thank you for your response – that is interesting to know.

I can confirm that I am just collecting scintillation photons. Here is how I am collecting dose and scintillation photons in the scintillator sensitive detector:

 // Bragg curve
G4double zmax = 128*3;
G4StepPoint* prePoint  = step->GetPreStepPoint();
G4StepPoint* postPoint = step->GetPostStepPoint();
G4double z1 = prePoint->GetPosition().z() + zmax*0.5;
G4double z2 = postPoint->GetPosition().z() + zmax*0.5;
if (z1 >= 0.0 && z2 <= zmax) {
    G4double z = z1 + G4UniformRand()*(z2-z1);
    if (step->GetTrack()->GetDefinition()->GetPDGCharge() == 0.) z = z2;
    G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
    analysisManager->FillH1(384, z, edep);
    // Quenched Bragg Curve
    G4Track* track = step->GetTrack();
    if (track->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) {  // particle is optical photon
        if (track->GetParentID() > 0) {  // particle is secondary
            if (track->GetCreatorProcess()->GetProcessName() == "Scintillation") {
                analysisManager->FillH1(385, z, 1);

It appears that the scintillation parameters and energy-dependent refractive index and absorption length data are necessary to give good results. I’ve attached two PDFs showing how the data changes with addition of each component, with both increasing energy and decreasing energy vectors, on G4 11.0.3. It seems that with increasing energy, the refractive index seems to give most of the “incorrect” behaviour. In particular, we seem to get divergences in the photon count corresponding to the thicknesses of the 3 mm square scintillator sheets that make up the segmented calorimeter.

Decreasing Energy.pdf (451.9 KB)
Increasing Energy.pdf (348.9 KB)

I’ve since emailed the manufacturer of the scintillator to see if they can provide a new set of data for these parameters, so that I can confirm their validity.

At this point, not sure what else I can try! I’ll at least do some more testing with G4 11.0.3 on Apple Silicon to see if I can get that working.

I’m not sure I understand the logic. I think you want to score the z coordinate of where the scintillation photons are created. However, the code scores the z position of the post-step point of every optical photon step.

A partial understanding of the plots: with only scintillation, the photons don’t interact within a volume. Most steps are transportation limited, ending on the boundary of a scintillator sheet. Optical photons don’t transmit into volumes with no RINDEX, so the plot doesn’t extend much deeper than the Bragg curve.

Once interactions are added, the volume boundary steps are washed out. With RINDEX defined, optical photons propagate deeper than the Bragg curve.

In order to score the origin position of the optical photons, you could iterate over secondaries as done in the OpNovice2 example (and others).


Thank you for your insight into how the properties were changing the graphs!

Yes, I am interested in scoring the location of where each scintillation photon is created or rather, the number of photons produced as a function of depth. At this point, I’m interested in collecting all scintillation photons, regardless of the process. So if I understand the OpNovice2 example, you suggest to tabulate this data like this:

static G4ParticleDefinition* opticalphoton = G4OpticalPhoton::OpticalPhotonDefinition();

G4AnalysisManager* analysisMan = G4AnalysisManager::Instance();
Run* run = static_cast<Run*>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());

G4Track* track         = step->GetTrack();
G4StepPoint* endPoint  = step->GetPostStepPoint();
G4StepPoint* startPoint = step->GetPreStepPoint();
G4double z1 = prePoint->GetPosition().z() + zmax*0.5;
G4double z2 = postPoint->GetPosition().z() + zmax*0.5;
G4double z = z1 + G4UniformRand()*(z2-z1);

const G4DynamicParticle* theParticle = track->GetDynamicParticle();
const G4ParticleDefinition* particleDef = theParticle->GetParticleDefinition();

TrackInformation* trackInfo = (TrackInformation*) (track->GetUserInformation());

if(particleDef == opticalphoton)
  analysisMan->FillH1(385, z, 1);

The math that I am doing with the pre- and post- step is from the way z is calculated in example TestEm7 when tabulating dose as a function of depth, so I am not recording the post-step per se. Maybe there is a better way to do this? e.g. with track->GetVertexPosition() as you mentioned above?

Is it necessary to loop over the secondaries, if I am only interested in the creation process (i.e. scintillation)? Do I need to check that these optical photons are scintillation photons? Also, I’m not sure how this way of doing things is different to how I originally did this: are we not looking at each optical photon step in the same way?

Thank you again for your help, it’s really appreciated!

My main point is: the optical photon may take many steps in the sensitive detector. Every time the photon interacts, or reaches a geometry boundary, will be a new step. Energy deposit is inherently a step-by-step mechanism, so that part of the code is fine. However, any given scintillation photon is only created once. The code needs to be written such that any photon is only counted once. Maybe I’m simply misunderstanding (only looking at a code snippet), but it appears that FillH1 will be called for every optical photon step, and thus counted multiple times (with different z each time).

There are several things you could do:

  • move the counting to UserTrackingAction or UserStackingAction, where it will only happen once per track
  • kill the photon once it is scored (assuming you don’t want to follow it)
  • add a UserTrackInformation to indicate whether a photon has been scored
  • simplest would be to add if (step->IsFirstStepInTrack()) to the conditions (but this doesn’t exist)
  • loop over secondaries created in each step. Each secondary is created once, so there will be no double counting.

For looping over secondaries, I meant the OpNovice2 code here:
where energy and time are scored. You’d want to add z coordinate.

I think z is set to the post-step coordinate for optical photons, giving you those spikes in the plots at the volume boundaries.

Hi @dsawkey,

We have some luck! It appears that it was the photon collection method that was causing the problems. Implementing the photon collection using what you suggested from OpNovice2 gives better results that I can immediately see are more physical, particularly in the photon tail after the Bragg peak. As you said, that condition setting z = z2 for photons (that I was overlooking) was causing the issues with the divergences at the boundaries of scintillator sheets.

Now it seems that, as you predicted, the refractive index and bulk absorption data actually has little impact on the final results. I am able to run the simulation in version 11.1 with no problems, regardless if I choose to import this data (with increasing energy).

One minor detail that I noticed is that I now see this message, only about once in a simulation of 10^5 or so protons:

-------- WWWW ------- G4Exception-START -------- WWWW -------
*** G4Exception : OpBoun06
issued by : G4OpBoundaryProcess
G4OpBoundaryProcess: Opticalphoton step length: 3.3133e-09 mm.
This is larger than the threshold 1e-09 mm to set status StepTooSmall.
Boundary scattering may be incorrect.

*** This is just a warning message. ***
-------- WWWW -------- G4Exception-END --------- WWWW -------

This occurs with or without refractive index/absorption length data. Not sure if this is something to be concerned about?

I suspect this change to the simulation will help with our issues running the simulation on Apple Silicon – I will test this and post what I find on the other post. For now, I’ll mark this topic as resolved – thank you for your help!

Great! Glad to be of help.

Re: the new warning. This can be ignored if it happens infrequently. More details: The G4OpBoundaryProcess uses the step length internally, as a flag to set material properties when a photon is reflected back into the volume it is already in. The idea is that if the step length is zero, the photon is reflected. However, sometimes “zero” isn’t precisely 0.0, and I don’t know what the tolerance should be. This particular photon may be sent into the wrong volume.

1 Like

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.