Filling ntuples inside IF statements

leonardoghizoni · February 18, 2024, 1:03pm

Dear all,

In my code I’m trying to record any e- or e+ and respective energies that enter a CCD volume in my geometry. My simulation makes usage of GPS, and I shoot omnidirectional monoenergetic electrons, getting the record of secondaries that reach my volume. For that I use G4SteppingManager, and try to record everything in a ROOT file via g4tools.
In my RunAction.cc file I set the following to create ntuples:

    // Create analysis manager
    auto analysisManager = G4AnalysisManager::Instance();
    G4cout << "Using " << analysisManager->GetType() << G4endl;

    analysisManager->SetVerboseLevel(1);

#ifdef G4MULTITHREADED
    analysisManager->SetNtupleMerging(true);
#endif

    // Create Ntuples - ROOT
    // Electron ntuple (id = 0)
    analysisManager->CreateNtuple("e-", "e- quantities");
    analysisManager->CreateNtupleIColumn("count"); // column id = 0
    analysisManager->CreateNtupleDColumn("eKin");  // column id = 1
    analysisManager->CreateNtupleDColumn("eTot");  // column id = 2
    analysisManager->FinishNtuple();

    // Positron ntuple (id = 1)
    analysisManager->CreateNtuple("e+", "e+ quantities");
    analysisManager->CreateNtupleIColumn("count"); // column id = 0
    analysisManager->CreateNtupleDColumn("eKin");  // column id = 1
    analysisManager->CreateNtupleDColumn("eTot");  // column id = 2
    analysisManager->FinishNtuple();

And in my SteppingManager.cc I add the following to record e-/e+ with respective energies:

    auto analysisManager = G4AnalysisManager::Instance();

    auto volPreStep = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume();
    auto volPostStep = step->GetPostStepPoint()->GetTouchableHandle()->GetVolume();

    if(volPreStep != fGeometry->GetCCD() && volPostStep == fGeometry->GetCCD()){
        if(step->GetTrack()->GetDefinition()->GetParticleName() == "e-"){
            analysisManager->FillNtupleIColumn(0, 0, 1);
            analysisManager->FillNtupleDColumn(0, 1, step->GetTrack()->GetKineticEnergy());
            analysisManager->FillNtupleDColumn(0, 2, step->GetTrack()->GetTotalEnergy());
            analysisManager->AddNtupleRow(0);
        }
        else if(step->GetTrack()->GetDefinition()->GetParticleName() == "e+"){
            analysisManager->FillNtupleIColumn(1, 0, 1);
            analysisManager->FillNtupleDColumn(1, 1, step->GetTrack()->GetKineticEnergy());
            analysisManager->FillNtupleDColumn(1, 2, step->GetTrack()->GetTotalEnergy());
            analysisManager->AddNtupleRow(1);
        }
    }

The simulation runs with no problem, but my ntuples are completely empty afterwards. I of course add the Write() and CloseFile() methods in EndOfRunAction of my RunAction.cc.
I tried adding the AddNtupleRow methods outside the if statements like below, but then the ntuples are filled at every step (with zeros when my if statements are not fulfilled), and my ROOT file becomes obscenely huge.

    if(volPreStep != fGeometry->GetCCD() && volPostStep == fGeometry->GetCCD()){
        if(step->GetTrack()->GetDefinition()->GetParticleName() == "e-"){
            analysisManager->FillNtupleIColumn(0, 0, 1);
            analysisManager->FillNtupleDColumn(0, 1, step->GetTrack()->GetKineticEnergy());
            analysisManager->FillNtupleDColumn(0, 2, step->GetTrack()->GetTotalEnergy());
        }
        else if(step->GetTrack()->GetDefinition()->GetParticleName() == "e+"){
            analysisManager->FillNtupleIColumn(1, 0, 1);
            analysisManager->FillNtupleDColumn(1, 1, step->GetTrack()->GetKineticEnergy());
            analysisManager->FillNtupleDColumn(1, 2, step->GetTrack()->GetTotalEnergy());
        }
    }
    analysisManager->AddNtupleRow(0);
    analysisManager->AddNtupleRow(1);

Any help on how to fill my ntuples inside the if statements only would be more than appreciated.

Thanks in advance.

_Geant4 Version:_10.7.3
_Operating System:_GNU/Linux Slackware 15 (conda environment)
_Compiler/Version:_GCC 10.4
_CMake Version:_3.21.4

JDecunha · February 22, 2024, 8:01pm

Hi @leonardoghizoni,

