Provide string to the SteppingAction class in Hadr04 to set a filename

Hello! In Hadr04, I want to save some extra output. I don’t think this is the smartest way to do it but I want to record the time and positions of the neutron captures, and was thinking of a text based output (I do not know how to add an additional histogram to Hadr04 but if that could be explained that would be awesome!).

For example

   // Positions of nCapture

  G4String process_name = aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();

   if(!process_name.compare("nCapture")){
     G4double time_finish=aStep->GetPostStepPoint()->GetGlobalTime();
     G4ThreeVector position_finish=aStep->GetPostStepPoint()->GetPosition();
     G4int eventid=G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
     ofstream capturesfile;
     capturesfile.open (<some_file_name>, ios_base::app);
     capturesfile << time_finish << ",";
     capturesfile << position_finish[0] << ",";
     capturesfile << position_finish[1] << ",";
     capturesfile << position_finish[2] << ",";
     capturesfile << eventid <<"\n";

     capturesfile.close();

How do I allow the user to set the <some_file_name> equal some some string.
I.e., for a given run, all steps for all events will add to the same run file called <some_file_name>? Thanks!

The geant4-way would be to use AnalysisManager, and fill the values into an ntuple:
https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Analysis/managers.html#ntuples

With ntuple-merging you get them in one file, and you also have no issues with race conditions for file access in case you are using multithreading.

Thanks for the reply and the advice! Just to clarify, I should add the

analysisManager->SetNtupleMerging(true);

to make sure that the Ntuple sits in one file? If so, where to add? E.g. If I define for example in the HistoManager.cc (it appears analysisManager is already in there)

analysisManager -> CreateNtuple("NeutronCapture","Neutron Capture N tuple");
analysisManager->SetNtupleMerging(True);
analysisManager->SetNtupleFileName(0,"NeutronCaptures.root");
analysisManager->CreateNtupleIColumn("RunID");
analysisManager->CreateNtupleIColumn("EventID");
analysisManager->CreateNtupleDColumn("t");
analysisManager->CreateNtupleDColumn("x");
analysisManager->CreateNtupleDColumn("y");
analysisManager->CreateNtupleDColumn("z");
analysisManager->CreateNtupleDColumn("captureEnergy");
analysisManager->FinishNtuple();

does this work? How do I then tell SteppingAction.cc to fill these entries? Thanks again!

I guess this should be in the UserSteppingAction? Not sure how to get each run ID though

// if it is a capture, get the step end time
  G4double t_finish = aStep->GetPostStepPoint()->GetGlobalTime();
  G4ThreeVector position_finish = aStep->GetPostStepPoint()->GetPosition();
  G4int eventid = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
  //Fill the Run ID column
  G4AnalysisManager::Instance()->FillNtupleIColumn(0,);
  //Fill the event column for the run
  G4AnalysisManager::Instance()->FillNtupleIColumn(1,eventid);
  //Fill the capture time
  G4AnalysisManager::Instance()->FillNtupleDColumn(2,t_finish);
  //Fill the x position of the Capture
  G4AnalysisManager::Instance()->FillNtupleDColumn(3,position_finish[0]);
  //Fill the y position of the capture
  G4AnalysisManager::Instance()->FillNtupleDColumn(4,position_finish[1]);
  //Fill the z position of the capture
  G4AnalysisManager::Instance()->FillNtupleDColumn(5,position_finish[2]);
  //Fill the energy of the capture
  G4AnalysisManager::Instance()->FillNtupleDColumn(6,ekin);

yes, I believe so. Try what happens if you dont add this line :slight_smile:
I have placed it into the constructor of the G4UserRunAction.

does it?

your next post answers this, but don’t forget the line

  G4AnalysisManager* man = G4AnalysisManager::Instance();
 ...
  man->AddNtupleRow();

and the first column is missing the runid.

Also it could improve the code to store the G4AnalysisManager::Instance() like in the example as you access it multiple times in a row. Don’t know if that makes a difference in performance, but it is more aesthetically, codewise :wink:

OK I’m gonna try implementing (spooky times ahead!), Back soon! Thank you so much!

To get the runID do I just say

G4RunManager::Instance()->GetRunID();

or something like this? I read in the manual it looks like the user can set it but it wasn’t clear how to do this in the C++ or via macro (Sorry for my confusion! :slight_smile: )

#include "G4Run.hh"
...
G4RunManager::GetRunManager()->GetCurrentRun()->GetRunID();

Legendary - OK, gonna go try it out and back soon, thanks so much @weller !!

So something interesting happened! Though I’m a bit unsure of the units?





E.g. why is the capture energy 0?

print the values in the run to check. in the snipped you posted there was no definition for „ekin“?

Ah my bad, ekin is this:


 if (aStep->GetTrack()->GetTrackID() == 1) {
    G4double ekin  = endPoint->GetKineticEnergy();

default unit is MeV, maybe the values are thus very close to zero on that scale. System of units — Book For Application Developers 11.1 documentation

ahh gotcha! I’m just trying to figure out how to load the values xD

Since the Geant4 units are not carried into the N-tuple, you might find it convenient to explicitly choose units for storing the values. In our experiment (a low-energy dark matter search), we have decided that our analysis code will use units of eV for energy, m for length, and ns for timestamps. So in our N-tuples, we save out

   aMan->FillNtupleDColumn(5, ekin/eV);
   aMan->FillNtupleDColumn(10, hitTime/ns);
   aMan->FillNtupleDColumn(11,hitX/meter);

and so on. I have also seen people who do the above, and have a separate set of columns where they store the unit value as well, e.g.,

  aMan->FillNtupleDColumn(23, eV);

Oh wow, that is great information - let me try to implement it and I’ll be back soon with new plots! Thanks!

So I gave it a go (and in Hadr04.cc I added the line G4VModularPhysicsList* phys = new FTFP_BERT_HP; to try and include high energy neutron beam) and I got this… it looks kinda weird





It gave messages appearing like it couldn’t find carbon scattering data for FTFP_BERT_HP (i changed my detector to include some fraction of carbon)

NeutronHP: /Capture file for Z = 6, A = 12 is not found and NeutronHP will use BLAH
NeutronHP: /Elastic file for Z = 6, A = 12 is not found and NeutronHP will use BLAH
NeutronHP: /Inelastic file for Z = 6, A = 12 is not found and NeutronHP will use BLAH
etc

where BLAH just refers to my directory variables (geant found them ok :slight_smile: )

Two things:

Those messages mean it didn’t find separate specific-isotope files (for C-12, in this case). I believe it falls back to using the generic “natural mix” files, which should have been mentioned in the BLAH strings.

You had asked earlier

and your plot shows that it’s identically zero, not just “small.” If you’re filling it with

where “endPoint” is the PostStepPoint, then I think it will always be zero for a capture (capture means “stopped and killed”). If you want to know the energy it came into the nucleus with, use the preStepPoint.

Oh got it - I guess those files cannot go above 10s of MeV right? Is there a way to use the higher energy? Thanks again!