Support for G4Accumulable<std::vector<type>>

I’ve been looking into G4Accumulables. If I understand correctly, G4Accumulable<T> only supports type T that are single objects. It seems to me, that the utility of a G4Accumulable could be greatly increased if it could support std::vectors of objects. For example, one could accumulate a vector of G4doubles that represent the amount of energy deposited at each location in a scoring grid.

I have written code (see below) that implements a std::vector of accumulables. The code is incomplete because it does not implement the different possible G4MergeModes. The actual use of G4VectorAccumulable is a bit messy and requires a static_cast every time it is pulled from the G4AccumulableManager. This messiness could be tidied up if G4AccumulableManager were modified to support G4VectorAccumulable natively.

My questions are as follows:
1.) Am I right to understand that the existing functionality of G4Accumulable does not support the use of std::vectors or other container types?
2.) If my understanding of point 1 is correct, would it be of value to the community to incorporate such functionality into Geant4?

G4VectorAccumulable.hh

#pragma once

//Geant4
#include "globals.hh"
#include "G4VAccumulable.hh"
//STD
#include <vector>

// Not sure if Geant4 is using C++20 features yet. But we could use a concept to restrict the G4VectorAccumulable to only numeric types. 
// template<typename T>
// concept AccumulableType = std::integral<T> || std::floating_point<T>;

template <typename AccumulableType>
class G4VectorAccumulable : public G4VAccumulable
{
	public:
	    G4VectorAccumulable(const G4String& name, const G4int& size) 
	    	: G4VAccumulable(name), _accumulableVector(size, AccumulableType()) //Default initialize the vector to required size
	    	{ }
	    		 
	    virtual ~G4VectorAccumulable() override {}

	    virtual void Merge(const G4VAccumulable& other) override;
	    virtual void Reset() override;

	    void SetEntry(const G4int& entryNumber, const AccumulableType& entryValue);
	    void AddToEntry(const G4int& entryNumber, const AccumulableType& entryValue);
	    const std::vector<AccumulableType>& GetVector() const { return _accumulableVector; } 

	private:
		std::vector<AccumulableType> _accumulableVector;
};

#include "G4VectorAccumulable.imp"

G4VectorAccumulable.imp

//This project
#include "G4VectorAccumulable.hh"

template <typename AccumulableType>
void G4VectorAccumulable<AccumulableType>::SetEntry(const G4int& entryNumber, const AccumulableType& entryValue)
{
	_accumulableVector[entryNumber] = entryValue;
}

template <typename AccumulableType>
void G4VectorAccumulable<AccumulableType>::AddToEntry(const G4int& entryNumber, const AccumulableType& entryValue)
{
	_accumulableVector[entryNumber] = _accumulableVector[entryNumber] + entryValue;
}

template <typename AccumulableType>
void G4VectorAccumulable<AccumulableType>::Merge(const G4VAccumulable& other)
{
	const G4VectorAccumulable<AccumulableType>& otherVector = static_cast<const G4VectorAccumulable<AccumulableType>&>(other);

	for (int i = 0; i < _accumulableVector.size(); i++) //This assumes vectors being merged have the same size. It shouldn't be possible for them to have different sizes.
	{
		_accumulableVector[i] = _accumulableVector[i] + (otherVector.GetVector())[i];
	}
}

template <typename AccumulableType>
void G4VectorAccumulable<AccumulableType>::Reset()
{
	for (auto& itemInVector : _accumulableVector) 
	{
		itemInVector = AccumulableType(); //Default initialize each item in the vector (i.e. zero it). But don't resize the vector itself.
	}
}

Though Ideas and Requirements is a fine place to post this, I’ve moved it to the persistency category just to get a bit more visibility (e.g. to @ivana)

Dear @JDecunha,

I asked @ivana once something similar and I know how to get an array of G4Accumulables which may not be what you want but in case it is please find the differences that you would need to implement to the example B1 to get this to work.

diff --git a/examples/basic/B1/include/RunAction.hh b/examples/basic/B1/include/RunAction.hh
index 3289905c29..e777357c54 100644
--- a/examples/basic/B1/include/RunAction.hh
+++ b/examples/basic/B1/include/RunAction.hh
@@ -34,6 +34,8 @@
 #include "G4Accumulable.hh"
 #include "globals.hh"
 