Have you confirmed that any particles are meeting the requirement of?

if(volPreStep != fGeometry->GetCCD() && volPostStep == fGeometry->GetCCD())

My initial suspicion looking at your code is that the “e-” and “e+” filtering if statements are fine, but the if statement before it might not be working as intended. You could use g4cout/std::cout to verify that particles are passing through the first if statement.

If you find that it is the 1st if statement that isn’t working another possibility would be to try calling:

auto volPreStepName = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetName();
auto volPostStepName = step->GetPostStepPoint()->GetTouchableHandle()->GetVolume()->GetName();
if (volPreStepName != “YourCCDPhysicalVolumeName” && volPostStepName == “YourCCDPhysicalVolumeName”)

I’ve just suggested the above approach because I don’t know how you’re going about getting the fGeometry pointer and if there may be some issues there.

I hope this is helpful for you,

Joseph

leonardoghizoni · February 23, 2024, 10:31am

Dear @JDecunha,

Thanks for your reply and input.

I can confirm that the particles are meeting the requirement of the first IF filtering. I’ve added a G4cout just to make sure via:

    if(volPreStep != fGeometry->GetCCD() && volPostStep == fGeometry->GetCCD()){
        if(step->GetTrack()->GetDefinition()->GetParticleName() == "e-"){
            G4cout << "Electron eKin = " << G4BestUnit(step->GetTrack()->GetKineticEnergy(), "Energy") << G4endl;
            G4cout << "Electron eTot = " << G4BestUnit(step->GetTrack()->GetTotalEnergy(), "Energy") << G4endl;
            analysisManager->FillNtupleIColumn(0, 0, 1);
            analysisManager->FillNtupleDColumn(0, 1, step->GetTrack()->GetKineticEnergy());
            analysisManager->FillNtupleDColumn(0, 2, step->GetTrack()->GetTotalEnergy());
            analysisManager->AddNtupleRow(0);
        }
        else if(step->GetTrack()->GetDefinition()->GetParticleName() == "e+"){
            G4cout << "Positron eKin = " << G4BestUnit(step->GetTrack()->GetKineticEnergy(), "Energy") << G4endl;
            G4cout << "Positron eTot = " << G4BestUnit(step->GetTrack()->GetTotalEnergy(), "Energy") << G4endl;
            analysisManager->FillNtupleIColumn(1, 0, 1);
            analysisManager->FillNtupleDColumn(1, 1, step->GetTrack()->GetKineticEnergy());
            analysisManager->FillNtupleDColumn(1, 2, step->GetTrack()->GetTotalEnergy());
            analysisManager->AddNtupleRow(1);
        }
    }

And I get what I expected:

### Run 0 start.
... create file : /home/leogh/Desktop/MIS.root - done
... open main analysis file : /home/leogh/Desktop/MIS.root - done
G4WT5 > ### Run 0 start.
G4WT3 > ### Run 0 start.
G4WT2 > ### Run 0 start.
G4WT1 > ### Run 0 start.
G4WT0 > ### Run 0 start.
G4WT4 > ### Run 0 start.
G4WT2 > Electron eKin = 256.969 MeV
G4WT2 > Electron eTot = 257.48 MeV
G4WT2 > Electron eKin = 148.012 keV
G4WT2 > Electron eTot = 659.011 keV
G4WT5 > Electron eKin = 2.29518 MeV
G4WT5 > Electron eTot = 2.80618 MeV
G4WT5 > Positron eKin = 41.94 MeV
G4WT5 > Positron eTot = 42.451 MeV
G4WT5 > Positron eKin = 13.2154 MeV
G4WT5 > Positron eTot = 13.7264 MeV
G4WT5 > Electron eKin = 30.0738 MeV
G4WT5 > Electron eTot = 30.5848 MeV
G4WT4 > Electron eKin = 336.677 MeV
G4WT4 > Electron eTot = 337.188 MeV
G4WT4 > Electron eKin = 4.50274 keV
G4WT4 > Electron eTot = 515.502 keV
G4WT4 > Electron eKin = 3.66182 keV
G4WT4 > Electron eTot = 514.661 keV
G4WT4 > Positron eKin = 113.339 MeV
G4WT4 > Positron eTot = 113.85 MeV

And so on…

As I mentioned, the ntuples are filled if I place analysisManager->AddNtupleRow(0) and analysisManager->AddNtupleRow(1) outside the IFs, but then they are filled at every step, and with zero in case the IF statements are not fulfilled. Since I need to simulate 1e8 or 1e9 for good statistics, my ROOT file becomes ridiculously huge (I tried a single run, and it got ~300GB!), and it’s not feasible
to follow that approach. I simply don’t understand why the analysisManager->AddNtupleRow() doesn’t work inside the IF statement.

