Problem with radioactive source in optical simulation

Recording, Visualizing and Persisting Data
1

I made a scintillator detector and incorporated all the optical properties. When I am simulating with radioactive source like Co-60 I am getting photopeak at 1.33 MeV anf 1.55 MeV which is not correct energy deposition profile. But after excluding optical properties I am getting correct energy spectrum( photopeak at 1.17 MeV anf 1.33 MeV) with Co-60 source .

2

Not only with Co-60 source but also with mono energetic source I am getting higher energy deposition.

3

It sounds like you have some incorrect optical parameters, e.g., in scintillation yield or an error in converting the number of generated optical photons back into an equivalent MeVee.

1 Like
4

Hello @John_McFee
My physics list have
G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
RegisterPhysics(opticalPhysics);
opticalPhysics->SetScintillationYieldFactor(0.001);

Earlier SetScintillationYieldFactor was set 0.08 , Now I changed to 0.001 to get nearly ~1.17 and ~1.33 .
Now my question is, should we always vary this value to get the energy deposition accordingly?

5

What version of Geant4 are you running (as far as I know, SetScintillationYieldFactor is not used for Geant4 versions greater than 10.7)? The parameter is used in particle-dependent scintillation to set the ratio of scintillation yield for a particle relative to the that of electrons. For which particle are you doing SetScintillationYieldFactor?

6

@John_McFee I am using Geant4 version 10.7.3. I am doing doing the simulation with gamma. Since I am using Co-60 ,the most common gamma rays emitted by Co-60 are two primary energies: 1.17 MeV and 1.33 MeV. I can send you the full code, if you say.

7

I am familiar with Co60. I asked because your definitions are puzzling. For gamma rays, one usually sets the scintillation yield factor to 1.0 in the physics list definitions and in the material definitions part of the detector construction module define ELECTRONSCINTILLATIONYIELD to be the appropriate vector for the scintillator in terms of number of scintillation photons per MeV, e.g.,

G4double energies_electron[2] = {1.0*keV, 20.0*MeV};
G4double yield_electron[2]    = {49., 980000.};
LaCl3_mt->AddProperty("ELECTRONSCINTILLATIONYIELD", energies_electron, yield_electron, 2);

What scintillator material are you using and what values do you have for the quantities I am calling ‘energies_electron’ and ‘yield_electron’?

1 Like
9

