add segementation to the calorimeter

draw_simhits.C

@@ -5,14 +5,15 @@
 
 /** Helper function for radius
  */
-float r(float x,float y) { return sqrt(x*x+y*y) ; }
+float radius(float x,float y) { return sqrt(x*x+y*y) ; }
 
 
 /** helper function for layer number
  */
 int layer(int cellid){
-  // so far only layer used in cellid
-  return ( cellid ) ;
+  //  encoding used in SimpleCalo:
+  //  1000000 * (_xID+1) +  1000 * (_yID+1)  + (_zID+1)
+  return ( cellid % 1000 - 1 ) ;
 }
 
 
@@ -25,10 +26,10 @@ void draw_simhits( const char* FILEN ) {
   auto* events = (TTree*) f->Get("events") ;
   
   // define some aliases for simpler plotting
-  events->SetAlias("x","SimCalorimeterHits.position.x");
-  events->SetAlias("y","SimCalorimeterHits.position.y");
-  events->SetAlias("z","SimCalorimeterHits.position.z");
-  events->SetAlias("l","SimCalorimeterHits.cellID");
+  events->SetAlias("xx","SimCalorimeterHits.position.x");
+  events->SetAlias("yy","SimCalorimeterHits.position.y");
+  events->SetAlias("zz","SimCalorimeterHits.position.z");
+  events->SetAlias("id","SimCalorimeterHits.cellID");
   events->SetAlias("e","SimCalorimeterHits.energy");
 
 
@@ -43,15 +44,15 @@ void draw_simhits( const char* FILEN ) {
   
   c1->cd(2) ;
   // hits scatter plot in x-y
-  events->Draw("x:y","abs(x)<50&abs(y)<50");
+  events->Draw("xx:yy","abs(xx)<50&abs(yy)<50");
 
   c1->cd(3) ;
   // longitudenal profile
-  events->Draw("l","e") ;
+  events->Draw("layer(id)","e") ;
 
   c1->cd(4) ;
   // radial profile
-  events->Draw("sqrt(x*x+y*y)","(sqrt(x*x+y*y)<50)*e") ;
+  events->Draw("radius(xx,yy)","(radius(xx,yy)<50)*e") ;
   
   
   std::string pdfFile( FILEN ) ;

include/SimpleCaloConfig.hh

@@ -39,6 +39,9 @@ public:
   /// number of calorimeter layers
   G4int nLayers =  100 ;
 
+  /// cell size
+  G4double cellSize = 5.0 * mm ;
+  
   /// width (in x-y) of the calorimter
   G4double width    =  100.*cm ; 
 

src/SimpleCaloEDM4hepEventAction.cc

@@ -31,6 +31,8 @@
 
 #include "SimpleCaloEDM4hepEventAction.hh"
 #include "SimpleCaloHit.hh"
+#include "SimpleCaloConfig.hh"
+
 #include "G4EventManager.hh"
 #include "G4SDManager.hh"
 #include "G4UImanager.hh"
@@ -40,6 +42,8 @@
 //std
 #include <iostream>
 #include <algorithm>
+#include <cmath>
+
 //Gflash
 using namespace std;
 
@@ -52,6 +56,54 @@ using namespace std;
 #include "edm4hep/SimCalorimeterHitCollection.h"
 
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+// helper class for calorimeter segmentation (cartesian x-y )
+
+class SimpleCaloSegementation {
+
+  G4ThreeVector _hitPosition = {} ;
+  edm4hep::Vector3f _cellPosition = {} ;
+  int _nCellsX = 0 ;
+  int _nCellsY = 0 ;
+
+  int _xID = 0 ;
+  int _yID = 0 ;
+  int _zID = 0 ;
+
+  int _nLayers = SimpleCaloConfig::get().nLayers ;
+  double _absTh = SimpleCaloConfig::get().absThickness ;
+  double _senTh = SimpleCaloConfig::get().sensThickness ;
+  
+public:
+  SimpleCaloSegementation(){
+
+    _nCellsX = round( SimpleCaloConfig::get().width / SimpleCaloConfig::get().cellSize );
+    _nCellsY = round( SimpleCaloConfig::get().width / SimpleCaloConfig::get().cellSize );
+
+  }
+  
+  void setPositionAndLayer( const G4ThreeVector& pos, int layer ) {
+    _hitPosition = pos ;
+    _xID = int( _nCellsX / 2 ) + _hitPosition.x() / SimpleCaloConfig::get().cellSize ;
+    _yID = int( _nCellsY / 2 ) + _hitPosition.y() / SimpleCaloConfig::get().cellSize ;
+    _zID = layer ;
+
+    _cellPosition = {
+      float( ( _xID -  int( _nCellsX / 2 ) + 0.5 ) * SimpleCaloConfig::get().cellSize ),
+      float( ( _yID -  int( _nCellsY / 2 ) + 0.5 ) * SimpleCaloConfig::get().cellSize ),
+      float( SimpleCaloConfig::get().zStart + _zID * ( _absTh + _senTh ) + _absTh + _senTh/2.) 
+    } ;
+  }
+
+  unsigned long getCellID(){ return  1000000 * (_xID+1) +  1000 * (_yID+1)  + (_zID+1) ; }
+
+  edm4hep::Vector3f getCellPosition() { return _cellPosition ; }
+
+};
+
+
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 SimpleCaloEDM4hepEventAction::SimpleCaloEDM4hepEventAction()
@@ -118,6 +170,9 @@ void SimpleCaloEDM4hepEventAction::EndOfEventAction(const G4Event *evt)
   if (fCalorimeterCollectionId<0) return;
   G4HCofThisEvent * HCE = evt->GetHCofThisEvent();
   SimpleCaloHitsCollection* THC = 0;
+
+  SimpleCaloSegementation seg ;
+
   G4double totE = 0;
   // Read out of the calo
   THC=(SimpleCaloHitsCollection *)(HCE->GetHC(fCalorimeterCollectionId));
@@ -152,15 +207,18 @@ void SimpleCaloEDM4hepEventAction::EndOfEventAction(const G4Event *evt)
 	
 	auto hit = fSHC->create() ;
 	
-	hit.setCellID( (*THC)[i]->GetLayerNum() ) ;  // fixme: add x/y coordinates ....
+	int layerNum = (*THC)[i]->GetLayerNum()  ;
+	G4ThreeVector hitpos=(*THC)[i]->GetPos();            
+
+	seg.setPositionAndLayer( hitpos , layerNum ) ;
+
+	hit.setCellID( seg.getCellID() ) ;
 	
 	totE += estep ;
 	hit.setEnergy( estep ) ;
 	
-	G4ThreeVector hitpos=(*THC)[i]->GetPos();            
-	G4ThreeVector l (hitpos.x(), hitpos.y(), hitpos.z());
-	
-	hit.setPosition( { (float) hitpos.x(), (float) hitpos.y(), (float) hitpos.z() } ) ;
+//	hit.setPosition( { (float) hitpos.x(), (float) hitpos.y(), (float) hitpos.z() } ) ;
+	hit.setPosition( seg.getCellPosition() ) ;
 	
       }
     }