But thanks again for the help.

JDecunha · February 23, 2024, 8:51pm

I agree that it is extremely puzzling that placing the AddNtupleRow commands within the if statement is preventing the data from outputting. Because this behavior is strange I am trying to imagine any obscure reasons it might be occuring. For completeness, would you mind sharing your macro file if it includes any analysis commands, and the code you use to write and close the file?

I believe that ROOT (and Geant’s implementation of AnalysisManager) has functionality that will automatically write an output file after an amount of data has been filled which exceeds the buffer size. If, for some reason your code to write the file was not being called properly, it might be the case that when you put the AddNtupleRow statements outside of the if statement, then the amount of data written is enormous and exceeds the buffer size, causing your data to be automatically written. A way you might test if you write command is being called is to call AnalysisManager::Write() within the if statement at every loop of the SteppingAction (not an efficient solution, but we are troubleshooting).

Joseph

leonardoghizoni · February 24, 2024, 12:12pm

@JDecunha thanks again for all the help. Really much appreciated.

I’ve tried attaching a zip with a minimally reproducible example so you can run it, but it seems I’m not allowed to attach zip here. So I added it in my google Drive, in the following link: https://drive.google.com/file/d/13J8OvlwvVSxg1I9rt_gCBY7MHqoBCTFQ/view?usp=drive_link

I had to change some things in regard to the geometry, as I can’t share it outside, so you’ll find some GDML function leftovers there. Nevertheless, it is fully reproducible.
You need to run it like this: “./mis steppingManager_800MeV.mac /path/pathTo”, where path is where you want to store the ROOT file.
I’ve added the commented AddNtupleRow outside the IF statements, so you can uncomment them and see that the ROOT file is filled, but then it is also filled with zeroes. In this last situation, you can check that the ROOT file correctly stores the data by comparing the G4cout outputs from the simulation with the following python macro:

import ROOT
import numpy as np 

rdf = ROOT.RDataFrame("e-", "MIS.root")
eKin = np.array(rdf.Take['double']("eKin").GetValue())

print(len(eKin))

eKin = eKin[eKin>0]
print(eKin)

As you’ll be able to confirm, I need to ignore all zeroes in the numpy array that were filled when the AddNtupleRow is placed outside the IF statement. When the AddNtupleRow is placed inside the IF statements, the arrays are completely empty.
In regard to calling AnalysisManager::Write() in each IF statement loop, I experienced a bunch of errors
Thank you so much again for the help.

JDecunha · February 28, 2024, 12:08am

I’m running on Geant4 11.2, so I had to modify the #include calls since some of the header files for G4AnalysisManager have changed names compared to the Geant verison you’re using. I’m happy to share the modified code with you if you need. I believe I have your software running correctly.

I believe that the main issue was the:

if(!IsMaster()) { break; }

call in the EndofRunAction. This seems like it is preventing the G4AnalysisManager instances running on each thread from outputting their data. I modified the section as follows:

void MISRunAction::EndOfRunAction(const G4Run* aRun)
{
auto analysisManager = G4AnalysisManager::Instance();

    if(IsMaster())
    {
        //All of the MISRun Output code here that should only run on the main thread
    }

analysisManager->Write();
analysisManager->CloseFile();
}

I hope this is helpful and gets the software running on your end as well!

Joseph

mkelsey · February 28, 2024, 1:23am

@JDecunha If your embedded comment is correct, shouldn’t that be “if (IsMaster())”?

JDecunha · February 28, 2024, 1:47am

@mkelsey You’re absolutely right. I’ve corrected the typo above.

Additionally, @leonardoghizoni, if you are seeing segfaults at the termination of your code you can remove the statement in the destructor of MISRunAction. I don’t believe you need to explicitly free the memory held by G4AnalysisManager.

leonardoghizoni · February 28, 2024, 7:55am

@JDecunha,

Oh my…!! Thank you so much!

I must say I’d probably never find the issue by myself. The IsMaster() function is a leftover from when I started building the code months ago trying to get global quantities based on some of the Geant4 examples, and just left it there for eventual future references.
The issue is solved, and now I have a ROOT file containing only what I wanted. Much appreciated all your help.
Also, thanks @mkelsey for the input as well

Cheers,
Leonardo.

system · March 6, 2024, 7:55am

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