From 90975a7c24076edd417f1801198399a83f94be9c Mon Sep 17 00:00:00 2001
From: Jan Kieseler <jkiesele@cern.ch>
Date: Fri, 8 Sep 2023 14:33:33 +0200
Subject: [PATCH] snapshot with output

---
 exampleB4a.cc                     | 19 ++++---
 include/ConstructionWrapper.hh    | 57 ++++++++++++++++++-
 include/DetectorConstruction.hh   | 31 ++++-------
 include/EventAction.hh            | 33 ++++++-----
 include/PrimaryGeneratorAction.hh | 10 +++-
 include/RunAction.hh              | 14 +++++
 include/SteppingAction.hh         |  2 +
 src/ActionInitialization.cc       | 13 ++++-
 src/ConstructionWrapper.cc        | 67 ++++++++++++++++++++++
 src/DetectorConstruction.cc       | 41 ++++++++++++--
 src/EventAction.cc                | 92 ++++++++++++++++++++-----------
 src/PrimaryGeneratorAction.cc     | 26 +++++++--
 src/RunAction.cc                  | 59 +++++++-------------
 src/SteppingAction.cc             | 26 ++++++---
 14 files changed, 347 insertions(+), 143 deletions(-)

diff --git a/exampleB4a.cc b/exampleB4a.cc
index ec929cd..7a819f2 100644
--- a/exampleB4a.cc
+++ b/exampleB4a.cc
@@ -39,6 +39,8 @@
 #include "FTFP_BERT.hh"
 #include "Randomize.hh"
 
+#include "ConstructionWrapper.hh"
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 namespace {
@@ -112,18 +114,18 @@ int main(int argc,char** argv)
   }
 #endif
 
-  ConstructionWrapper cw;
-  cw.addLayer(1, "G4_Pb", false, 2, 2);
-  cw.addLayer(10, "G4_Si", true,  2, 2);
-  cw.addLayer(1, "G4_Pb", false, 2, 2);
-  cw.addLayer(20, "G4_Si", true,  2, 2);
-  cw.addLayer(4, "G4_Pb", false, 2, 2);
-  cw.addLayer(3, "G4_Si", true,  2, 2);
+  ConstructionWrapper * cw = new ConstructionWrapper();
+  cw->addLayer(1, "G4_Pb", false);
+  cw->addLayer(10, "G4_Si", true, 9);
+  cw->addLayer(1, "G4_Pb", false);
+  cw->addLayer(20, "G4_Si", true, 15);
+  cw->addLayer(4, "G4_Pb", false);
+  cw->addLayer(3, "G4_Si", true,  23);
+
 
   // Set mandatory initialization classes
   //
   auto detConstruction = new B4::DetectorConstruction(cw);
-  // FOR LATER: assign constructionwrapper here (will be passed from other function)
 
   runManager->SetUserInitialization(detConstruction);
 
@@ -167,6 +169,7 @@ int main(int argc,char** argv)
 
   delete visManager;
   delete runManager;
+  delete cw;
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
diff --git a/include/ConstructionWrapper.hh b/include/ConstructionWrapper.hh
index b93aae5..ba716ee 100644
--- a/include/ConstructionWrapper.hh
+++ b/include/ConstructionWrapper.hh
@@ -1,7 +1,34 @@
 
+#ifndef CONSTRUCTIONWRAPPER_HH
+#define CONSTRUCTIONWRAPPER_HH
+
+#include "G4ThreeVector.hh"
 #include <vector>
 #include <string>
-#include "G4VPhysicalVolume.hh"
+//#include "G4VPhysicalVolume.hh"
+class G4VPhysicalVolume;
+
+
+class Sensor{
+public:
+    Sensor(){};
+    ~Sensor(){};
+
+    const G4double getEnergy()const{
+        return energy;
+    }
+    G4ThreeVector getPos()const{
+        return position;
+    }
+    G4ThreeVector getSize()const{
+        return size;
+    }
+
+    
+    G4ThreeVector position;
+    G4ThreeVector size;
+    mutable G4double energy;
+};
 
 class Layer{
     
@@ -25,6 +52,8 @@ public:
     bool isActive;
     G4VPhysicalVolume* physicalVolume;
     std::string name;
+
+    std::vector<Sensor> sensors;
 };
 
 