+#include <array>
+
 class G4Run;
 
 /// Run action class
@@ -59,6 +61,7 @@ class RunAction : public G4UserRunAction
   private:
     G4Accumulable<G4double> fEdep = 0.;
     G4Accumulable<G4double> fEdep2 = 0.;
+    std::array<G4Accumulable<G4double>, 2> fEdepArray = {0., 0.};
 };
 
 }
diff --git a/examples/basic/B1/src/RunAction.cc b/examples/basic/B1/src/RunAction.cc
index 9fa4a5212b..aae34a232d 100644
--- a/examples/basic/B1/src/RunAction.cc
+++ b/examples/basic/B1/src/RunAction.cc
@@ -63,6 +63,10 @@ RunAction::RunAction()
   G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
   accumulableManager->RegisterAccumulable(fEdep);
   accumulableManager->RegisterAccumulable(fEdep2);
+
+  for (auto& elem : fEdepArray) {
+    accumulableManager->RegisterAccumulable(elem);
+  }
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -89,6 +93,14 @@ void RunAction::EndOfRunAction(const G4Run* run)
   G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
   accumulableManager->Merge();
 
+  // Print accumulable values
+  G4cout << fEdep.GetValue() << " " << fEdep2.GetValue() << G4endl;
+
+  for (const auto& elem : fEdepArray) {
+    G4cout << elem.GetValue() << " ";
+  }
+  G4cout << G4endl;
+
   // Compute dose = total energy deposit in a run and its variance
   //
   G4double edep  = fEdep.GetValue();
@@ -149,6 +161,8 @@ void RunAction::AddEdep(G4double edep)
 {
   fEdep  += edep;
   fEdep2 += edep*edep;
+  fEdepArray[0] += edep;
+  fEdepArray[1] += edep*edep;
 }
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

I hope this is somehow helpful

/Pico

Hello,

My questions are as follows:
1.) Am I right to understand that the existing functionality of G4Accumulable does not support the use of std::vectors or other container types?

Yes, G4Accumulable is intended to be used with fundamental (primitive) types like int, double, float only

2.) If my understanding of point 1 is correct, would it be of value to the community to incorporate such functionality into Geant4?

There is an implementation of G4VectorAccummulable available in examples/extended/analysis/B1Con, and an example of a map implementation in the Application Developers Guide see section Analysis/Accummulables

It would be of possible to provide the implementation for the standard library collections, like vector, array or map if there is a user request.

Hello @pico and @ivana,

Thank you both for your helpful replies.

I feel a little silly for not thinking of the approach to implement a std::array of accumulables.

The only disadvantage I can see of the std::array approach, is that the user has to hold onto a pointer to the std::array and pass it between the various user actions where the array is accessed. For example, an array of accumulables initialized in a RunAction would have to have its address be passed to the SteppingAction or SensitiveDetector in order for the elements to be filled correctly. Alternatively, if a std::vector or other standard library data member were stored directly they could be pulled and indexed from the G4AccumulableManager without the need to pass pointers around.

Please correct me if my understanding above is incorrect,

Joseph

Hello,

We have added the support for the most frequent std collections (array, vector) in the work plan for 2024, to make it easier for the users.

Best regards,

dear ivana:
What if fEdepArray has 1,000 numbers, do they all need to be written as fEdepArray = {0., 0.}?

Have a look at Array initialization - cppreference.com.

cppreference.com is always a good place to look for such information. I just googled “C++ array initialisation”.

Thank you, dear @allison.
I have tried these methods, for example, std::array<G4Accumulable, 2> fEdepArray = {0.}. I encountered errors when using G4Accumulate . It seems that specifying all members of the array explicitly is necessary.
However, when i tried std::array<G4double, 2> fEdepArray = {0.}, no errors occurred. This discrepancy makes me wonder if there are additional intricacies to the G4Accumulate template that I’m not aware of.
In the end, I decided to abandon the use of G4Accumulate and instead opted to output the results of all threads. I then processed the data using MATLAB.

What if fEdepArray has 1,000 numbers, do they all need to be written as fEdepArray = {0., 0.}?

Hello,

The problem here is caused by the deleted default constructor in G4Accummulable<T>, which prevents form declaring the array without an explicit initialization. I am going to remove this limitation in the development version.

Thank you for pointing at this issue,