Scintillation by Particle not recognizing ELECTRONSCINTILLATIONYIELD from material properties table

krisogren · October 12, 2023, 3:31pm

I am trying to build a model of a plastic scintillator with scintillation by particle type (after poring over the application developers book), but I’m running into an issue where I have defined the ELECTRONSCINTILLATIONYIELD in my material properties table, yet somehow Geant4 is kicking an error that this property is not defined.

This is the error:

-------- 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 ***
G4WT0 > G4Track (0x164c592d0) - track ID = 4, parent ID = 1
G4WT0 >  Particle type : C12 - creator process : conv, creator model : model_AugerElectron
G4WT0 >  Kinetic energy : 0 eV  - Momentum direction : (-0.114823,0.779454,-0.615847)
G4WT0 >  Step length : 10.1354 nm  - total energy deposit : 81.697 eV
G4WT0 >  Pre-step point : (9.64211,-1.06817,5.49988) - Physical volume : CherenkovDet (BC422)
G4WT0 >  - defined by : not available
G4WT0 >  Post-step point : (9.64211,-1.06817,5.49987) - Physical volume : CherenkovDet (BC422)
G4WT0 >  - defined by : ionIoni - step status : 3
G4WT0 >  *** Note: Step information might not be properly updated.
G4WT0 >
-------- EEEE -------- G4Exception-END --------- EEEE -------

Here is my material definition:

	//BC-422 Plastic Scintillator---------------------------------------------------------------------
	matBC422 = new G4Material("BC422", density=1.032*g/cm3, ncomponents=2);
	matBC422->AddElement(C, natoms=10);
	matBC422->AddElement(H, natoms=11);
	
	G4MaterialPropertiesTable* BC422_MPT = new G4MaterialPropertiesTable();
	
	G4double BC422_rise_time = 0.35*ns; //Source: St. Gobain data sheet
	G4double BC422_decay_time = 1.6*ns;
	
	vector<G4double> BC422_rindex = {1.573,1.575,1.58,1.595,1.605,1.63,1.644}; //Source: Nakamura 2015
	
	vector<G4double> BC422_energy_r = {WavelengthToEnergy(700)*eV,WavelengthToEnergy(656)*eV,WavelengthToEnergy(590)*eV,WavelengthToEnergy(487)*eV,
		                                  WavelengthToEnergy(436)*eV,WavelengthToEnergy(365)*eV,WavelengthToEnergy(340)*eV};
	
	
	vector<G4double> BC422_abs_length = {8.0*cm,8.0*cm,8.0*cm,8.0*cm,8.0*cm,8.0*cm,8.0*cm}; //Source: St. Gobain
	
	vector<G4double> BC422_emission = {2.37,7.19,11.09,15.17,19.81,25.28,30.95,37.63,43.38,50.44,57.30,64.92,71.22,79.02,86.64,96.10,100.,98.89,87.0,66.59,46.73,27.61,9.42,2.92};
	
	vector<G4double> BC422_energy_em = {WavelengthToEnergy(458)*eV,WavelengthToEnergy(450)*eV,WavelengthToEnergy(445)*eV,WavelengthToEnergy(440)*eV,
				WavelengthToEnergy(435)*eV,WavelengthToEnergy(430)*eV,WavelengthToEnergy(425)*eV,WavelengthToEnergy(420)*eV,WavelengthToEnergy(415)*eV,
			  WavelengthToEnergy(410)*eV,WavelengthToEnergy(405)*eV,WavelengthToEnergy(400)*eV,WavelengthToEnergy(395)*eV,WavelengthToEnergy(390)*eV,
			  WavelengthToEnergy(385)*eV,WavelengthToEnergy(380)*eV,WavelengthToEnergy(376.5)*eV,WavelengthToEnergy(375)*eV,WavelengthToEnergy(370)*eV,
			  WavelengthToEnergy(365)*eV,WavelengthToEnergy(360)*eV,WavelengthToEnergy(355)*eV,WavelengthToEnergy(350)*eV,WavelengthToEnergy(347)*eV};
	
	//energy in MeV
	vector<G4double> BC422_LY_energy = {0.01*MeV,1.*MeV,5.*MeV,10.*MeV,20.*MeV,30.*MeV,40.*MeV,50.*MeV,60.*MeV,70.*MeV,80.*MeV,90.*MeV,100.*MeV,110.*MeV,120.*MeV,130.*MeV};
	
	//Light yield in number of photons, Source: St. Gobain BC422 relative light yield plot
	vector<G4double> BC422_electron_LY = {110.,11000.,55000.,110000.,220000.,330000.,440000.,550000.,660000.,770000.,880000.,990000.,1100000.,1210000.,1320000.,1430000.};
	vector<G4double> BC422_proton_LY = {90.,9041.,28630.,63288.,141644.,226027.,317945.,408356.,504795.,607260.,703699.,800137.,902603.,1005068.,1109041.,1211507.};
	vector<G4double> BC422_alpha_LY = {10.,1506.,10548.,24110.,51233.,81370.,120548.,159726.,204932.,250137.,301370.,355616.,414384.,474658.,537945.,602740.};
	
	//apply properties to material table
	BC422_MPT->AddProperty("RINDEX",BC422_energy_r,BC422_rindex,false,true);
	BC422_MPT->AddProperty("ABSLENGTH",BC422_energy_r,BC422_abs_length,false,true);
	BC422_MPT->AddProperty("SCINTILLATIONCOMPONENT1",BC422_energy_em,BC422_emission,false,true);
	BC422_MPT->AddProperty("ELECTRONSCINTILLATIONYIELD",BC422_LY_energy,BC422_electron_LY,false,true);
	BC422_MPT->AddProperty("PROTONSCINTILLATIONYIELD",BC422_LY_energy,BC422_proton_LY,false,true);
	BC422_MPT->AddProperty("ALPHASCINTILLATIONYIELD",BC422_LY_energy,BC422_alpha_LY,false,true);
	BC422_MPT->AddConstProperty("SCINTILLATIONTIMECONSTANT1",BC422_decay_time,false);
	BC422_MPT->AddConstProperty("SCINTILLATIONRISETIME1",BC422_rise_time,false);
	//BC422_MPT->AddConstProperty("ELECTRONSCINTILLATIONYIELD1",1.0,false);
	BC422_MPT->AddConstProperty("RESOLUTIONSCALE",1.0,false);
	
	matBC422->SetMaterialPropertiesTable(BC422_MPT);
	BC422_MPT->DumpTable();