@@ -34,7 +63,7 @@ public:
   ConstructionWrapper(double xy_width=50);
   ~ConstructionWrapper() {};
 
-  void addLayer(double thickness_cm, std::string material, bool isActive=true, int nx=1, int ny=1);
+  void addLayer(double thickness_cm, std::string material, bool isActive=true, int nx=1, int ny=-1);
   
    std::vector<Layer>& getLayers() ;
    const std::vector<Layer>& getLayers() const;
@@ -42,7 +71,29 @@ public:
    double getXYWidth() const{
        return xywidth;
    }
+    void resetSensorEnergies()const;//energies are mutable
+    int getNSensors()const{ 
+        int n_sensors = 0;
+        for(const auto& layer: layers){
+            n_sensors += layer.sensors.size();
+        }
+        return n_sensors;
+    }
+
+     Layer* getLayerByVolume(G4VPhysicalVolume* volume);
+
+     const Layer* getLayerByVolume(G4VPhysicalVolume* volume)const;
+
+     Sensor* getSensorByVolume(G4VPhysicalVolume* volume);
+
+     const Sensor* getSensorByVolume(G4VPhysicalVolume* volume)const;
+
+    void printSensorEnergies()const;
+
+
 private:
     double xywidth;
     std::vector<Layer> layers;
-};
\ No newline at end of file
+};
+
+#endif
\ No newline at end of file
diff --git a/include/DetectorConstruction.hh b/include/DetectorConstruction.hh
index 9db1e56..e32ae51 100644
--- a/include/DetectorConstruction.hh
+++ b/include/DetectorConstruction.hh
@@ -32,10 +32,10 @@
 
 #include "G4VUserDetectorConstruction.hh"
 #include "globals.hh"
-#include "ConstructionWrapper.hh"
 
 class G4VPhysicalVolume;
 class G4GlobalMagFieldMessenger;
+class ConstructionWrapper;
 
 namespace B4
 {
@@ -59,19 +59,20 @@ namespace B4
 class DetectorConstruction : public G4VUserDetectorConstruction
 {
   public:
-    DetectorConstruction(const ConstructionWrapper& cw){
-        this->cw = cw;
-    }
+    DetectorConstruction( ConstructionWrapper* Cw): G4VUserDetectorConstruction(), cw(Cw){}
     ~DetectorConstruction() override = default;
 
   public:
     G4VPhysicalVolume* Construct() override;
     void ConstructSDandField() override;
 
-    // get methods
-    //
-    const G4VPhysicalVolume* GetAbsorberPV() const;
-    const G4VPhysicalVolume* GetGapPV() const;
+
+    ConstructionWrapper* getConstructionWrapper(){
+        return cw;
+    }
+    const ConstructionWrapper* getConstructionWrapper()const{
+        return cw;
+    }
 
   private:
     // methods
@@ -84,23 +85,11 @@ class DetectorConstruction : public G4VUserDetectorConstruction
     static G4ThreadLocal G4GlobalMagFieldMessenger*  fMagFieldMessenger;
                                       // magnetic field messenger
 
-    G4VPhysicalVolume* fAbsorberPV = nullptr; // the absorber physical volume
-    G4VPhysicalVolume* fGapPV = nullptr;      // the gap physical volume
-
     G4bool fCheckOverlaps = true; // option to activate checking of volumes overlaps
 
-    ConstructionWrapper cw;
+    ConstructionWrapper* cw=nullptr;
 };
 
-// inline functions
-
-inline const G4VPhysicalVolume* DetectorConstruction::GetAbsorberPV() const {
-  return fAbsorberPV;
-}
-
-inline const G4VPhysicalVolume* DetectorConstruction::GetGapPV() const  {
-  return fGapPV;
-}
 
 }
 
diff --git a/include/EventAction.hh b/include/EventAction.hh
index 0f734ed..d8e9160 100644
--- a/include/EventAction.hh
+++ b/include/EventAction.hh
@@ -31,7 +31,10 @@
 #define B4aEventAction_h 1
 
 #include "G4UserEventAction.hh"
+#include "RunAction.hh"
+#include "PrimaryGeneratorAction.hh"
 #include "globals.hh"
