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.1ns;
// 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->AddConstProperty("WLSTIMECONSTANT",12*ns);
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.);
Sheet_mat->SetMaterialPropertiesTable(scintMPT);
// Set Birks Constant and create an instance of Birks Law
Sheet_mat->GetIonisation()->SetBirksConstant(0.09*mm/MeV); </code>