Add stepping histograms generation to the Event Action

The stepping action create different histograms to track stepping quantities, such as step kinetic energy, step pseudorapidity,
step length and other quantities. The histograms are merged automatically by each worker thread at the end of the run
and saved in a ROOT file. No dependency on ROOT is introduced.

FullSimLight/include/MySteppingAction.hh

@@ -4,21 +4,126 @@
 
 #include "globals.hh"
 #include "G4UserSteppingAction.hh"
+#include "MyRunAction.hh"
 
 class MyEventAction;
 class G4Step;
 
+
 class MySteppingAction : public G4UserSteppingAction {
 
 public:
 
-  MySteppingAction(MyEventAction* evtact);
+  MySteppingAction(MyEventAction* evtact, MyRunAction* run);
   virtual ~MySteppingAction();
 
   virtual void UserSteppingAction(const G4Step*);
+    
+  // maps to hold info per volume/process/material per particle type
+  //typedef std::map<G4String, TH1*> HistoMap_t;
+  typedef std::map<G4String, int> HistoMap_t;
+  typedef std::map<G4String, HistoMap_t> HistoMapMap_t;
+
+    
+  /// this holds all the data from individual threads that needs to be merged at EoR
+  struct Report
+  {
+    // distributions per volume per particle type
+    HistoMapMap_t histoMapMap_vol_stepSize;
+    HistoMapMap_t histoMapMap_vol_stepKineticEnergy;
+    HistoMapMap_t histoMapMap_vol_postStepKineticEnergy;
+    HistoMapMap_t histoMapMap_vol_stepPseudorapidity;
+    HistoMapMap_t histoMapMap_vol_stepEnergyDeposit;
+    HistoMapMap_t histoMapMap_vol_stepEnergyNonIonDeposit;
+    HistoMapMap_t histoMapMap_vol_stepSecondaryKinetic;
+    HistoMapMap_t histoMapMap_vol_numberOfSteps;
+    HistoMapMap_t histoMapMap_vol_numberOfStepsPerInitialE;
+    HistoMapMap_t histoMapMap_vol_trackLengthPerInitialE;
+    HistoMapMap_t histoMapMap_vol_InitialE;
+
+    // distributions per material per particle type
+    HistoMapMap_t histoMapMap_mat_stepSize;
+    HistoMapMap_t histoMapMap_mat_stepKineticEnergy;
+    HistoMapMap_t histoMapMap_mat_postStepKineticEnergy;
+    HistoMapMap_t histoMapMap_mat_stepPseudorapidity;
+    HistoMapMap_t histoMapMap_mat_stepEnergyDeposit;
+    HistoMapMap_t histoMapMap_mat_stepEnergyNonIonDeposit;
+    HistoMapMap_t histoMapMap_mat_stepSecondaryKinetic;
+    HistoMapMap_t histoMapMap_mat_numberOfSteps;
+    HistoMapMap_t histoMapMap_mat_numberOfStepsPerInitialE;
+    HistoMapMap_t histoMapMap_mat_trackLengthPerInitialE;
+    HistoMapMap_t histoMapMap_mat_InitialE;
+
+    // distributions per process per particle type
+    HistoMapMap_t histoMapMap_prc_stepSize;
+    HistoMapMap_t histoMapMap_prc_stepKineticEnergy;
+    HistoMapMap_t histoMapMap_prc_postStepKineticEnergy;
+    HistoMapMap_t histoMapMap_prc_stepPseudorapidity;
+    HistoMapMap_t histoMapMap_prc_stepEnergyDeposit;
+    HistoMapMap_t histoMapMap_prc_stepEnergyNonIonDeposit;
+    HistoMapMap_t histoMapMap_prc_stepSecondaryKinetic;
+    HistoMapMap_t histoMapMap_prc_numberOfSteps;
+    HistoMapMap_t histoMapMap_prc_numberOfStepsPerInitialE;
+    HistoMapMap_t histoMapMap_prc_trackLengthPerInitialE;
+    HistoMapMap_t histoMapMap_prc_InitialE;
+
+    // all atlas
+    HistoMapMap_t histoMapMap_numberOfSteps;
+    HistoMapMap_t histoMapMap_numberOfStepsPerInitialE;
+    HistoMapMap_t histoMapMap_trackLengthPerInitialE;
+    HistoMapMap_t histoMapMap_InitialE;
+    HistoMapMap_t histoMapMap_stepKinetic;
+    HistoMapMap_t histoMapMap_postStepKinetic;
+
+    // 2D maps
+    HistoMapMap_t histoMapMap2D_vol_RZ;
+    HistoMapMap_t histoMapMap2D_mat_RZ;
+    HistoMapMap_t histoMapMap2D_prc_RZ;
+    HistoMapMap_t histoMapMap2D_vol_RZ_E;
+    HistoMapMap_t histoMapMap2D_mat_RZ_E;
+
+    // rather complicated function that merges two maps
+    void mergeMaps(HistoMapMap_t &selfMap, const HistoMapMap_t& refMap);
+
+    // function needed by ActionToolBaseReport base class
+    void merge(const Report & rep);
+  };
+  
+  /// ctor
+  //StepHistogram(const Config&);
+
+  const Report& getReport() const { return m_report; }
+
+
 
 private:
+  //Pointer to the run action
+  MyRunAction* m_run;
   MyEventAction*   fEventAction;
+  // report
+  Report m_report;
+
+  // initialize and fill histogram in a map
+  void InitializeFillHistogram2D(HistoMapMap_t &hMapMap, const char* suffix,
+                                              const G4String& pdgId, const G4String& vol,
+                                              int nbinsx, double xmin, double xmax,
+                                              int nbinsy, double ymin, double ymax,
+                                              double valuex, double valuey, double weight);
+
+  void InitializeFillHistogram(HistoMapMap_t &hMapMap, const char* suffix,
+                               const G4String& pdgId, const G4String& vol,
+                               int nbins, double xmin, double xmax, double value, double weight);
+
+  void InitializeFillHistogram(HistoMapMap_t &hMapMap, const char* suffix,
+                               const G4String& pdgId, const G4String& vol,
+                               int nbins, double *edges, double value, double weight);
+
+  float m_initialKineticEnergyOfStep;
+  G4String m_initialVolume;
+  G4String m_initialMaterial;
+  G4String m_initialProcess;
+  int m_trackID;
+
 };
 
 #endif