+class  ConstructionWrapper;
 
 namespace B4a
 {
@@ -51,30 +54,26 @@ class EventAction : public G4UserEventAction
     ~EventAction() override = default;
 
     void  BeginOfEventAction(const G4Event* event) override;
-    void    EndOfEventAction(const G4Event* event) override;
+    void  EndOfEventAction(const G4Event* event) override;
 
-    void AddAbs(G4double de, G4double dl);
-    void AddGap(G4double de, G4double dl);
+    void setConstructionWrapper(const ConstructionWrapper * Cw){
+        this->cw = Cw;
+    }
+    void setPrimaryGeneratorAction(const B4::PrimaryGeneratorAction * Gen){
+        this->gen = Gen;
+    }
+    void setRunAction(const B4::RunAction * Run){
+        this->run = Run;
+    }
 
   private:
-    G4double  fEnergyAbs = 0.;
-    G4double  fEnergyGap = 0.;
-    G4double  fTrackLAbs = 0.;
-    G4double  fTrackLGap = 0.;
+    const ConstructionWrapper * cw=nullptr;
+    const B4::PrimaryGeneratorAction * gen=nullptr;
+    const B4::RunAction * run=nullptr;
 };
 
 // inline functions
 
-inline void EventAction::AddAbs(G4double de, G4double dl) {
-  fEnergyAbs += de;
-  fTrackLAbs += dl;
-}
-
-inline void EventAction::AddGap(G4double de, G4double dl) {
-  fEnergyGap += de;
-  fTrackLGap += dl;
-}
-
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/include/PrimaryGeneratorAction.hh b/include/PrimaryGeneratorAction.hh
index 7240a83..1b84a98 100644
--- a/include/PrimaryGeneratorAction.hh
+++ b/include/PrimaryGeneratorAction.hh
@@ -56,9 +56,17 @@ public:
 
   // set methods
   void SetRandomFlag(G4bool value);
-
+  G4double getPartEnergy() const{return partEnergy;}
+  void setMinPartEnergy(G4double value){minPartEnergy = value;}
+  void setMaxPartEnergy(G4double value){maxPartEnergy = value;}
+  void setPartSpecies(G4String value){partSpecies = value;}
 private:
   G4ParticleGun* fParticleGun = nullptr; // G4 particle gun
+  G4double partEnergy;
+  G4double minPartEnergy; //MeV
+  G4double maxPartEnergy;
+  G4String partSpecies;
+
 };
 
 }
diff --git a/include/RunAction.hh b/include/RunAction.hh
index 97cae21..7324984 100644
--- a/include/RunAction.hh
+++ b/include/RunAction.hh
@@ -64,6 +64,20 @@ class RunAction : public G4UserRunAction
 
     void BeginOfRunAction(const G4Run*) override;
     void   EndOfRunAction(const G4Run*) override;
+
+//now this is dumb but the way this is setup it has to be done like this
+
+//yeah mutable fine whatever
+  mutable std::vector<double> hitEnergy;
+  mutable std::vector<double> hitX;
+  mutable std::vector<double> hitY;
+  mutable std::vector<double> hitZ;
+  mutable std::vector<double> hitDX;
+  mutable std::vector<double> hitDY;
+  mutable std::vector<double> hitDZ;
+  mutable std::vector<int> hitLayer;
+  mutable std::vector<int> hitCopyNumber;
+
 };
 
 }
diff --git a/include/SteppingAction.hh b/include/SteppingAction.hh
index 8b2ddd1..40485a6 100644
--- a/include/SteppingAction.hh
+++ b/include/SteppingAction.hh
@@ -31,6 +31,7 @@
 #define B4aSteppingAction_h 1
 
 #include "G4UserSteppingAction.hh"