You are clearly doing particle-independent scintillation (i.e., just gamma rays in this case). So get rid of the stuff you added (after 'Now I added"). You only need to define ScintillationYieldFactor if you are doing particle-dependent scintillation. Let optical physics use the default value (which is 1.0, I believe). Btw, 340000 should be 204000.

The other optical parameters seem reasonable. However, make the ‘ScinEnergy’ vector the same as the ‘energy’ vector. There is no good reason to make them different and there can be lots of problems sometimes if they are different (although likely not in this case).

Lastly, how are you evaluating the energy of the produced optical photons? That is, you count x scintillation photons somewhere after each incident gamma ray and convert x to an equivalent gamma (electron) energy. How and where do you do that?

1 Like
10

Yes! You are right @John_McFee . When I ignore energy deposition due to optical photon I am getting correct spectra. Wondering in real experiment why don’t we get photopeak at 1.33 MeV anf 1.55 MeV.?

12

The monoenergetic gamma ray source emits the gamma ray in an event at t=0 (in a sense instantaneous). So the time profiel will follow the decay constant of the scintillator (~20ns). The Co60 source has a half life of 5.27 years. That likely explains the difference.

1 Like
13

Ok. Thank you @John_McFee for the clarification. I want to get the arrival time of optical photon to the Cathode with the condition that the clock starts at zero when optical photon created.
Now I used ,
if (currentPhysicalName == “Cathode”){ **
** const G4String particleName
** = aStep->GetTrack()->GetDefinition()->GetParticleName();**
** if (particleName == “opticalphoton”){ **
** Time = aStep → GetPreStepPoint() → GetGlobalTime(); **
** } **
** }**

We know global time refers the time since the beginning of the event of the track after the completion of the current step. And GetDeltaTime() gives time of flight spent by a step.
But for my case what to do? Can @John_McFee you please tell ?

14

I want to get the arrival time of optical photon to the Cathode with the condition that the clock starts at zero when optical photon created.

Try using local time (time since a track is created), e.g., GetLocalTime(), on each optical photon track.

1 Like
15

Hello @John_McFee !
To get the time since the track of optical photon was created I did the following .
if (currentPhysicalName == “Cathode”){
const G4String particleName
= aStep->GetTrack()->GetDefinition()->GetParticleName();
if (particleName == “opticalphoton”)
{
G4double time = aStep->GetPostStepPoint()->GetLocalTime();

	} 
}

But I don’t understand whether I should put the condition of currentPhysicalName == “Cathode”
My next query is,
Now there are many optical photon in each event, so this time is an average value?

16

Hi,
The condition to select the “Cathode” ensures that you are evaluating the time taken by the optical photon from its creation until it reaches the photocathode. That is the time for the particular optical photon associated with the step. It is one optical photon among the many generated per scintillation. So there will be a time distribution of optical photons reaching the cathode. How to process those times is up to you. You could save them in a .root file or .csv file and post process either in Geant4 or externally. If you just want an average time per event, you can calculate the average using accumulables (seeSteppingAction in example B1).

1 Like
17

@John_McFee thanks again. I saved this G4double time = aStep->GetPostStepPoint()->GetLocalTime(); in a .csv file. Now I got value for each event.
I am confused that does this value giving the time for any one optical photon out of many in each event?
If so then how it decides which time to store out of many in each event?

I could see in RunAction of example B1 this line are attached,
void B1RunAction::BeginOfRunAction(const G4Run*)
{
// inform the runManager to save random number seed
G4RunManager::GetRunManager()->SetRandomNumberStore(false);
// reset accumulables to their initial values
G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
accumulableManager->Reset();

}
void B1RunAction::EndOfRunAction(const G4Run* run)
{
G4int nofEvents = run->GetNumberOfEvent();
if (nofEvents == 0) return;

// Merge accumulables
G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
accumulableManager->Merge();

So, if I only incorporate above lines I would get average time? And do I have to set SetRandomNumberStore(true) ?

19

When you say event, do you mean a radiation absorption → scintillation? Or a single scintillation photon?

For each radiation absorption you should produce many photons. You can store all of these photon-to-cathode time in a file. You can then histogram these times to get a decay profile, which should match closely to the decay parameters you gave your scintillator.

As for “Which time to store out of many event” That depends on your logic. You need to write it in a way that a new time will be written for each photon. If you are only getting one photon per event, you have written your code in a way that perhaps overwrites the time each photon, and only the last one is saved

20

@loydms As you know each radiation absorption would produce many optical photons. So each optical photon has different local time. For a single event(radiation absorption) I am getting a distribution of local time. When I do /run/beamOn 1000 and print G4double time = aStep->GetPostStepPoint()->GetLocalTime(); I am getting a list of 1000 local time.
So my question was is this the average time am I printing ? Because there are many optical photon produced in each event ,in that sense I quoted average.

If this is not the average time then which time it is?
As you said “perhaps overwrites the time each photon, and only the last one is saved”. But I want to print the average time of those optical photon which produce in a single event(radiation absorption).

This is my code :
const G4String particleName
= aStep->GetTrack()->GetDefinition()->GetParticleName();
if (particleName == “opticalphoton”){
G4double time = aStep->GetPostStepPoint()->GetLocalTime();
fEventAction->absTime = time;

21

So if you run 1000 events (electrons? gammas?) you should get something like 100,000+ photons, depending on your energy, as you know. Since you are printing 1000 values, it seems you are only printing 1 time per event, not a time per photon.

If it were an average time, you would see all of the times about the same, while if it was just the last photon absorbed it would be a distribution similar to the decay time of the scintillator.

const G4String particleName
= aStep->GetTrack()->GetDefinition()->GetParticleName();
if (particleName == “opticalphoton”){
G4double time = aStep->GetPostStepPoint()->GetLocalTime();
fEventAction->absTime = time;

So at each step of an optical photon, you are giving your event action 1 time. In fact, you may be sending your event action multiple times per photon, as you don’t have a condition that specifies that your photon is in your detector. I am guessing that you have your event action write the time in the deconstructor. So for a single event, you have your stepping action sending your variable time to event action many many times, but only the last value is printed to your file.

The best way to do this is to histogram the times when the photons are detected. You can also have an array and print all your values, but you will have a file with 100,000+ lines

1 Like
23

Based on your code you do make sure that your photon is in your sensitive volume

if (currentPhysicalName == “Cathode”)

I would also kill your photon if physical name is cathode. This ensures that you will not have a photon scatter in your cathode and then record more times.
I see that you set a your absTime variable in your fEventAction, but I believe EventAction is called when you crate your event, and the destructor is called after the event. So you are constantly changing the variable, and then presumably saving the value to a file after your event (I’d have to see your eventAction to see how this is handled)

What you want to do is store all of these values. You can print them via stepping action, but you will have a very long file. The correct way to do this is histogram the data. This can be done via ROOT

You can also create a histogram yourself. Create an array with the number of bins you want (i.e. if you want 1ns resolution over 200 ns you will have 200 bins) and then check which bin each time falls into using floor division, and increment the array index.

1 Like
25

I am expecting to have like this

Then you will have to write your code to output this. You will have to have something in your stepping action to write to the file each time a photon is absorbed