FullSimLight/src/MyActionInitialization.cc

@@ -61,7 +61,7 @@ void MyActionInitialization::Build() const {
           MyEventAction* evtact = new MyEventAction();
           SetUserAction(evtact);
           SetUserAction(new MyTrackingAction(evtact));
-          SetUserAction(new MySteppingAction(evtact));
+          SetUserAction(new MySteppingAction(evtact, runact));
           
       }
       else

FullSimLight/src/MyRunAction.cc

@@ -23,7 +23,7 @@ MyRunAction::MyRunAction(bool isGeantino, G4String geantinoMapFilename) : G4User
 }
 
 MyRunAction::~MyRunAction() {
-    if(fIsGeantino)
+    //if(fIsGeantino)
         delete G4AnalysisManager::Instance();
 }
 
@@ -97,11 +97,40 @@ void MyRunAction::BeginOfRunAction(const G4Run* /*aRun*/){
             analysisManager->SetP1YAxisTitle(fEtaInt_id,"Interaction Length #lambda");
         }
         
+    }
+    //Create step histograms only if not geantino simulation
+    else
+    {
+        //here create the "DEFAULT ONES"
+        // Create analysis manager
+        G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
+        if (isMaster) {
+        fMasterAnalysisManager = analysisManager;
+        G4cout<<"MyRunAction::BeginOfRunAction, created MASTER istance of the G4AnalysisManager: "<<fMasterAnalysisManager<<G4endl;
+                   
+        }
+            
+        G4cout<<"MyRunAction::BeginOfRunAction, created WORKER istance of the G4AnalysisManager: "<<analysisManager<<G4endl;
+            
+        G4cout << "Using G4AnalysisManager type: " << analysisManager->GetType() << G4endl;
+        analysisManager->SetVerboseLevel(1);
+           
+        G4String histoFileName = "Histo.root"; //TO DO : customizable
+        // Open output root file
+        G4cout<<"\n\nBeginOfRunAction: \n...create a root file using the G4AnalysisManager" << G4endl;
+        if (!analysisManager->OpenFile(histoFileName)){
+                G4cout<<"\nBeginOfRunAction ERROR: File cannot be opened!"<<G4endl;
+                exit(-1);
+            } else
+                G4cout<<"\n...output File "<<histoFileName<<" opened!"<<G4endl;
+            
+        ///There are no global/default histograms to create
+        
     }
     
     if (isMaster) {
         
-        //G4cout<<"\nBeginOfRunAction isMaster, and fMasterAnalysisManager: "<<fMasterAnalysisManager<<G4endl;
+        G4cout<<"\nBeginOfRunAction isMaster, and fMasterAnalysisManager: "<<fMasterAnalysisManager<<G4endl;
         
         std::vector<G4Region*>* regionVect =  G4RegionStore::GetInstance();
         int numRegions = regionVect->size();
@@ -141,7 +170,7 @@ void MyRunAction::BeginOfRunAction(const G4Run* /*aRun*/){
         fTimer = new G4Timer();
         fTimer->Start();
     }
-    //else G4cout<<"BeginOfRunAction isWorker, and fMasterAnalysisManager: "<<fMasterAnalysisManager<<G4endl;
+    else G4cout<<"BeginOfRunAction isWorker, and fMasterAnalysisManager: "<<fMasterAnalysisManager<<G4endl;
 }
 
 
