Issue with optical photons crossing a boundary

Please fill out the following information to help in answering your question, and also see tips for posting code snippets. If you don’t provide this information it will take more time to help with your problem!

Geant4 Version: 11.2.2
Operating System: AlmaLinux 9
Compiler/Version: 13.1.0
CMake Version: 3.26.2

Hi,

I’m simulating a rectangular scintillating crystal with two SiPMs attached to the bases of the crystal. I activated cerenkov and scintillation processes in the crystal, the SiPMs are simply two G4Box made of Silicon. I count the number of optical photons produced in the crystal and the deposited energy in the crystal to see if they are proportional and everything seems to work. I count the number of photons in the SteppingAction in this way:

void SteppingAction::UserSteppingAction(const G4Step* step){
  auto theTrack = step->GetTrack();
  G4ParticleDefinition* particleType = theTrack->GetDefinition();
  G4StepPoint *thePrePoint  = step->GetPreStepPoint();
  G4StepPoint *thePostPoint = step->GetPostStepPoint();
  G4VPhysicalVolume *thePrePV = thePrePoint->GetPhysicalVolume();
  G4VPhysicalVolume *thePostPV = thePostPoint->GetPhysicalVolume();

  if (particleType == G4OpticalPhoton::OpticalPhotonDefinition()) {
    if (nStep == 1 && G4StrUtil::contains(thePrePVName, "Calorimeter")) {
      if (processName == "Scintillation") {
        fEventAction->tot_phot_scint_calo += 1;
       }
     }
  }
}

Here’s the plot of the energy deposit and the number of photons

Then I want to count the number of photons the reach the SiPMs on the bases of the crystal, and in this case I seem to get problems. Perhaps it’s important to say that the base of the crystal is 12 x 12 mm, while the SiPMs are 3 x 3 mm and 6 x 6 mm. I count the number of photons in this way (for example the 6 x 6 mm SiPM is called SiPM_C):

if (particleType == G4OpticalPhoton::OpticalPhotonDefinition()) {
  if (G4StrUtil::contains(thePrePVName, "Calorimeter") && G4StrUtil::contains(thePostPVName, "SiPM_C")) {
    fEventAction->tot_phot_scint_SiPM_C ++;
  }
}

I expect to collect less photons than the total produced, but still to see a correlation with the deposited energy. Instead I get this:

The number of photons detected in each event follows this distribution that doesn’t make much sense to me:

All these results seem unphysical to me but I can’t understand where is the error in my code.
Edit: There are other collegues working with me on this so they can also provide details of the code in the replies if necessary.

I am not sure exactly what the issue is but one issue I could see is that you would start multiply counting any optical photons that are reflected or scattered between the two surfaces. That could cause issues like you are seeing.

One way to prevent this would be to add a flag to the photons track when it first enters the silicon region and pass this info to all of its secondaries. So when the particle or any of its descendants cross the boundary you can check if this flag is set, and only add photons that don’t have this flag set. Or, depending on what you are doing, you could just kill the track entirely once it enters the silicon. Either way this example might help.

Unrelated but you might want to just pass the logical volumes to the stepping action from the detector construction and then check if you are in those volumes. That will avoid having to use the string utilities that might slow down your simulation.

1 Like

hi, thank you for your reply.

I’m simulating a calorimeter, so I want the number of optical photons to be proportional to the deposited energy in the calorimeter. The issue is that photons reaching the sipms are completely uncorrelated to the energy deposit in the crystal, I can’t think of a physical reason why it should happen, so I assume there’s a bug in the code.

I tired to run the code killing the tracks that enter the silicon but nothing changed at all.

Thank you for the tip about the logical volumes, much appreciated.

If the count is right but the energy is wrong (and likely “missing”) then I’d explicitly set deexcitation, auger, and pixe on and off to see what impact that has since these can be very relevant for solid state detectors. Fast way might be:

G4Region* siliconLogReg = new G4Region(“siliconLogReg”);
siliconLogReg->AddRootLogicalVolume(yourSiliconBoxLogVolume);

G4EmParameters* emOptions = G4EmParameters::Instance();
emOptions->SetDeexActiveRegion(“siliconLogReg”, true, true, true); // deexc., auger, pixe

If nothing changes from activating (or deactivating) those three processes there might be a physics list issue.

2 Likes

One thing to check is how you insert the primary particle. If the incident beam is zero thickness and always the same angle, you may see miracles never occurring in the lab.