Some of the events generate 0 Scintillation output (memory issue?)

abhinabmuk · February 29, 2024, 12:06pm

11.1.1 GEANT4 version
Linus OS

Hi,

So I modified the WLS example that consists of WLS plugged in the plastic scintillator wrapped in tyvek with copied sensitive detector with copy number (from the WLS example). I am shooting a muon beam of 4 GeV right in front of the plastic scintillator and it seems to be working perfectly apart from some events that are “disappearing”. For example, I have attached the scintillation output file which shows after certain events (exactly 75) it gives 0 scintillation output and the program rushes through those events. I believe it may be due to stacking action issues or memory loss since I tested both the physics with muon decay off. I also tested with and without wrapping material and it yields the same results. I have included my Stacking Action.cc code below.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
/// \file optical/wls/src/WLSStackingAction.cc
/// \brief Implementation of the WLSStackingAction class
//
//
#include "WLSStackingAction.hh"

#include "G4OpticalPhoton.hh"
#include "G4RunManager.hh"
#include "G4Track.hh"


#include <iostream>
#include <fstream> // Required for file handling

#include "G4AnalysisManager.hh" // Ntuples update
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

WLSStackingAction::WLSStackingAction()
  : fPhotonCounter(0)
  , fScintillationCounter(0)
  , fCerenkovCounter(0)

{}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

