Scintillation by Particle Type

loydms · March 22, 2022, 8:25pm

I am trying to simulate scintillation by particle type, and am running into an error.

I set ScintbyParticle type to True using
G4OpticalParameters::Instance()->SetScintByParticleType(true);

but I get an error when I try to set
pt20->AddConstProperty(“ALPHASCINTILLATIONYIELD”, 1255.23 / MeV);

and the application will not run due to this keyword not being in the property table. It works when I use ALPHASCINTILLATIONYIELD1, but then when I try to run I get an error:

G4Scintillation::PostStepDoIt(): Request for scintillation yield for energy deposit and particle

type without correct entry in MaterialPropertiesTable.

ScintillationByParticleType requires at minimum that

ELECTRONSCINTILLATIONYIELD is set by the user

I set alpha, ion, and electron scintillation yield using 1 at the end

John_McFee · March 23, 2022, 2:30am

Hi. The comand pt20->AddConstProperty(“ALPHASCINTILLATIONYIELD", 125523/MeV); has a couple of things wrong with it. You have to do the following:

You need to set ELECTRONSCINTILLATIONYIELD like this:
pt20->AddProperty(“ELECTRONSCINTILLATIONYIELD”, energy, scintYieldElectron);
note that you use AddProperty, not AddConstProperty. energy and scintYieldElectron are vectors of G4double which give yield (for an electron) energy pairs. So for example you might have for a particular scintillator,
std::vector energy = {1.0keV, 20.0MeV} (should span the particle range you are interested in)
std::vector scintYieldElectron = {20., 400000.} (replace with actual numbers for your scintillator)
next you have set the alpha yield
std::vector scintYieldAlpha = {1.25523., 25106.6}
pt20->AddProperty(“ALPHASCINTILLATIONYIELD”, energy, scintYieldAlpha);

(ALPHASCINTILLATIONYIELD1,2,3 are factors that ALPHASCINTILLATIONYIELD is multiplied by if you have multiple decay time constants with different relative yields. If you are dealing with only one time constant, you do not have to set them.)

loydms · March 23, 2022, 12:40pm

Thanks John, That fixed my issue. However, now I am getting an error that I have not set the neutron scintillation yield, as seen below. However, If I try to set NEUTRONSCINTILLATIONYIELD it says it is not a property in the table. Also, Neutrons are not charged and should not induced scintillation directly I believe.

-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : Scint01
issued by : G4Scintillation::PostStepDoIt

G4Scintillation::PostStepDoIt(): Request for scintillation yield for energy deposit and particle
type without correct entry in MaterialPropertiesTable.
ScintillationByParticleType requires at minimum that
ELECTRONSCINTILLATIONYIELD is set by the user

Missing MaterialPropertiesTable entry - No correct entry in MaterialPropertiesTable

*** Fatal Exception *** core dump ***
G4Track (0x7fe410e82730) - track ID = 1, parent ID = 0
Particle type : neutron - creator process : not available
Kinetic energy : 5 eV - Momentum direction : (0.00401353,-0.00373058,0.999985)
Step length : 500.008 um - total energy deposit : 0 eV
Pre-step point : (8.23816,-22.6081,-2.25) - Physical volume : Scint1 (MPND_gs20)

defined by : Transportation - step status : 1
Post-step point : (8.24017,-22.6099,-1.75) - Physical volume : Spacer (diff)
defined by : Transportation - step status : 1
*** Note: Step information might not be properly updated.

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

*** G4Exception: Aborting execution ***

loydms · March 23, 2022, 12:44pm

I have fixed this by also setting proton, deuteron, and triton yields

loydms · March 23, 2022, 2:06pm

John,

I have this working now by changing the code as you suggested. However, the results I am getting are not realistic. If I set all of the individual particles the same as if I just used a single yield, it seems as if the scintillation process is producing gammas instead of optical photons. Do I have to change how I have specified the emission spectrum, decay time, etc for scintillation by particle?

Also, I am getting a warning about the scintillation energy above the range for what I have specified, even though I used 6 MeV as my maximum range and am using thermal neutrons as a source (the produced alpha and triton should have a combined energy < 5 MeV.

John_McFee · March 23, 2022, 3:36pm

More information is needed. What makes you think that the scintillation process is producing gammas? Individual emission spectra do not depend on the type of ionizing particle, although their ratios can

I am not familiar with the warning you mention. What is the energy range of your emission spectra (i.e., SCINTILLATIONCOMPONENT1, etc.) and your yield vectors (i.e., ELECTRONSCINTILLATIONYIELD, etc.)?

loydms · March 23, 2022, 4:33pm

At each step I output particle, process and prestep/poststep PV and energy. When I don’t set SetScintByParticleType to true, I see that triton ends with a process of scintillaiton and the next particles stepped are optical photons (also, my alphas have the processes ionioni/noprocess but still produce optical photons). But if I change nothing but SetScintByParticleType(true), then the triton scintillation process is followed by gammas and no optical photons are produced except for Cherenkov.

Here is the code I use to set my scintillation properties. I called SetScintByParticleType in my main code after adding my physics list

  auto pt20 = new G4MaterialPropertiesTable;
double energies[2] = {1. * eV, 6. * MeV};
G4double alpha[2] = {1255.23 / MeV, 1255.23 / MeV};
G4double triton[2] = {1255.23 / MeV, 1255.23 / MeV};
G4double electron[2] = {3765. / MeV, 3765. / MeV};

pt20->AddProperty("ALPHASCINTILLATIONYIELD", energies, alpha, 2);
pt20->AddProperty("TRITONSCINTILLATIONYIELD", energies, triton, 2);
pt20->AddProperty("ELECTRONSCINTILLATIONYIELD", energies, electron, 2);
pt20->AddProperty("IONSCINTILLATIONYIELD", energies, alpha, 2);
pt20->AddProperty("PROTONSCINTILLATIONYIELD", energies, alpha, 2);
pt20->AddProperty("DEUTERONSCINTILLATIONYIELD", energies, alpha, 2);

// No special Fano factor
pt20->AddConstProperty("SCINTILLATIONYIELD", 1255.23 / MeV);
pt20->AddConstProperty("RESOLUTIONSCALE", 1.);
// Decay time
pt20->AddConstProperty("SCINTILLATIONTIMECONSTANT1", 57. * ns);
// All fast component, no slow component
pt20->AddConstProperty("SCINTILLATIONYIELD1", 1.);

G4double ScintEnergy[] = {
    1.771147508*eV, 1.785693715*eV, 1.806886491*eV, 1.835205916*eV, 1.859861711*eV, 1.887526235 *eV,
    1.911215667*eV, 1.937982147*eV, 1.968038263*eV, 1.993817314*eV, 2.025656659*eV, 2.052965257 *eV,
    2.086731237*eV, 2.11868072 *eV, 2.14858044 *eV, 2.176213997*eV, 2.214200693*eV, 2.243574408 *eV,
    2.270347301*eV, 2.294304242*eV, 2.322309111*eV, 2.354650359*eV, 2.376721902*eV, 2.406788116 *eV,
    2.433745705*eV, 2.469279938*eV, 2.493591512*eV, 2.522559239*eV, 2.556474115*eV, 2.586920261 *eV,
    2.609157128*eV, 2.627217637*eV, 2.650155755*eV, 2.673476749*eV, 2.692441907*eV, 2.711710765 *eV,
    2.726355401*eV, 2.746080968*eV, 2.761100218*eV, 2.77121984 *eV, 2.796776839*eV, 2.817562006 *eV,
    2.828136098*eV, 2.8547708  *eV, 2.871018238*eV, 2.892938281*eV, 2.898532973*eV, 2.915258618 *eV,
    2.926567493*eV, 2.937990045*eV, 2.955188453*eV, 2.966796809*eV, 2.978496723*eV, 2.990276024 *eV,
    2.996253925*eV, 3.008174378*eV, 3.020203583*eV, 3.038367271*eV, 3.050653453*eV, 3.069200332 *eV,
    3.081737652*eV, 3.100651336*eV, 3.139096176*eV, 3.191743043*eV, 3.225507311*eV, 3.239154774 *eV,
    3.245982849*eV, 3.259772991*eV, 3.273696694*eV, 3.2806075  *eV, 3.294726235*eV, 3.301758535 *eV,
    3.308837152*eV, 3.323134953*eV, 3.330272736*eV, 3.337474277*eV, 3.344707031*eV, 3.351971202 *eV,
    3.359216802*eV, 3.366443383*eV, 3.366544205*eV, 3.381295329*eV, 3.396090789*eV, 3.396227597 *eV,
    3.411206068*eV, 3.418692972*eV, 3.426230217*eV, 3.441598022*eV, 3.449219077*eV, 3.464740871 *eV,
    3.480313198*eV, 3.495989894*eV, 3.496098625*eV, 3.520093022*eV, 3.528102993*eV, 3.544381799 *eV,
    3.560849122*eV, 3.569083637*eV, 3.60277172 *eV, 3.654582691*eV, 3.735258252*eV, 3.829204446 *eV,
    3.927975404*eV, 4.010757254*eV, 4.130449603*eV, 4.245659698*eV, 4.367453005*eV, 4.470036937 *eV,
    4.577524753*eV, 4.763720755*eV, 4.886025126*eV, 5.06482313 *eV, 5.293804937*eV, 5.464820789 *eV,
    5.64730209 *eV, 5.864795368*eV, 6.149039422*eV
};

const G4int nPoints = sizeof(ScintEnergy)/sizeof(G4double);

G4double ScintFast[nPoints] = {
    0.004514673, 0.006772009, 0.006772009, 0.009029345, 0.006772009, 0.004514673,
    0.002257336, 0.004514673, 0.002257336, 0.004514673, 0.006772009, 0.004514673,
    0.004514673, 0.006772009, 0.006772009, 0.004514673, 0.006772009, 0.009029345,
    0.011286682, 0.013544018, 0.015801354, 0.020316027, 0.0248307  , 0.027088036,
    0.033860045, 0.036117381, 0.047404063, 0.058690745, 0.065462754, 0.074492099,
    0.090293454, 0.101580135, 0.11738149 , 0.128668172, 0.139954853, 0.158013544,
    0.173814898, 0.18510158 , 0.200902935, 0.214446953, 0.23476298 , 0.250564334,
    0.270880361, 0.293453725, 0.311512415, 0.329571106, 0.34537246 , 0.358916479,
    0.376975169, 0.399548533, 0.415349887, 0.431151242, 0.446952596, 0.460496614,
    0.476297968, 0.489841986, 0.505643341, 0.519187359, 0.53724605 , 0.553047404,
    0.571106095, 0.584650113, 0.598194131, 0.598194131, 0.591422122, 0.58013544,
    0.568848758, 0.553047404, 0.539503386, 0.519187359, 0.507900677, 0.492099323,
    0.478555305, 0.458239278, 0.440180587, 0.426636569, 0.413092551, 0.399548533,
    0.379232506, 0.35214447 , 0.365688488, 0.338600451, 0.300225734, 0.318284424,
    0.286681716, 0.264108352, 0.243792325, 0.227990971, 0.205417607, 0.182844244,
    0.148984199, 0.110609481, 0.124153499, 0.09255079 , 0.074492099, 0.056433409,
    0.042889391, 0.029345372, 0.018058691, 0.009029345, 0.006772009, 0.006772009,
    0.004514673, 0.004514673, 0.004514673, 0.006772009, 0.006772009, 0.006772009,
    0.004514673, 0.004514673, 0.004514673, 0.004514673, 0.006772009, 0.006772009,
    0.009029345, 0.006772009, 0.006772009,
};

G4double r_ind[nPoints] = {
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55, 1.55, 1.55,
    1.55, 1.55, 1.55, 1.55, 1.55
};

G4double abs[nPoints] = {
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm, 100*cm, 100*cm, 100*cm,
    100*cm, 100*cm, 100*cm
  };
pt20->AddProperty("SCINTILLATIONCOMPONENT1", ScintEnergy, ScintFast, nPoints);
pt20->AddProperty("RINDEX", ScintEnergy, r_ind, nPoints);
pt20->AddProperty("ABSLENGTH", ScintEnergy, abs, nPoints);

gs20->SetMaterialPropertiesTable(pt20);
m_materials[Material::GS20] = gs20;

John_McFee · March 23, 2022, 5:23pm

Hi.
1.The line pt20->AddConstProperty(“SCINTILLATIONYIELD”, 1255.23 / MeV); is unnecessary. The parameter SCINTILLATIONYIELD is only used for particle independent scintillation.

2.The values for the elements of alpha, triton, etc. do not make sense. Each entry should be a number of photons for the corresponding energy of the vector ‘energies’, not an energy per MeV. For example, for energies[j]=6*MeV, electron[j]=7531.38, etc. (This implies that if your scintillator has a constant yield per MeV, the values of the yield entries should increase with energy.)

3.The first entry of vector ‘energies’, 1.*eV, is unrealistically small. The limits should be the range for which you are examining scintillation. It likely will not cause an error, but it dosen’t make sense. Usually entries are in the keV, MeV range.

loydms · March 23, 2022, 6:04pm

Ah I see. I thought the scintillation yield was set as photons/MeV at each energy (as is when set by scintillationyield), but you are saying that when set by particle type is it set as an array of photons at each energy.

I was planning on having the scintillation yield be linear with energy, so I used 1 ev and 6 MeV just to have it interpolate flat. I can do that by changing the yield for each energy.

Thanks!

system · March 30, 2022, 6:05pm

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