+#include "ConstructionWrapper.hh"
 
 namespace B4
 {
@@ -57,6 +58,7 @@ public:
 
   void UserSteppingAction(const G4Step* step) override;
 
+
 private:
   const B4::DetectorConstruction* fDetConstruction = nullptr;
   EventAction* fEventAction = nullptr;
diff --git a/src/ActionInitialization.cc b/src/ActionInitialization.cc
index 4f97a02..b71ba21 100644
--- a/src/ActionInitialization.cc
+++ b/src/ActionInitialization.cc
@@ -52,15 +52,22 @@ void ActionInitialization::BuildForMaster() const
   SetUserAction(new RunAction);
 }
 
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 void ActionInitialization::Build() const
 {
-  SetUserAction(new PrimaryGeneratorAction);
-  SetUserAction(new RunAction);
+  auto gen = new PrimaryGeneratorAction;
+  SetUserAction(gen);
+  auto runact = new RunAction;
+  SetUserAction(runact);
   auto eventAction = new EventAction;
+  eventAction->setConstructionWrapper(fDetConstruction->getConstructionWrapper());
+  eventAction->setRunAction(runact);
+  eventAction->setPrimaryGeneratorAction(gen);
   SetUserAction(eventAction);
-  SetUserAction(new SteppingAction(fDetConstruction,eventAction));
+  auto steppingAction = new SteppingAction(fDetConstruction,eventAction);
+  SetUserAction(steppingAction);
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/src/ConstructionWrapper.cc b/src/ConstructionWrapper.cc
index 4273c75..2505906 100644
--- a/src/ConstructionWrapper.cc
+++ b/src/ConstructionWrapper.cc
@@ -1,4 +1,7 @@
+
+
 #include "ConstructionWrapper.hh"
+#include "G4VPhysicalVolume.hh"
 
 Layer::Layer(){
     thickness = 0;
@@ -38,6 +41,10 @@ ConstructionWrapper::ConstructionWrapper(double xy_width){
 }
 
 void ConstructionWrapper::addLayer(double thickness, std::string material, bool isActive, int nx, int ny){
+    if(ny<0){
+        ny = nx;
+    }
+
     Layer layer;
     layer.setThickness(thickness);
     layer.setMaterial(material);
@@ -60,3 +67,63 @@ const std::vector<Layer> & ConstructionWrapper::getLayers() const{
 std::vector<Layer> & ConstructionWrapper::getLayers(){
     return layers;
 }
+
+void ConstructionWrapper::resetSensorEnergies()const {
+    for(const auto & layer : layers){
+        for(const auto & sensor : layer.sensors){
+            sensor.energy = 0;
+        }
+    }
+}
+
+Layer* ConstructionWrapper::getLayerByVolume(G4VPhysicalVolume* volume){
+    for(auto & layer : layers){
+        if(layer.physicalVolume == volume){
+            return &layer;
+        }
+    }
+    return nullptr;
+}
+
+const Layer* ConstructionWrapper::getLayerByVolume(G4VPhysicalVolume* volume)const{
+    for(auto & layer : layers){
+        if(layer.physicalVolume == volume){
+            return &layer;
+        }
+    }
+    return nullptr;
+}
+
+Sensor* ConstructionWrapper::getSensorByVolume(G4VPhysicalVolume* volume){
+    int copyNo = volume->GetCopyNo();
+    auto layer = getLayerByVolume(volume);
+    if(layer == nullptr){
+        return nullptr;
+    }
+    if(copyNo >= layer->sensors.size()){
+        return nullptr;
+    }
+    return &layer->sensors[copyNo];
+}
+
+const Sensor* ConstructionWrapper::getSensorByVolume(G4VPhysicalVolume* volume)const{
+    int copyNo = volume->GetCopyNo();
+    auto layer = getLayerByVolume(volume);
+    if(layer == nullptr){
+        return nullptr;
+    }
+    if(copyNo >= layer->sensors.size()){
+        return nullptr;
+    }
+    return &layer->sensors[copyNo];
+}
+
+
+void ConstructionWrapper::printSensorEnergies()const{
+    for(const auto& layer: layers){
+        for(const auto& sensor: layer.sensors){
+            G4cout << sensor.energy << " ";
+        }
+        G4cout << G4endl;
+    }
+}
\ No newline at end of file
diff --git a/src/DetectorConstruction.cc b/src/DetectorConstruction.cc
index 44c437d..56e8a3b 100644
--- a/src/DetectorConstruction.cc
+++ b/src/DetectorConstruction.cc
@@ -28,6 +28,7 @@
 /// \brief Implementation of the B4::DetectorConstruction class
 
 #include "DetectorConstruction.hh"
+#include "ConstructionWrapper.hh"
 
 #include "G4Material.hh"
 #include "G4NistManager.hh"
@@ -100,6 +101,7 @@ G4GlobalMagFieldMessenger* DetectorConstruction::fMagFieldMessenger = nullptr;
 
 G4VPhysicalVolume* DetectorConstruction::Construct()
 {
+  fCheckOverlaps=true;
   // Define materials
   DefineMaterials();
 
@@ -114,7 +116,7 @@ void DetectorConstruction::DefineMaterials()
   // Lead material defined using NIST Manager
   auto nistManager = G4NistManager::Instance();
 nistManager->FindOrBuildMaterial("G4_AIR");
-  auto cwLayers = cw.getLayers();
+  auto cwLayers = cw->getLayers();
   for(auto layer : cwLayers){
       nistManager->FindOrBuildMaterial(layer.material);
   }
@@ -128,14 +130,14 @@ nistManager->FindOrBuildMaterial("G4_AIR");
 G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
 {
   // Geometry parameters
-  auto & cwLayers = cw.getLayers();
+  auto & cwLayers = cw->getLayers();
   G4int nofLayers = cwLayers.size();
   
   G4double caloLength = 0;
   for(auto layer : cwLayers){
       caloLength += layer.thickness * cm;
   }
-  G4double calorSizeXY = cw.getXYWidth() * cm;
+  G4double calorSizeXY = cw->getXYWidth() * cm;
 
   auto worldSizeXY = 1.2 * calorSizeXY;
   auto worldSizeZ  = 1.2 * caloLength;
@@ -210,6 +212,21 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
                      defaultMaterial, //G4Material::GetMaterial(layer.material),  // its material
                      layer.name);         // its name
 
+      //set layerLV to be visible and green if active and transparent otherwise
+      if(layer.isActive){
+          auto layerVisAtt = new G4VisAttributes(G4Colour(0.0,1.0,0.0));
+          layerVisAtt->SetVisibility(true);
+          layerVisAtt->SetForceSolid(true);
+          layerLV->SetVisAttributes(layerVisAtt);
+      } //else set the absorbers to be transparent
+      else{
+          auto layerVisAtt = new G4VisAttributes(G4Colour(1.0,1.0,1.0,0.1));
+          layerVisAtt->SetVisibility(true);
+          layerVisAtt->SetForceSolid(true);
+          layerLV->SetVisAttributes(layerVisAtt);
+      }
+
+
       auto sensorLV = new G4LogicalVolume(
           new G4Box("sensor", layer.sens_xwidth/2*cm, layer.sens_ywidth/2*cm, layer.thickness/2*cm),
           G4Material::GetMaterial(layer.material),
@@ -231,7 +248,20 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
         0,                        // copy number
         fCheckOverlaps);          // checking overlaps
       
-      layer.assignPhysicalVolume(pv);
+      layer.assignPhysicalVolume(ppv); //maybe this needs to be ppv - check
+
+      // assign sensors to layer by copyNumber; access the copyNumber of the volumes
+      if(layer.isActive){
+         for(int i = 0; i < layer.nx*layer.ny; i++){
+             Sensor sensor;
+             sensor.position = G4ThreeVector(layer.sens_xwidth * (i % layer.nx) - layer.sens_xwidth * ((float)layer.nx-1) / 2.*cm, 
+             layer.sens_ywidth * (i / layer.nx) - layer.sens_ywidth * ((float)layer.ny-1) / 2.*cm, position);
+             sensor.size = G4ThreeVector(layer.sens_xwidth, layer.sens_ywidth, layer.thickness);
+             sensor.energy = 0;
+             layer.sensors.push_back(sensor);//this should now be aligned with copy number
+         }
+      }
+
       position += layer.thickness / 2 *cm; //assign the physical volume to the layer
   }
   //
@@ -239,10 +269,11 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
   //
   worldLV->SetVisAttributes (G4VisAttributes::GetInvisible());
 
-  auto simpleBoxVisAtt= new G4VisAttributes(G4Colour(1.0,1.0,1.0));
+  auto simpleBoxVisAtt= new G4VisAttributes(G4Colour(1.0,1.0,1.0,0.));
   simpleBoxVisAtt->SetVisibility(true);
   calorLV->SetVisAttributes(simpleBoxVisAtt);
 
+
   //
   // Always return the physical World
   //
diff --git a/src/EventAction.cc b/src/EventAction.cc
index 30b804e..5bbf3a7 100644
--- a/src/EventAction.cc
+++ b/src/EventAction.cc
@@ -27,6 +27,8 @@
 /// \file B4/B4a/src/EventAction.cc
 /// \brief Implementation of the B4a::EventAction class
 
+#include "ConstructionWrapper.hh"
+
 #include "EventAction.hh"
 #include "RunAction.hh"
 
@@ -34,6 +36,7 @@
 #include "G4RunManager.hh"
 #include "G4Event.hh"
 #include "G4UnitsTable.hh"
+#include "PrimaryGeneratorAction.hh"
 
 #include "Randomize.hh"
 #include <iomanip>
@@ -46,10 +49,28 @@ namespace B4a
 void EventAction::BeginOfEventAction(const G4Event* /*event*/)
 {
   // initialisation per event
-  fEnergyAbs = 0.;
-  fEnergyGap = 0.;
-  fTrackLAbs = 0.;
-  fTrackLGap = 0.;
+  if(cw == nullptr){
+    G4cout << "ConstructionWrapper not found" << G4endl;
+    throw std::runtime_error("ConstructionWrapper not found");
+  }
+  if(gen == nullptr){
+    G4cout << "PrimaryGeneratorAction not found" << G4endl;
+    throw std::runtime_error("PrimaryGeneratorAction not found");
+  }
+  cw->resetSensorEnergies();
+
+  if(run->hitEnergy.size() != cw->getNSensors()){
+      run->hitEnergy.resize(cw->getNSensors(),0);
+      run->hitX.resize(cw->getNSensors(),0);
+      run->hitY.resize(cw->getNSensors(),0);  
+      run->hitZ.resize(cw->getNSensors(),0);  
+      run->hitDX.resize(cw->getNSensors(),0);
+      run->hitDY.resize(cw->getNSensors(),0); 
+      run->hitDZ.resize(cw->getNSensors(),0);
+      run->hitLayer.resize(cw->getNSensors(),0);
+      run->hitCopyNumber.resize(cw->getNSensors(),0);
+  }
+
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -62,38 +83,47 @@ void EventAction::EndOfEventAction(const G4Event* event)
   // get analysis manager
   auto analysisManager = G4AnalysisManager::Instance();
 
-  // fill histograms
-  analysisManager->FillH1(0, fEnergyAbs);
-  analysisManager->FillH1(1, fEnergyGap);
-  analysisManager->FillH1(2, fTrackLAbs);
-  analysisManager->FillH1(3, fTrackLGap);
+
+  //cw->printSensorEnergies();//DEBUG
+
+  std::vector<double> sensorEnergies;
+  double tot_energy=0;
+  int sensorNumber=0;
+  int layerNumber=0;
+  for(const auto& layer: cw->getLayers()){
+    for(const auto& sensor: layer.sensors){
+      sensorEnergies.push_back(sensor.getEnergy());
+      tot_energy+= sensor.getEnergy();
+      
+      //fill the vectors
+      run->hitEnergy[sensorNumber] = sensor.getEnergy();
+      run->hitX[sensorNumber] = sensor.getPos().x();
+      run->hitY[sensorNumber] = sensor.getPos().y();
+      run->hitZ[sensorNumber] = sensor.getPos().z();
+      run->hitDX[sensorNumber] = sensor.getSize().x();
+      run->hitDY[sensorNumber] = sensor.getSize().y();
+      run->hitDZ[sensorNumber] = sensor.getSize().z();
+      run->hitLayer[sensorNumber] = layerNumber;
+      run->hitCopyNumber[sensorNumber] = sensorNumber;
+       sensorNumber++;
+    }
+    layerNumber++;
+  }
+  //assign the vectors
+
 
   // fill ntuple
-  analysisManager->FillNtupleDColumn(0, fEnergyAbs);
-  analysisManager->FillNtupleDColumn(1, fEnergyGap);
-  analysisManager->FillNtupleDColumn(2, fTrackLAbs);
-  analysisManager->FillNtupleDColumn(3, fTrackLGap);
+  analysisManager->FillNtupleDColumn(0, gen->getPartEnergy()); //needs to be true energy DEBUG
+  analysisManager->FillNtupleDColumn(1, tot_energy);
+  analysisManager->FillNtupleIColumn(2, (int)cw->getLayers().size());
+  analysisManager->FillNtupleIColumn(3, cw->getNSensors());
+
+  //vectors are stored automatically
   analysisManager->AddNtupleRow();
 
-  // Print per event (modulo n)
-  //
-  auto eventID = event->GetEventID();
-  auto printModulo = G4RunManager::GetRunManager()->GetPrintProgress();
-  if ( ( printModulo > 0 ) && ( eventID % printModulo == 0 ) ) {
-    G4cout << "---> End of event: " << eventID << G4endl;
 
-    G4cout
-       << "   Absorber: total energy: " << std::setw(7)
-                                        << G4BestUnit(fEnergyAbs,"Energy")
-       << "       total track length: " << std::setw(7)
-                                        << G4BestUnit(fTrackLAbs,"Length")
-       << G4endl
-       << "        Gap: total energy: " << std::setw(7)
-                                        << G4BestUnit(fEnergyGap,"Energy")
-       << "       total track length: " << std::setw(7)
-                                        << G4BestUnit(fTrackLGap,"Length")
-       << G4endl;
-  }
+
+
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/src/PrimaryGeneratorAction.cc b/src/PrimaryGeneratorAction.cc
index 793ee4e..88b0477 100644
--- a/src/PrimaryGeneratorAction.cc
+++ b/src/PrimaryGeneratorAction.cc
@@ -39,6 +39,8 @@
 #include "G4ParticleDefinition.hh"
 #include "G4SystemOfUnits.hh"
 #include "Randomize.hh"
+#include "G4INCLRandom.hh"
+#include <G4INCLGeant4Random.hh>
 
 namespace B4
 {
@@ -50,13 +52,15 @@ PrimaryGeneratorAction::PrimaryGeneratorAction()
   G4int nofParticles = 1;
   fParticleGun = new G4ParticleGun(nofParticles);
 
+  partSpecies = "e-";
+  minPartEnergy = 1000;    
+  maxPartEnergy = 1000;
+
+  G4INCL::Random::setGenerator( new G4INCL::Geant4RandomGenerator());
+
+
   // default particle kinematic
   //
-  auto particleDefinition
-    = G4ParticleTable::GetParticleTable()->FindParticle("e-");
-  fParticleGun->SetParticleDefinition(particleDefinition);
-  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
-  fParticleGun->SetParticleEnergy(50.*MeV);
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -70,8 +74,19 @@ PrimaryGeneratorAction::~PrimaryGeneratorAction()
 
 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
 {
+
+
   // This function is called at the begining of event
 
+  G4double rand =  G4INCL::Random::shoot();
+  partEnergy = rand*(maxPartEnergy-minPartEnergy)+minPartEnergy;
+
+  auto particleDefinition
+    = G4ParticleTable::GetParticleTable()->FindParticle(partSpecies);
+  fParticleGun->SetParticleDefinition(particleDefinition);
+  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
+  fParticleGun->SetParticleEnergy(partEnergy*MeV);
+
   // In order to avoid dependence of PrimaryGeneratorAction
   // on DetectorConstruction class we get world volume
   // from G4LogicalVolumeStore
@@ -102,6 +117,7 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
     ->SetParticlePosition(G4ThreeVector(0., 0., -worldZHalfLength));
 
   fParticleGun->GeneratePrimaryVertex(anEvent);
+
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/src/RunAction.cc b/src/RunAction.cc
index ba1bcb3..d39bd6f 100644
--- a/src/RunAction.cc
+++ b/src/RunAction.cc
@@ -60,19 +60,26 @@ RunAction::RunAction()
   // Book histograms, ntuple
   //
 
-  // Creating histograms
-  analysisManager->CreateH1("Eabs","Edep in absorber", 100, 0., 800*MeV);
-  analysisManager->CreateH1("Egap","Edep in gap", 100, 0., 100*MeV);
-  analysisManager->CreateH1("Labs","trackL in absorber", 100, 0., 1*m);
-  analysisManager->CreateH1("Lgap","trackL in gap", 100, 0., 50*cm);
-
   // Creating ntuple
   //
-  analysisManager->CreateNtuple("B4", "Edep and TrackL");
-  analysisManager->CreateNtupleDColumn("Eabs");
-  analysisManager->CreateNtupleDColumn("Egap");
-  analysisManager->CreateNtupleDColumn("Labs");
-  analysisManager->CreateNtupleDColumn("Lgap");
+  analysisManager->CreateNtuple("Hits", "Hits");
+
+  analysisManager->CreateNtupleDColumn("true_energy");
+  analysisManager->CreateNtupleDColumn("total_dep_energy");
+  analysisManager->CreateNtupleIColumn("N_layers");
+  analysisManager->CreateNtupleIColumn("N_sensors");
+
+  analysisManager->CreateNtupleDColumn("sensor_energy",hitEnergy);
+  analysisManager->CreateNtupleDColumn("sensor_x",hitX);
+  analysisManager->CreateNtupleDColumn("sensor_y",hitY);
+  analysisManager->CreateNtupleDColumn("sensor_z",hitZ);
+  analysisManager->CreateNtupleDColumn("sensor_dx",hitDX);
+  analysisManager->CreateNtupleDColumn("sensor_dy",hitDY);
+  analysisManager->CreateNtupleDColumn("sensor_dz",hitDZ);
+  analysisManager->CreateNtupleIColumn("sensor_layer",hitLayer);
+  analysisManager->CreateNtupleIColumn("sensor_copy_number",hitCopyNumber);
+
+
   analysisManager->FinishNtuple();
 }
 
@@ -104,35 +111,7 @@ void RunAction::EndOfRunAction(const G4Run* /*run*/)
   // print histogram statistics
   //
   auto analysisManager = G4AnalysisManager::Instance();
-  if ( analysisManager->GetH1(1) ) {
-    G4cout << G4endl << " ----> print histograms statistic ";
-    if(isMaster) {
-      G4cout << "for the entire run " << G4endl << G4endl;
-    }
-    else {
-      G4cout << "for the local thread " << G4endl << G4endl;
-    }
-
-    G4cout << " EAbs : mean = "
-       << G4BestUnit(analysisManager->GetH1(0)->mean(), "Energy")
-       << " rms = "
-       << G4BestUnit(analysisManager->GetH1(0)->rms(),  "Energy") << G4endl;
-
-    G4cout << " EGap : mean = "
-       << G4BestUnit(analysisManager->GetH1(1)->mean(), "Energy")
-       << " rms = "
-       << G4BestUnit(analysisManager->GetH1(1)->rms(),  "Energy") << G4endl;
-
-    G4cout << " LAbs : mean = "
-      << G4BestUnit(analysisManager->GetH1(2)->mean(), "Length")
-      << " rms = "
-      << G4BestUnit(analysisManager->GetH1(2)->rms(),  "Length") << G4endl;
-
-    G4cout << " LGap : mean = "
-      << G4BestUnit(analysisManager->GetH1(3)->mean(), "Length")
-      << " rms = "
-      << G4BestUnit(analysisManager->GetH1(3)->rms(),  "Length") << G4endl;
-  }
+  
 
   // save histograms & ntuple
   //
diff --git a/src/SteppingAction.cc b/src/SteppingAction.cc
index 0678b73..eeabfff 100644
--- a/src/SteppingAction.cc
+++ b/src/SteppingAction.cc
@@ -59,19 +59,27 @@ void SteppingAction::UserSteppingAction(const G4Step* step)
   // energy deposit
   auto edep = step->GetTotalEnergyDeposit();
 
-  // step length
-  G4double stepLength = 0.;
-  if ( step->GetTrack()->GetDefinition()->GetPDGCharge() != 0. ) {
-    stepLength = step->GetStepLength();
+  //return;
+  
+  //find the layer through the constructionwrapper that corresponds to the volume
+  auto cw = fDetConstruction->getConstructionWrapper();
+  if(cw == nullptr){
+    G4cout << "ConstructionWrapper not found" << G4endl;
+    throw std::runtime_error("ConstructionWrapper not found");
   }
 
-  if ( volume == fDetConstruction->GetAbsorberPV() ) {
-    fEventAction->AddAbs(edep,stepLength);
+  
+  auto sensor = fDetConstruction->getConstructionWrapper()->getSensorByVolume(volume);
+  if(sensor == nullptr){
+    //can simply be not an active volume, so no need to throw an error
+    //G4cout << "Sensor not found in" << volume->GetName() << G4endl; //DEBUG
+    return;
   }
+  //G4cout << "Sensor found in " << volume->GetName() << " with copy number " << volume->GetCopyNo()<< G4endl; //DEBUG
+  sensor->energy += edep;
+  //find sensor by volume copy number
+
 
-  if ( volume == fDetConstruction->GetGapPV() ) {
-    fEventAction->AddGap(edep,stepLength);
-  }
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......