Hi everyone

I wrote a simulation program to calculate the coincident time resolution of LYSO detector. There is one point source emitting 511 keV radiation isotropically and two LYSO detectors of size 2.7x2.7x18 mm. I defined a scorer as a cathode at the end of the scintillator and saved the global time when optical photons entered the scoring volume. I am attaching the code snippet below. I was expecting the FWHM of difference of time (photonTime-photonTime1) captured to be in range of 100s of picosecond but rather I was getting it to be around 4 ns, which is quite high. Am I doing it correctly or what is wrong here ? Kindly let me know.

```
if(stepstatus==fGeomBoundary){
// check if the photon is absorbed in the sensitive volume
if (currentPhysicalName == "Cathode"){
if (particleName == "opticalphoton"){
photonTime = aStep -> GetPreStepPoint() -> GetGlobalTime();
outFile.open("cathodelnew.data", std::ofstream::app);
outFile <<std::setprecision(15)<<photonTime/s << " "<<evt1<<std::endl;
outFile.close();
}
}
}
if(stepstatus==fGeomBoundary){
if (currentPhysicalName == "Cathode1"){
if (particleName == "opticalphoton"){
photonTime1 = aStep -> GetPreStepPoint() -> GetGlobalTime();
outFile.open("cathodel1new.data", std::ofstream::app);
outFile <<std::setprecision(15)<<photonTime1/s << " "<<evt1<<std::endl;
outFile.close();
}
}
}}
```

Hello @dsawkey

Thank you for your response.

- Here are my parameters:

`G4double ScintEnergy[nEntries] = {2.034 * eV, 4.136 * eV};

G4double ScintFast[nEntries] = {1.0,1.0};

```
G4double LYSORefractionIndex[nEntries] = {1.9,1.9};
G4double LYSOAbsorptionLength[nEntries] = {1.2*cm,1.2*cm};
G4MaterialPropertiesTable* LYSOMPT = new G4MaterialPropertiesTable();
LYSOMPT->AddProperty("RINDEX",ScintEnergy,LYSORefractionIndex,
nEntries);
LYSOMPT->AddProperty("ABSLENGTH",ScintEnergy, LYSOAbsorptionLength,
nEntries);
LYSOMPT-> AddProperty("SCINTILLATIONCOMPONENT1",ScintEnergy,ScintFast,nEntries);
LYSOMPT->AddConstProperty("SCINTILLATIONYIELD",33./keV);
LYSOMPT->AddConstProperty("RESOLUTIONSCALE",1.);
LYSOMPT->AddConstProperty("SCINTILLATIONTIMECONSTANT1",36.*ns);
LYSOMPT->AddConstProperty("SCINTILLATIONYIELD1",1.);`
```

- I am not sure about verbose tracking yet. But as shown in the code snippet, I find individual time and then subtract it. The individual time looks like this. Just including 10 values. The value after space is event number.

```
1.13011928033542e-07 314
4.79966003650184e-08 314
6.29807737831337e-08 314
5.45306872701844e-09 314
5.08149316159875e-08 314
2.46186757674192e-08 314
5.90541163313098e-09 314
2.79142089216495e-08 314
2.44534502117008e-08 314
4.34485277607281e-08 314
```

The time the optical photons are emitted is drawn from an exponential distribution. This line:

sets the mean value to be 36 ns. If you put your time values in a histogram, you should see that exponential distribution.

1 Like

Oh I see this. But then, what should I put here on this property? This is the actual value of the decay constant of LYSO. If thatâ€™s the case then how I can get picosecond scale of resolution ?

Kindly help

I also request @John_McFee , @weller , @bmorgan, @mkelsey and other experts in the field to provide some insights if possible. I have been stuck on this for a while and unable to figure it out. Your help will be highly appreciated.

You are measuring the arrival time at the photocathode of the optical photons from scintillation. That is dictated by the scintillator decay constant and the optical geometry. So it its distribution should be roughly an expoential with a time constant similar to the decay constant. The â€śrange of 100 psâ€ť that you expect is for the FWHM from a timing coincidence experiment. Methods vary, but they all involve triggering very close to the beginning of a pulse or at a fixed fraction of the pulse height (e.g., search â€śconstant fraction discriminatorâ€ť). As such, the time distribution is much narrower than the decay time distribution. I suppose you could simulate that in Geant4, but off line analysis is likely easier.

Thank you so much for the reply.

I actually read it in an article about calculating the average time when the photons reach the scorer to calculate intrinsic time resolution.

"

To obtain the ITR of each configuration, we obtained the ATvalue, this value represents the average time at which the optical photons reach to the Scorer and it is obtained event byevent

link for article: View of Intrinsic time resolution and efficiency study for simulated scintillators plastics with Geant4

"

So is doing from this approach not correct ? I calculated the average time of photons sorted by event number and then subtracted corresponding time for same event number for two detectors.

Iâ€™ll have to read the paper in detail, but at first glance it seems to me that the method you describe will yield a time difference with uncertainty of the order of the decay time constant. The response from each detector is exponentially distributed (P=1/lambda*exp(-lambda*t) where lambda = 1/35ns). The variance of the time of arrival of one detector is 1/lambda^2 and the variance of the time difference between the two detectors (assuming they are the same) is 2/lambda^2. So the standard deviation of the time difference between the two detectors should be ~1.4*36ns.

1 Like

I understood your point. Thanks for the detailed description, How does people finding CTR perform the simulation, this is still a little unclear to me. But thank you for your response.

There is a digitization package in Geant4 (See Digitization in the Users Guides) which can be used to simulate some of the aspects of a signal processing chain (ADCs. TDCs, readouts, trigger logic, etc.), but I am not familiar with it. Iâ€™ll let the experts comment on it further. You can analyze the time response outside Geant4, which is generally my preference.

If you save the arrival time of optical photons that arrive at the PMT front face and their corresponding event number in a file you have exactly what a PMT photocathode would â€śseeâ€ť event by event. One should be able to use fast prototyping software like Matlab or Octave to analyse that to simulate pretty much any signal processing technique ypu can think of (e.g., single photon detection, threshold detection, CFD, etc.). A number of those will have significantly better time resolution than one would get from the raw scintillation decay signal.

1 Like