Interestingly, the DumpTable function seems to show that ELECTRONSCINTILLATIONYIELD is set, yet somehow it cannot be accessed when the particle is being tracked? DumpTable output below:

0: RINDEX
1.7712e-06   1.573
1.89e-06   1.575
2.10143e-06   1.58
2.54588e-06   1.595
2.84367e-06   1.605
3.39683e-06   1.63
3.64659e-06   1.644
9: GROUPVEL
1.7712e-06   186.925
1.8306e-06   186.809
1.99571e-06   184.527
2.32365e-06   179.982
2.69478e-06   177.35
3.12025e-06   170.516
3.64659e-06   162.814
16: ABSLENGTH
1.7712e-06   80
1.89e-06   80
2.10143e-06   80
2.54588e-06   80
2.84367e-06   80
3.39683e-06   80
3.64659e-06   80
17: PROTONSCINTILLATIONYIELD
0.01   90
1   9041
5   28630
10   63288
20   141644
30   226027
40   317945
50   408356
60   504795
70   607260
80   703699
90   800137
100   902603
110   1.00507e+06
120   1.10904e+06
130   1.21151e+06
20: ALPHASCINTILLATIONYIELD
0.01   10
1   1506
5   10548
10   24110
20   51233
30   81370
40   120548
50   159726
60   204932
70   250137
80   301370
90   355616
100   414384
110   474658
120   537945
130   602740
22: ELECTRONSCINTILLATIONYIELD
0.01   110
1   11000
5   55000
10   110000
20   220000
30   330000
40   440000
50   550000
60   660000
70   770000
80   880000
90   990000
100   1.1e+06
110   1.21e+06
120   1.32e+06
130   1.43e+06
23: SCINTILLATIONCOMPONENT1
2.70708e-06   2.37
2.7552e-06   7.19
2.78616e-06   11.09
2.81782e-06   15.17
2.85021e-06   19.81
2.88335e-06   25.28
2.91727e-06   30.95
2.952e-06   37.63
2.98757e-06   43.38
3.024e-06   50.44
3.06134e-06   57.3
3.0996e-06   64.92
3.13884e-06   71.22
3.17908e-06   79.02
3.22037e-06   86.64
3.26274e-06   96.1
3.29307e-06   100
3.30624e-06   98.89
3.35092e-06   87
3.39683e-06   66.59
3.444e-06   46.73
3.49251e-06   27.61
3.54241e-06   9.42
3.57303e-06   2.92
11: RESOLUTIONSCALE 1
30: SCINTILLATIONTIMECONSTANT1 1.6
33: SCINTILLATIONRISETIME1 0.35

I am able to run the extended optical examples LXe and OpNovice2 without problems, so I’m sure I’ve done something wrong, but I just can’t spot it. Also, if I turn off scintillation by particle and just specify a constant property for the scintillation yield, it works fine. Problem is I don’t have a Birks’ constant for the material, so scintillation by particle is the way to get the response to neutrons, etc.

As an aside, it was not clear to me from the application book that RESOLUTIONSCALE is a required parameter and does not just default to 1. I was previously getting errors of:

-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : mat202
      issued by : G4MaterialPropertiesTable::GetConstProperty()
Constant Material Property Index 11 not found.
*** Fatal Exception *** core dump ***

and it turned out that Index 11 is RESOLUTIONSCALE. Just a side note, but maybe that information will help someone.

Geant4 Version: 11.0.2
Operating System: MacOS 12.6.9
Compiler/Version: Apple clang 14.0.0
CMake Version: 3.24.1

dsawkey · October 12, 2023, 7:03pm

The error message is misleading. The particle is C12 so IONSCINTILLATIONYIELD needs to be set.

krisogren · October 12, 2023, 10:58pm

Thank you! Setting IONSCINTILLATIONYIELD solved my problem.

As an aside, does IONSCINTILLATIONYIELD apply to all heavy ions, or just C12? (Not clear in application book)

loydms · October 13, 2023, 1:58pm

I believe all ions other than protons and alpha particles

dsawkey · October 13, 2023, 2:02pm

I’m not sure. Protons, deuterons, tritons, alpha all get their own yields; nucleus and neutron use ION (see code below); and everything else uses ELECTRON.

  // Ions (particles derived from G4VIon and G4Ions) and recoil ions
  // below the production cut from neutrons after hElastic
  else if(pDef->GetParticleType() == "nucleus" ||
          pDef == G4Neutron::NeutronDefinition())
  {
    yieldVector = MPT->GetProperty(kIONSCINTILLATIONYIELD);

system · October 20, 2023, 2:03pm

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