@@ -159,6 +188,20 @@ void MyRunAction::EndOfRunAction(const G4Run*) {
             G4cout<<"Output File successfully saved and closed! " << G4endl;
         }
         
+    }
+    //Histogram stepping action TO DO: remove duplication of code?
+    else {
+        
+        auto analysisManager= G4AnalysisManager::Instance();
+        //Finalize analysisManager and Write out file
+        if (analysisManager->IsOpenFile()){
+            G4cout<<"\n\n EndOfRunAction, writing output file: "<<analysisManager->GetFileName() << G4endl;
+            analysisManager->Write();
+            G4cout<<"... closing histogramming step file: "<<analysisManager->GetFileName() << G4endl;
+            analysisManager->CloseFile();
+            G4cout<<"Output File successfully saved and closed! " << G4endl;
+        }
+        
     }
     if (isMaster) {
         fTimer->Stop();

FullSimLight/src/MySteppingAction.cc

@@ -7,17 +7,417 @@
 #include "G4ParticleDefinition.hh"
 
 
-MySteppingAction::MySteppingAction(MyEventAction* evtact) : G4UserSteppingAction(), fEventAction(evtact) {}
+//#include "G4StepPoint.hh"
+#include "G4VProcess.hh"
+#include "G4DebuggingHelper.hh"
+
+
+MySteppingAction::MySteppingAction(MyEventAction* evtact, MyRunAction* run) : G4UserSteppingAction(),
+m_run(run),
+fEventAction(evtact),
+m_initialKineticEnergyOfStep(),
+m_initialVolume(),
+m_initialMaterial(),
+m_initialProcess()
+{}
 
 
 MySteppingAction::~MySteppingAction() {}
 
 
 void MySteppingAction::UserSteppingAction(const G4Step* theStep) {
+  std::cout<<"************"<<std::endl;
+  std::cout<<"User stepping action: START"<<std::endl;
   G4double edep      = theStep->GetTotalEnergyDeposit();
   G4double stepl     = theStep->GetStepLength();
   G4Track* theTrack  = theStep->GetTrack();
   G4bool   isCharged = (theTrack->GetDefinition()->GetPDGCharge()!=0.);
   G4int primTrackID  = static_cast<MyTrackInformation*>(theTrack->GetUserInformation())->GetPrimaryTrackID();
   fEventAction->AddData(edep, stepl, isCharged, primTrackID-1);
+    
+
+  // track
+  G4Track *tr = theStep->GetTrack();
+
+  // pre-step point
+  G4StepPoint *PreStepPoint = theStep->GetPreStepPoint();
+
+  // post-step point
+  G4StepPoint *PostStepPoint = theStep->GetPostStepPoint();
+
+  // pre-step kinetic energy
+  double stepKinetic = PreStepPoint->GetKineticEnergy();
+
+  // post-step kinetic energy
+  double postStepKinetic = PostStepPoint->GetKineticEnergy();
+
+  // pre-step position
+  const G4ThreeVector& myPos = PreStepPoint->GetPosition();
+
+  // particle name
+  G4String particleName = "nucleus";
+  
+  if (!(tr->GetDefinition()->GetParticleType() == "nucleus"))
+    particleName = G4DebuggingHelpers::ClassifyParticle(tr->GetParticleDefinition());
+
+  // pre-step volume
+  G4String volumeName = PreStepPoint->GetPhysicalVolume()->GetName();
+  volumeName= G4DebuggingHelpers::ClassifyVolume(volumeName);
+
+
+  // pre-step material
+  G4String materialName = PreStepPoint->GetMaterial()->GetName();
+  materialName = G4DebuggingHelpers::ClassifyMaterial(materialName);
+
+  // process name (uses post-step point)
+  G4String processName = PostStepPoint->GetProcessDefinedStep() ?
+                         PostStepPoint->GetProcessDefinedStep()->GetProcessName() : "Unknown";
+
+  // secondaries
+  const std::vector<const G4Track*>* secondaries = theStep->GetSecondaryInCurrentStep();
+  
+    
+  std::cout<<"particleName: "<<particleName<<std::endl;
+  std::cout<<"volumeName: "<<volumeName<<std::endl;
+  std::cout<<"materialName: "<<materialName<<std::endl;
+  std::cout<<"processName: "<<processName<<std::endl;
+
+  // 2D map
+  bool do2DHistograms = false; //TO DO - configurable member
+  if (do2DHistograms) {
+    float DepositedE = theStep->GetTotalEnergyDeposit();
+    InitializeFillHistogram2D(m_report.histoMapMap2D_vol_RZ, "vol_RZ", particleName, volumeName,
+                            2000, -20000, 20000, 1000, 0, 5000, myPos.getZ(), myPos.perp(), 1.);
+    InitializeFillHistogram2D(m_report.histoMapMap2D_mat_RZ, "mat_RZ", particleName, materialName,
+                            2000, -20000, 20000, 1000, 0, 5000, myPos.getZ(), myPos.perp(), 1.);
+    InitializeFillHistogram2D(m_report.histoMapMap2D_prc_RZ, "prc_RZ", particleName, processName,
+                            2000, -20000, 20000, 1000, 0, 5000, myPos.getZ(), myPos.perp(), 1.);
+    InitializeFillHistogram2D(m_report.histoMapMap2D_vol_RZ_E, "vol_RZ_E", particleName, volumeName,
+                            2000, -20000, 20000, 1000, 0, 5000, myPos.getZ(), myPos.perp(), DepositedE);
+    InitializeFillHistogram2D(m_report.histoMapMap2D_mat_RZ_E, "mat_RZ_E", particleName, materialName,
+                            2000, -20000, 20000, 1000, 0, 5000, myPos.getZ(), myPos.perp(), DepositedE);
+  }
+
+  // step length
+  InitializeFillHistogram(m_report.histoMapMap_vol_stepSize, "vol_stepLength", particleName, volumeName,
+                          1000, -12, 4, log10(theStep->GetStepLength()), 1.);
+  InitializeFillHistogram(m_report.histoMapMap_mat_stepSize, "mat_stepLength", particleName, materialName,
+                          1000, -12, 4, log10(theStep->GetStepLength()), 1.);
+  InitializeFillHistogram(m_report.histoMapMap_prc_stepSize, "prc_stepLength", particleName, processName,
+                          1000, -12, 4, log10(theStep->GetStepLength()), 1.);
+
+ 
+  // step pseudorapidity
+  InitializeFillHistogram(m_report.histoMapMap_vol_stepPseudorapidity, "vol_stepPseudorapidity", particleName, volumeName,
+                          200, -10, 10, myPos.eta(), 1.);
+  InitializeFillHistogram(m_report.histoMapMap_mat_stepPseudorapidity, "mat_stepPseudorapidity", particleName, materialName,
+                          200, -10, 10, myPos.eta(), 1.);
+  InitializeFillHistogram(m_report.histoMapMap_prc_stepPseudorapidity, "prc_stepPseudorapidity", particleName, processName,
+                          200, -10, 10, myPos.eta(), 1.);
+
+  // step kinetic energy
+  InitializeFillHistogram(m_report.histoMapMap_vol_stepKineticEnergy, "vol_stepKineticEnergy", particleName, volumeName,
+                          1000, -9, 7, log10(stepKinetic), 1.);
+  InitializeFillHistogram(m_report.histoMapMap_mat_stepKineticEnergy, "mat_stepKineticEnergy", particleName, materialName,
+                          1000, -9, 7, log10(stepKinetic), 1.);
+  InitializeFillHistogram(m_report.histoMapMap_prc_stepKineticEnergy, "prc_stepKineticEnergy", particleName, processName,
+                          1000, -9, 7, log10(stepKinetic), 1.);
+  InitializeFillHistogram(m_report.histoMapMap_stepKinetic, "stepKineticEnergy", particleName, "AllATLAS",
+                          1000, -9, 7, log10(stepKinetic), 1.);
+
+  // post step kinetic energy
+  InitializeFillHistogram(m_report.histoMapMap_vol_postStepKineticEnergy, "vol_postStepKineticEnergy", particleName, volumeName,
+                          1000, -9, 7, log10(postStepKinetic), 1.);
+  InitializeFillHistogram(m_report.histoMapMap_mat_postStepKineticEnergy, "mat_postStepKineticEnergy", particleName, materialName,
+                          1000, -9, 7, log10(postStepKinetic), 1.);
+  InitializeFillHistogram(m_report.histoMapMap_prc_postStepKineticEnergy, "prc_postStepKineticEnergy", particleName, processName,
+                          1000, -9, 7, log10(postStepKinetic), 1.);
+  InitializeFillHistogram(m_report.histoMapMap_postStepKinetic, "postStepKineticEnergy", particleName, "AllATLAS",
+                          1000, -9, 7, log10(postStepKinetic), 1.);
+
+  // step energy deposit
+  InitializeFillHistogram(m_report.histoMapMap_vol_stepEnergyDeposit, "vol_stepEnergyDeposit", particleName, volumeName,
+                          1000, -11, 3, log10(theStep->GetTotalEnergyDeposit()), 1.);
+  InitializeFillHistogram(m_report.histoMapMap_mat_stepEnergyDeposit, "mat_stepEnergyDeposit", particleName, materialName,
+                          1000, -11, 3, log10(theStep->GetTotalEnergyDeposit()), 1.);
+  InitializeFillHistogram(m_report.histoMapMap_prc_stepEnergyDeposit, "prc_stepEnergyDeposit", particleName, processName,
+                          1000, -11, 3, log10(theStep->GetTotalEnergyDeposit()), 1.);
+
+  // step non-ionizing energy deposit
+  InitializeFillHistogram(m_report.histoMapMap_vol_stepEnergyNonIonDeposit, "vol_stepEnergyNonIonDeposit", particleName, volumeName,
+                          1000, -11, 1, log10(theStep->GetNonIonizingEnergyDeposit()), 1.);
+  InitializeFillHistogram(m_report.histoMapMap_mat_stepEnergyNonIonDeposit, "mat_stepEnergyNonIonDeposit", particleName, materialName,
+                          1000, -11, 1, log10(theStep->GetNonIonizingEnergyDeposit()), 1.);
+  InitializeFillHistogram(m_report.histoMapMap_prc_stepEnergyNonIonDeposit, "prc_stepEnergyNonIonDeposit", particleName, processName,
+                          1000, -11, 1, log10(theStep->GetNonIonizingEnergyDeposit()), 1.);
+
+  // secondary kinetic energy
+  for (const auto &track : *secondaries) {
+    G4String secondary_particleName = G4DebuggingHelpers::ClassifyParticle(track->GetParticleDefinition());
+    InitializeFillHistogram(m_report.histoMapMap_vol_stepSecondaryKinetic, "vol_stepSecondaryKinetic", secondary_particleName, volumeName,
+                            1000, -7, 5, log10(track->GetKineticEnergy()), 1.);
+    InitializeFillHistogram(m_report.histoMapMap_mat_stepSecondaryKinetic, "mat_stepSecondaryKinetic", secondary_particleName, materialName,
+                            1000, -7, 5, log10(track->GetKineticEnergy()), 1.);
+    InitializeFillHistogram(m_report.histoMapMap_prc_stepSecondaryKinetic, "prc_stepSecondaryKinetic", secondary_particleName, processName,
+                            1000, -7, 5, log10(track->GetKineticEnergy()), 1.);
+  }
+
+  bool doGeneralHistograms = false; //TO DO: configurable
+  // stop here if 'general' histograms not activated
+  // _______________________________________________
+  if (!doGeneralHistograms)
+    return;
+    
+  // first step (after initial step)
+    if (tr->GetCurrentStepNumber()==1) {
+      // initial kinetic energy
+      m_initialKineticEnergyOfStep = stepKinetic;
+      
+      // initial volume/material/processes
+      m_initialVolume = volumeName;
+      m_initialMaterial = materialName;
+      m_initialProcess = processName;
+
+      // save track ID for checking if we later have the same track
+      m_trackID = tr->GetTrackID();
+
+      // initial energy
+      InitializeFillHistogram(m_report.histoMapMap_InitialE, "InitialE", particleName, "AllATLAS",
+                              1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), 1.0);
+      InitializeFillHistogram(m_report.histoMapMap_vol_InitialE, "vol_InitialE", particleName, m_initialVolume,
+                              1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), 1.0);
+      InitializeFillHistogram(m_report.histoMapMap_mat_InitialE, "mat_InitialE", particleName, m_initialMaterial,
+                              1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), 1.0);
+      InitializeFillHistogram(m_report.histoMapMap_prc_InitialE, "prc_InitialE", particleName, m_initialProcess,
+                              1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), 1.0);
+    }
+
+    // last step
+    if ( tr->GetTrackStatus() == 2 ) {
+      // assert to check if we have the correct track
+      if (not (tr->GetTrackID() == m_trackID)) {
+        std::cerr<<"Track ID changed between the assumed first step and the last."<<std::endl;
+        throw std::exception();
+      }
+      // number of steps
+      int nSteps = tr->GetCurrentStepNumber() + 1;
+      InitializeFillHistogram(m_report.histoMapMap_numberOfSteps, "numberOfSteps", particleName, "AllATLAS",
+                              10000, 0.5, 10000.5, nSteps, 1.);
+      InitializeFillHistogram(m_report.histoMapMap_vol_numberOfSteps, "vol_numberOfSteps", particleName, m_initialVolume,
+                              10000, 0.5, 10000.5, nSteps, 1.);
+      InitializeFillHistogram(m_report.histoMapMap_mat_numberOfSteps, "mat_numberOfSteps", particleName, m_initialMaterial,
+                              10000, 0.5, 10000.5, nSteps, 1.);
+      InitializeFillHistogram(m_report.histoMapMap_prc_numberOfSteps, "prc_numberOfSteps", particleName, m_initialProcess,
+                              10000, 0.5, 10000.5, nSteps, 1.);
+      // number of steps vs initial energy
+      InitializeFillHistogram(m_report.histoMapMap_numberOfStepsPerInitialE, "numberOfStepsPerInitialE", particleName, "AllATLAS",
+                              1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), nSteps);
+      InitializeFillHistogram(m_report.histoMapMap_vol_numberOfStepsPerInitialE, "vol_numberOfStepsPerInitialE", particleName, m_initialVolume,
+                              1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), nSteps);
+      InitializeFillHistogram(m_report.histoMapMap_mat_numberOfStepsPerInitialE, "mat_numberOfStepsPerInitialE", particleName, m_initialMaterial,
+                              1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), nSteps);
+      InitializeFillHistogram(m_report.histoMapMap_prc_numberOfStepsPerInitialE, "prc_numberOfStepsPerInitialE", particleName, m_initialProcess,
+                              1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), nSteps);
+      // track length vs initial energy
+      InitializeFillHistogram(m_report.histoMapMap_trackLengthPerInitialE, "trackLengthPerInitialE", particleName, "AllATLAS",
+                              1000, -9, 7, log10(tr->GetTrackLength()), 1.);
+      InitializeFillHistogram(m_report.histoMapMap_vol_trackLengthPerInitialE, "vol_trackLengthPerInitialE", particleName, m_initialVolume,
+                              1000, -9, 7, log10(tr->GetTrackLength()), 1.);
+      InitializeFillHistogram(m_report.histoMapMap_mat_trackLengthPerInitialE, "mat_trackLengthPerInitialE", particleName, m_initialMaterial,
+                              1000, -9, 7, log10(tr->GetTrackLength()), 1.);
+      InitializeFillHistogram(m_report.histoMapMap_prc_trackLengthPerInitialE, "prc_trackLengthPerInitialE", particleName, m_initialProcess,
+                              1000, -9, 7, log10(tr->GetTrackLength()), 1.);
+    }
+    std::cout<<"User Stepping action DONE"<<std::endl;
+    std::cout<<"************"<<std::endl;
+}
+
+void MySteppingAction::InitializeFillHistogram2D(HistoMapMap_t &hMapMap, const char* suffix,
+                                              const G4String& particleName, const G4String& vol,
+                                              int nbinsx, double xmin, double xmax,
+                                              int nbinsy, double ymin, double ymax,
+                                              double valuex, double valuey, double weight)
+{
+  if ( hMapMap.find(vol) == hMapMap.end() ) {
+    // initialize HistoMap_t if not yet exist
+    hMapMap.emplace(vol,HistoMap_t());
+  }
+  HistoMap_t &hMap = hMapMap[vol];
+    
+  auto analysisManager = G4AnalysisManager::Instance();
+  
+  if ( hMap.find(particleName) == hMap.end() ) {
+    // initialize histogram if not yet exist
+    std::ostringstream stringStream;
+    stringStream << vol << "_" << particleName << "_" << suffix;
+    //HERE I create the histogram if it doesn't exist
+    //DONE in what follows
+    //hMap[particleName] = new TH2F(stringStream.str().c_str(), stringStream.str().c_str(), nbinsx, xmin, xmax, nbinsy, ymin, ymax);
+    int id = m_run->fMasterAnalysisManager->GetH2Id(stringStream.str().c_str(), false);
+    if(id<0) {
+          //G4cout<<"Geant4: Eta rad profile of "<<pathEtaRL<<" didn't exist, creating P1 now on MASTER:  "<<m_run->fMasterAnalysisManager<<G4endl;
+          //const std::string name(detName+"_RL");
+          //hMap[particleName] =
+          hMap[particleName] = m_run->fMasterAnalysisManager->CreateH2(stringStream.str().c_str(), stringStream.str().c_str(), nbinsx, xmin, xmax, nbinsy, ymin, ymax);
+//          m_run->fMasterAnalysisManager->SetH2XAxisTitle(hMap[particleName], "#eta");
+//          m_run->fMasterAnalysisManager->SetP1YAxisTitle(hMap[particleName], "%X0");
+          
+      }
+      id=analysisManager->GetH2Id(stringStream.str().c_str(), false);
+      if(id<0)
+      {
+          //G4cout<<"Geant4: Eta rad profile of "<<pathEtaRL<<" didn't exist, creating P1 now on WORKER: "<<analysisManager<<G4endl;
+          hMap[particleName] = analysisManager->CreateH2(stringStream.str().c_str(), stringStream.str().c_str(), nbinsx, xmin, xmax, nbinsy, ymin, ymax);
+//          analysisManager->SetP1XAxisTitle(hMap[particleName], "#eta");
+//          analysisManager->SetP1YAxisTitle(hMap[particleName], "%X0");
+          
+      }
+    
+  }
+  //TO DO ((TH2*)hMap[particleName])->Fill(valuex, valuey, weight);
+  analysisManager->FillH2(hMap[particleName], valuex, valuey, weight);
+    
+}
+
+void MySteppingAction::InitializeFillHistogram(HistoMapMap_t &hMapMap, const char* suffix,
+                                            const G4String& particleName, const G4String& vol,
+                                            int nbins, double xmin, double xmax, double value, double weight)
+{
+  if ( hMapMap.find(vol) == hMapMap.end() ) {
+    // initialize HistoMap_t if not yet exist
+    hMapMap.emplace(vol,HistoMap_t());
+  }
+  HistoMap_t &hMap = hMapMap[vol];
+  auto analysisManager = G4AnalysisManager::Instance();
+  if ( hMap.find(particleName) == hMap.end() ) {
+    // initialize histogram if not yet exist
+    std::ostringstream stringStream;
+    stringStream << vol << "_" << particleName << "_" << suffix;
+    
+    //hMap[particleName] = new TH1F(stringStream.str().c_str(), stringStream.str().c_str(), nbins, xmin, xmax);
+    
+    int id = m_run->fMasterAnalysisManager->GetH1Id(stringStream.str().c_str(), false);
+    if(id<0) {
+        hMap[particleName] = m_run->fMasterAnalysisManager->CreateH1(stringStream.str().c_str(), stringStream.str().c_str(), nbins, xmin, xmax);
+        
+    }
+    id=analysisManager->GetH1Id(stringStream.str().c_str(), false);
+    if(id<0)
+    {
+        hMap[particleName] = analysisManager->CreateH1(stringStream.str().c_str(), stringStream.str().c_str(), nbins, xmin, xmax);
+    }
+    
+  }
+  //TO DO hMap[particleName]->Fill(value, weight);
+  analysisManager->FillH1(hMap[particleName], value, weight);
+}
+
+void MySteppingAction::InitializeFillHistogram(HistoMapMap_t &hMapMap, const char* suffix,
+                                            const G4String& particleName, const G4String& vol,
+                                            int nbins, double *edges, double value, double weight)
+{
+  if ( hMapMap.find(vol) == hMapMap.end() ) {
+    // initialize HistoMap_t if not yet exist
+    hMapMap.emplace(vol,HistoMap_t());
+  }
+  HistoMap_t &hMap = hMapMap[vol];
+  auto analysisManager = G4AnalysisManager::Instance();
+  if ( hMap.find(particleName) == hMap.end() ) {
+    // initialize histogram if not yet exist
+    std::ostringstream stringStream;
+    stringStream << vol << "_" << particleName << "_" << suffix;
+      
+    //hMap[particleName] = new TH1F(stringStream.str().c_str(), stringStream.str().c_str(), nbins, edges);
+    int id = m_run->fMasterAnalysisManager->GetH1Id(stringStream.str().c_str(), false);
+    if(id<0) {
+        hMap[particleName] = m_run->fMasterAnalysisManager->CreateH1(stringStream.str().c_str(), stringStream.str().c_str(), nbins, edges[0], edges[nbins]); //TO DO : check if this is proper and actually used
+        
+    }
+    id=analysisManager->GetH2Id(stringStream.str().c_str(), false);
+    if(id<0)
+    {
+        hMap[particleName] = analysisManager->CreateH1(stringStream.str().c_str(), stringStream.str().c_str(), nbins, edges[0], edges[nbins]);
+    }
+    
+  }
+  //hMap[particleName]->Fill(value, weight);
+  analysisManager->FillH1(hMap[particleName], value, weight);
+  
 }