WLSStackingAction::~WLSStackingAction() {}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4ClassificationOfNewTrack WLSStackingAction::ClassifyNewTrack(
  const G4Track* aTrack)
{
  G4ParticleDefinition* particleType = aTrack->GetDefinition();

  // keep primary particle
  if (aTrack->GetParentID() == 0) { return fUrgent; }
  
  if(particleType == G4OpticalPhoton::OpticalPhotonDefinition())
  {
    // keep optical photon
    ++fPhotonCounter;
    return fUrgent;
  }


 //stack and count optical photons 
  if(aTrack->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition())
  {  // particle is optical photon
 //   ++fPhotonCounter; // total optical photons produced
    if(aTrack->GetParentID() > 0)
    {  // particle is secondary
      if(aTrack->GetCreatorProcess()->GetProcessName() == "Scintillation")
        ++fScintillationCounter;
      else if(aTrack->GetCreatorProcess()->GetProcessName() == "Cerenkov")
        ++fCerenkovCounter;
    }
  }


  return fUrgent;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void WLSStackingAction::NewStage()
{
  // G4cout << "Number of optical photons produces in this event : "
  //        << fPhotonCounter << G4endl;

// Open a file stream for writing
std::ofstream outputFile("scintillation_output.txt", std::ios::app);

// Check if the file is successfully opened
if (outputFile.is_open()) {
    // Redirect the output to the file stream
    outputFile << fPhotonCounter << std::endl;

    // Close the file stream
    outputFile.close();
} else {
    // Print an error message if the file cannot be opened
    G4cout << "Error: Unable to open the output file." << G4endl;
}
  
 auto analysisManager = G4AnalysisManager::Instance();
 analysisManager->FillNtupleDColumn(4,fPhotonCounter);
 analysisManager->AddNtupleRow();

         
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void WLSStackingAction::PrepareNewEvent() { 
  fPhotonCounter = 0; 
  fScintillationCounter = 0; 
  fCerenkovCounter =0;
  }

I ran 100 events for 3 runs and it seems to get stuck randomly at certain events and then run again. I am shooting the muon beam pretty close to the plastic scintillator tile wrapped in tyvek.

I have attached the DetectorConstruction.cc file
WLSDetectorConstruction.cc (36.4 KB)
for reference as well.

scintillation_output.txt (1.5 KB)

dsawkey · March 11, 2024, 9:41pm

First, it’s not random–you’ve done 3 runs of 100 events, each with the same initial seed, and the results (including zeros) are the same each time.

Are you using multiple threads? The code to write to the file is not thread-safe. Try running only one thread and see if you get the same result. Also, the analysis manager is creating a histogram of number of photons–what does that look like?

Next, try with verbose on:
/tracking/verbose 1
And make sure there are no photons generated. If there is too much output, determine the random number seed for each event, then restart the simulation with the seed that gives a bad event (either wls or LXe optical example does this, and many others).

abhinabmuk · March 12, 2024, 7:20pm

Hello dsawkey,

Thank you for the feedback. I tested both multi-thread and single thread although single-thread events tend to ‘last longer’ till until I get 0s again and the program rushes through the events. I also merged the Ntuples to overcome the thread-safety following the GEANT4 documentation using

analysisManager->AddNtupleRow(0);

Instead of histograms, I have just ROOT files and it shows 0 scintillations in events that GEANT4 quickly processes/skips through. I have attached it below.

I removed the WLS and holes and just had the Scintillator wrapped in Tyvek and ran 1000 events, although it took time it didn’t crash leading me to speculate maybe optical interactions in WLS are causing issues.

Another clarification question this code line in wlsdetectorconstruction.cc:

logicWLSfiber->SetUserLimits(new G4UserLimits(DBL_MAX, DBL_MAX, 10. * ms));

Does this mean the maximum possible step and the maximum total track length are both DBL_MAX? If I were to reduce these numbers would this hamper the accuracy of the model?

I will also do the tracking as you suggested, I will edit out the Optical Physics addition for that since I don’t want to deal with optical photons.

abhinabmuk · March 13, 2024, 12:28pm

I tried with /tracking/verbose 1 and for the seeds it was giving error, it ran smoothly. I did notice a lot of electrons being tracked/generated from 4 GeV muon hitting materials … I am strongly suspecting optical process overloading.

abhinabmuk · March 19, 2024, 4:15pm

Hi, I think I might have narrowed the issue and it is associated with a sensitive detector. If I don’t place the sensitive detector, events run perfectly.

Once I place the sensitive detector it causes issues with scintillator photons produced mainly the issue seems to be resolved with these couple lines of code:

std::vector<G4double> refl_mppc(60, fMPPCReflectivity);  
photonDetSurfaceProperty->AddProperty("REFLECTIVITY", p_mppc, refl_mppc);

If I remove this it works great but I can only see 1 output of duplicated volume of sensitive detector i.e it only shows 1 detector ID.

Furthermore, If I try removing

new G4LogicalSkinSurface("PhotonDetSurface", logicPhotonDet,photonDetSurface);

This causes issues as it is unable to detect photons in sensitive detector [see image below].

Am I duplicating sensitive detector correctly ? :

 auto BialkaliCathode = FindMaterial("BialkaliCathode");


// Physical Construction
G4VSolid* solidPhotonDet = nullptr;

solidPhotonDet = new G4Box("PhotonDet", fMPPCHalfL, fMPPCHalfL, fMPPCZ);
G4LogicalVolume* logicPhotonDet = new G4LogicalVolume(solidPhotonDet, BialkaliCathode, "PhotonDet_LV");

// G4double xPMT = -12.5 * cm;
// G4double yPMT = 12.5 * cm; 
// G4double zPMT = fPDPos;
// G4PVPlacement* physicalDetector = new G4PVPlacement(0, G4ThreeVector(xPMT, yPMT, zPMT), logicPhotonDet, "PhotonDet_LV", fLogicWorld, false, 0); // copy number of PMTs 

G4PVPlacement* physicalDetector;

for (G4int i = 0; i < 2; ++i) {
    for (G4int j = 0; j < 2; ++j) {
        G4double xPMT = -12.5 * cm + pmtSpacing * i;
        G4double yPMT = 12.5 * cm - pmtSpacing * j;
        G4double zPMT = fPDPos;
        physicalDetector = new G4PVPlacement(0, G4ThreeVector(xPMT, yPMT, zPMT), logicPhotonDet, "PhotonDet", fLogicWorld, false, i * 2 + j); // copy number of PMTs 
    }
}

dsawkey · March 20, 2024, 2:47pm

No. There’s only one G4PVPlacement*, which gets over-written each time through the loops. You want an array of G4PVPlacement pointers.

abhinabmuk · March 21, 2024, 4:06pm

Thanks a lot dsawkey but the problem is still persistent even if I use 1 detector. It seems something associated with optical photons interacting with sensitive detectors. I have attached an image below I got a message in the terminal.

I also noticed I had a lot of urgent stacks if I ran fewer events [and it runs fine] ?

abhinabmuk · March 22, 2024, 3:37pm

Hello,

Another update. It seems the issue originates from tyvek wrapping/coating which is unlcear to me. If I use the default coating it works:

  density = 1.52 * g / cm3;

  fCoating = new G4Material("Coating", density, ncomponents = 2);

  fCoating->AddMaterial(TiO2, fractionmass = 15 * perCent);
  fCoating->AddMaterial(fPolystyrene, fractionmass = 85 * perCent);

but as soon as I replace the wrapping with tyvek it seems to be causing issues :

  //------------------------------------------------------
  //-------------Tyvek Wrapping--------------------------
  //------------------------------------------------------

  elements.push_back("C");
  natoms.push_back(2);
  elements.push_back("H");
  natoms.push_back(4);

  density =0.36 * g / cm3;

  fTyvek = fNistMan->ConstructNewMaterial("Tyvek", elements, natoms, density);

  elements.clear();
  natoms.clear();

Any hints as to why?

dsawkey · March 22, 2024, 8:29pm

No idea. It doesn’t look like the problem originates with optical processes, though, as there are so many things going wrong. I’d say, go back to the original WLS example and modify it step by step. When you find exactly what change causes the break, post back.