+
+void MySteppingAction::Report::mergeMaps(HistoMapMap_t &selfMap, const HistoMapMap_t& refMap) {
+  for (auto const& ref : refMap)
+  {
+    if ( selfMap.find(ref.first) == selfMap.end() ) {
+      // HistoMap_t does not yet exist
+      selfMap.emplace(ref.first, ref.second);
+    }
+    else {
+      HistoMap_t &target = selfMap[ref.first];
+      for (auto const& hm : ref.second)
+      {
+        if ( target.find(hm.first) == target.end() ) {
+          // histogram does not yet exist
+          target.emplace(hm.first, hm.second);
+        }
+        else {
+          // add histograms
+          //TO DO target[hm.first]->Add(hm.second);
+        }
+      }
+    }
+  }
+}
+
+void MySteppingAction::Report::merge(const Report & rep) {
+  mergeMaps(histoMapMap_vol_stepSize, rep.histoMapMap_vol_stepSize);
+  mergeMaps(histoMapMap_vol_stepKineticEnergy, rep.histoMapMap_vol_stepKineticEnergy);
+  mergeMaps(histoMapMap_vol_postStepKineticEnergy, rep.histoMapMap_vol_postStepKineticEnergy);
+  mergeMaps(histoMapMap_vol_stepPseudorapidity, rep.histoMapMap_vol_stepPseudorapidity);
+  mergeMaps(histoMapMap_vol_stepEnergyDeposit, rep.histoMapMap_vol_stepEnergyDeposit);
+  mergeMaps(histoMapMap_vol_stepEnergyNonIonDeposit, rep.histoMapMap_vol_stepEnergyNonIonDeposit);
+  mergeMaps(histoMapMap_vol_stepSecondaryKinetic, rep.histoMapMap_vol_stepSecondaryKinetic);
+
+  mergeMaps(histoMapMap_mat_stepSize, rep.histoMapMap_mat_stepSize);
+  mergeMaps(histoMapMap_mat_stepKineticEnergy, rep.histoMapMap_mat_stepKineticEnergy);
+  mergeMaps(histoMapMap_mat_postStepKineticEnergy, rep.histoMapMap_mat_postStepKineticEnergy);
+  mergeMaps(histoMapMap_mat_stepPseudorapidity, rep.histoMapMap_mat_stepPseudorapidity);
+  mergeMaps(histoMapMap_mat_stepEnergyDeposit, rep.histoMapMap_mat_stepEnergyDeposit);
+  mergeMaps(histoMapMap_mat_stepEnergyNonIonDeposit, rep.histoMapMap_mat_stepEnergyNonIonDeposit);
+  mergeMaps(histoMapMap_mat_stepSecondaryKinetic, rep.histoMapMap_mat_stepSecondaryKinetic);
+
+  mergeMaps(histoMapMap_prc_stepSize, rep.histoMapMap_prc_stepSize);
+  mergeMaps(histoMapMap_prc_stepKineticEnergy, rep.histoMapMap_prc_stepKineticEnergy);
+  mergeMaps(histoMapMap_prc_postStepKineticEnergy, rep.histoMapMap_prc_postStepKineticEnergy);
+  mergeMaps(histoMapMap_prc_stepPseudorapidity, rep.histoMapMap_prc_stepPseudorapidity);
+  mergeMaps(histoMapMap_prc_stepEnergyDeposit, rep.histoMapMap_prc_stepEnergyDeposit);
+  mergeMaps(histoMapMap_prc_stepEnergyNonIonDeposit, rep.histoMapMap_prc_stepEnergyNonIonDeposit);
+  mergeMaps(histoMapMap_prc_stepSecondaryKinetic, rep.histoMapMap_prc_stepSecondaryKinetic);
+
+  mergeMaps(histoMapMap_numberOfSteps, rep.histoMapMap_numberOfSteps);
+  mergeMaps(histoMapMap_vol_numberOfSteps, rep.histoMapMap_vol_numberOfSteps);
+  mergeMaps(histoMapMap_mat_numberOfSteps, rep.histoMapMap_mat_numberOfSteps);
+  mergeMaps(histoMapMap_prc_numberOfSteps, rep.histoMapMap_prc_numberOfSteps);
+  mergeMaps(histoMapMap_numberOfStepsPerInitialE, rep.histoMapMap_numberOfStepsPerInitialE);
+  mergeMaps(histoMapMap_vol_numberOfStepsPerInitialE, rep.histoMapMap_vol_numberOfStepsPerInitialE);
+  mergeMaps(histoMapMap_mat_numberOfStepsPerInitialE, rep.histoMapMap_mat_numberOfStepsPerInitialE);
+  mergeMaps(histoMapMap_prc_numberOfStepsPerInitialE, rep.histoMapMap_prc_numberOfStepsPerInitialE);
+  mergeMaps(histoMapMap_trackLengthPerInitialE, rep.histoMapMap_trackLengthPerInitialE);
+  mergeMaps(histoMapMap_vol_trackLengthPerInitialE, rep.histoMapMap_vol_trackLengthPerInitialE);
+  mergeMaps(histoMapMap_mat_trackLengthPerInitialE, rep.histoMapMap_mat_trackLengthPerInitialE);
+  mergeMaps(histoMapMap_prc_trackLengthPerInitialE, rep.histoMapMap_prc_trackLengthPerInitialE);
+  mergeMaps(histoMapMap_InitialE, rep.histoMapMap_InitialE);
+  mergeMaps(histoMapMap_vol_InitialE, rep.histoMapMap_vol_InitialE);
+  mergeMaps(histoMapMap_mat_InitialE, rep.histoMapMap_mat_InitialE);
+  mergeMaps(histoMapMap_prc_InitialE, rep.histoMapMap_prc_InitialE);
+  mergeMaps(histoMapMap_stepKinetic, rep.histoMapMap_stepKinetic);
+  mergeMaps(histoMapMap_postStepKinetic, rep.histoMapMap_postStepKinetic);
+
+  mergeMaps(histoMapMap2D_vol_RZ, rep.histoMapMap2D_vol_RZ);
+  mergeMaps(histoMapMap2D_mat_RZ, rep.histoMapMap2D_mat_RZ);
+  mergeMaps(histoMapMap2D_prc_RZ, rep.histoMapMap2D_prc_RZ);
+  mergeMaps(histoMapMap2D_vol_RZ_E, rep.histoMapMap2D_vol_RZ_E);
+  mergeMaps(histoMapMap2D_mat_RZ_E, rep.histoMapMap2D_mat_RZ_E);
+}
+
+