Merge pull request #221 from fucd/dc

SIM: digi of cylinder sensor & combined sensistive detector 

Digitisers/SimpleDigi/CMakeLists.txt

@@ -3,6 +3,7 @@ gaudi_add_module(SimpleDigi
                  SOURCES src/PlanarDigiAlg.cpp
                          src/TPCDigiAlg.cpp
                          src/voxel.cpp
+                         src/CylinderDigiAlg.cpp
                  LINK GearSvc
                       EventSeeder
                       TrackSystemSvcLib

Digitisers/SimpleDigi/src/CylinderDigiAlg.cpp

+#include "CylinderDigiAlg.h"
+
+#include "edm4hep/Vector3f.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include "DD4hep/DD4hepUnits.h"
+#include "DDRec/Vector3D.h"
+
+#include "GaudiKernel/INTupleSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/IRndmGen.h"
+#include "GaudiKernel/IRndmGenSvc.h"
+#include "GaudiKernel/RndmGenerators.h"
+
+DECLARE_COMPONENT( CylinderDigiAlg )
+
+CylinderDigiAlg::CylinderDigiAlg(const std::string& name, ISvcLocator* svcLoc)
+: GaudiAlgorithm(name, svcLoc){
+  // Input collections
+  declareProperty("InputSimTrackerHitCollection", m_inputColHdls, "Handle of the Input SimTrackerHit collection");
+
+  // Output collections
+  declareProperty("OutputTrackerHitCollection", m_outputColHdls, "Handle of the output TrackerHit collection");
+  declareProperty("TrackerHitAssociationCollection", m_assColHdls, "Handle of the Association collection between SimTrackerHit and TrackerHit");
+}
+
+StatusCode CylinderDigiAlg::initialize(){
+  m_geosvc = service<IGeomSvc>("GeomSvc");
+  if(!m_geosvc){
+    error() << "Failed to get the GeomSvc" << endmsg;
+    return StatusCode::FAILURE;
+  }
+  auto detector = m_geosvc->lcdd();
+  if(!detector){
+    error() << "Failed to get the Detector from GeomSvc" << endmsg;
+    return StatusCode::FAILURE;
+  }
+  std::string name = m_inputColHdls.objKey(); 
+  debug() << "Readout name: " << name << endmsg;
+  m_decoder = m_geosvc->getDecoder(name);
+  if(!m_decoder){
+    error() << "Failed to get the decoder. " << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  info() << "CylinderDigiAlg::initialized" << endmsg;
+  return GaudiAlgorithm::initialize();
+}
+
+
+StatusCode CylinderDigiAlg::execute(){
+  auto trkhitVec = m_outputColHdls.createAndPut();
+  auto assVec = m_assColHdls.createAndPut();
+
+  const edm4hep::SimTrackerHitCollection* STHCol = nullptr;
+  try {
+    STHCol = m_inputColHdls.get();
+  }
+  catch(GaudiException &e){
+    debug() << "Collection " << m_inputColHdls.fullKey() << " is unavailable in event " << m_nEvt << endmsg;
+    return StatusCode::SUCCESS;
+  }
+  if(STHCol->size()==0) return StatusCode::SUCCESS;
+  debug() << m_inputColHdls.fullKey() << " has SimTrackerHit "<< STHCol->size() << endmsg;
+  
+  for(auto simhit : *STHCol){
+    auto particle = simhit.getMCParticle();
+    if(!particle.isAvailable()) continue;
+    
+    auto& mom0 = particle.getMomentum();
+    double pt = sqrt(mom0.x*mom0.x+mom0.y*mom0.y);
+    if(particle.isCreatedInSimulation()&&pt<0.01&&particle.isStopped()) continue;
+
+    auto cellId = simhit.getCellID();
+    int system  = m_decoder->get(cellId, "system");
+    int chamber = m_decoder->get(cellId, "chamber");
+    int layer   = m_decoder->get(cellId, "layer"  );
+    auto& pos   = simhit.getPosition();
+    auto& mom   = simhit.getMomentum();
+    
+    double phi = atan2(pos[1], pos[0]);
+    double r   = sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
+    double dphi = m_resRPhi/r;
+    phi += randSvc()->generator(Rndm::Gauss(0, dphi))->shoot();
+    double smearedX = r*cos(phi);
+    double smearedY = r*sin(phi);
+    double smearedZ = pos[2] + randSvc()->generator(Rndm::Gauss(0, m_resZ))->shoot();
+
+    auto trkHit = trkhitVec->create();
+    trkHit.setCellID(cellId);
+    trkHit.setTime(simhit.getTime());
+    trkHit.setEDep(simhit.getEDep());
+    trkHit.setPosition (edm4hep::Vector3d(smearedX, smearedY, smearedZ));
+    trkHit.setCovMatrix(std::array<float, 6>{m_resRPhi*m_resRPhi/2, 0, m_resRPhi*m_resRPhi/2, 0, 0, m_resZ*m_resZ});
+    //trkHit.setType(CEPC::CYLINDER);
+    trkHit.addToRawHits(simhit.getObjectID());
+    debug() << "Hit " << simhit.id() << ": " << pos << " -> " << trkHit.getPosition() << "s:" << system << " c:" << chamber << " l:" << layer
+	    << " pt = " << pt << " " << mom.x << " " << mom.y << " " << mom.z << endmsg;
+
+    auto ass = assVec->create();
+
+    float weight = 1.0;
+
+    debug() <<" Set relation between " << " sim hit " << simhit.id() << " to tracker hit " << trkHit.id() << " with a weight of " << weight << endmsg;
+    ass.setSim(simhit);
+    ass.setRec(trkHit);
+    ass.setWeight(weight);
+  }
+
+  m_nEvt++;
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode CylinderDigiAlg::finalize(){
+  info() << "Processed " << m_nEvt << " events " << endmsg;
+  return GaudiAlgorithm::finalize();
+}

Digitisers/SimpleDigi/src/CylinderDigiAlg.h

+#ifndef CylinderDigiAlg_H
+#define CylinderDigiAlg_H
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiKernel/NTuple.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+
+#include <DDRec/DetectorData.h>
+#include "DetInterface/IGeomSvc.h"
+
+class CylinderDigiAlg : public GaudiAlgorithm{
+ public:
+  
+  CylinderDigiAlg(const std::string& name, ISvcLocator* svcLoc);
+ 
+  virtual StatusCode initialize() ;
+  virtual StatusCode execute() ; 
+  virtual StatusCode finalize() ;
+ 
+protected:
+  
+  SmartIF<IGeomSvc> m_geosvc;
+  dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+
+  // Input collections
+  DataHandle<edm4hep::SimTrackerHitCollection>            m_inputColHdls{"DriftChamberHitsCollection", Gaudi::DataHandle::Reader, this};
+  // Output collections
+  DataHandle<edm4hep::TrackerHitCollection>               m_outputColHdls{"DCTrackerHits", Gaudi::DataHandle::Writer, this};
+  DataHandle<edm4hep::MCRecoTrackerAssociationCollection> m_assColHdls{"DCTrackerHitAssociationCollection", Gaudi::DataHandle::Writer, this};
+
+  Gaudi::Property<float> m_resRPhi{this, "ResRPhi", 0.1};
+  Gaudi::Property<float> m_resZ   {this, "ResZ", 2.828};
+
+  int m_nEvt=0;
+};
+#endif

Simulation/DetSimSD/CMakeLists.txt

@@ -15,6 +15,7 @@ gaudi_add_module(DetSimSD
 
                          src/GenericTrackerSensDetTool.cpp
                          src/GenericTrackerSensitiveDetector.cpp
+                         src/TrackerCombineSensitiveDetector.cpp
                  LINK DetSimInterface
                       DetInterface
                       ${DD4hep_COMPONENT_LIBRARIES} 

Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp

@@ -5,7 +5,7 @@
 #include "DD4hep/Detector.h"
 
 #include "DriftChamberSensitiveDetector.h"
-
+#include "TrackerCombineSensitiveDetector.h"
 DECLARE_COMPONENT(DriftChamberSensDetTool)
 
 StatusCode DriftChamberSensDetTool::initialize() {
@@ -39,12 +39,19 @@ DriftChamberSensDetTool::createSD(const std::string& name) {
     G4VSensitiveDetector* sd = nullptr;
 
     if (name == "DriftChamber") {
-        DriftChamberSensitiveDetector* dcsd = new DriftChamberSensitiveDetector(name, *dd4hep_geo);
-        dcsd->setDedxSimTool(m_dedx_simtool);
-
-        sd = dcsd;
+      auto sens = dd4hep_geo->sensitiveDetector(name);
+      if(!sens.combineHits()){
+	DriftChamberSensitiveDetector* dcsd = new DriftChamberSensitiveDetector(name, *dd4hep_geo);
+	dcsd->setDedxSimTool(m_dedx_simtool);
+	sd = dcsd;
+      }
+      else{
+	sd = new TrackerCombineSensitiveDetector(name, *dd4hep_geo);
+      }
+    }
+    else{
+      warning() << "The SD " << name << " want to use SD for DriftChamber" << endmsg;
     }
-
 
     return sd;
 }

Simulation/DetSimSD/src/GenericTrackerSensDetTool.cpp

@@ -5,6 +5,7 @@
 #include "DD4hep/Detector.h"
 
 #include "GenericTrackerSensitiveDetector.h"
+#include "TrackerCombineSensitiveDetector.h"
 
 #include "CLHEP/Units/SystemOfUnits.h"
 
@@ -33,8 +34,14 @@ G4VSensitiveDetector* GenericTrackerSensDetTool::createSD(const std::string& nam
   debug() << "createSD for " << name << endmsg;
   
   dd4hep::Detector* dd4hep_geo = m_geosvc->lcdd();
-
-  G4VSensitiveDetector* sd = new GenericTrackerSensitiveDetector(name, *dd4hep_geo);
+  auto sens = dd4hep_geo->sensitiveDetector(name);
+  G4VSensitiveDetector* sd = nullptr;
+  if(sens.combineHits()){
+    sd = new TrackerCombineSensitiveDetector(name, *dd4hep_geo);
+  }
+  else{
+    sd = new GenericTrackerSensitiveDetector(name, *dd4hep_geo);
+  }
 
   return sd;
 }

Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.cpp

+#include "TrackerCombineSensitiveDetector.h"
+
+#include "G4Step.hh"
+#include "G4VProcess.hh"
+#include "G4SDManager.hh"
+#include "DD4hep/DD4hepUnits.h"
+
+TrackerCombineSensitiveDetector::TrackerCombineSensitiveDetector(const std::string& name,
+								 dd4hep::Detector& description)
+  : DDG4SensitiveDetector(name, description),
+    m_hc(nullptr){
+  G4String CollName = m_sensitive.hitsCollection();
+  collectionName.insert(CollName);
+}
+
+void TrackerCombineSensitiveDetector::Initialize(G4HCofThisEvent* HCE){
+  userData.e_cut = m_sensitive.energyCutoff();
+
+  m_hc = new HitCollection(GetName(), collectionName[0]);
+  int HCID = G4SDManager::GetSDMpointer()->GetCollectionID(m_hc);
+  HCE->AddHitsCollection( HCID, m_hc ); 
+}
+
+G4bool TrackerCombineSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*){
+  dd4hep::sim::Geant4StepHandler h(step);
+  bool return_code = false;
+  if ( userData.current == -1 ) userData.start(getCellID(step), step, h.pre);
+  else if ( !userData.track || userData.current != h.track->GetTrackID() ) {
+    return_code = userData.extractHit(m_hc) != 0;
+    userData.start(getCellID(step), step, h.pre);
+  }
+
+  // ....update .....
+  userData.update(step);
+
+  void *prePV = h.volume(h.pre), *postPV = h.volume(h.post);
+  if ( prePV != postPV ) {
+    return_code = userData.extractHit(m_hc) != 0;
+    void* postSD = h.sd(h.post);
+    if ( 0 != postSD )   {
+      void* preSD = h.sd(h.pre);
+      if ( preSD == postSD ) {
+	// fucd: getCellID(step) for preVolume not postVolume, so should start at next step 
+	//userData.start(getCellID(step), step, h.post);
+      }
+    }
+  }
+  else if ( userData.track->GetTrackStatus() == fStopAndKill ) {
+    return_code = userData.extractHit(m_hc) != 0;
+  }
+  return return_code;
+}
+
+void TrackerCombineSensitiveDetector::EndOfEvent(G4HCofThisEvent* HCE){
+  userData.clear();
+}

Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.h

+// *********************************************************
+//
+// $Id: TrackerCombineSensitiveDetector.hh,v 1.0 2022/03/27
+
+#ifndef TrackerCombineSensitiveDetector_h
+#define TrackerCombineSensitiveDetector_h
+
+#include "DetSimSD/DDG4SensitiveDetector.h"
+#include "DDG4/Defs.h"
+
+namespace dd4hep {
+  namespace sim {
+    struct TrackerCombine {
+      Geant4TrackerHit  pre;
+      Geant4TrackerHit  post;
+      G4Track*          track;
+      double            e_cut;
+      int               current;
+      long long int     cellID;
+      TrackerCombine() : pre(), post(), track(0), e_cut(0.0), current(-1), cellID(0)  {}
+      void start(long long int cell, G4Step* step, G4StepPoint* point)   {
+	pre.storePoint(step,point);
+	current = pre.truth.trackID;
+	track   = step->GetTrack();
+	cellID  = cell;
+	post    = pre;
+	//std::cout << "start: " << cellID << " " << pre.position << " current = " << current << std::endl;
+      }
+      void update(G4Step* step) {
+	post.storePoint(step,step->GetPostStepPoint());
+	pre.truth.deposit += post.truth.deposit;
+	//std::cout << "update: " << cellID << " " << post.position << " current = " << current << std::endl;
+      }
+      void clear()   {
+	pre.truth.clear();
+	current = -1;
+	track = 0;
+      }
+      Geant4TrackerHit* extractHit(DDG4SensitiveDetector::HitCollection* c) {
+	//std::cout << "create Geant4TrackerHit: " << cellID << " current = " << current << " track = " << track << " de = " << pre.truth.deposit << std::endl;  
+	if ( current == -1 || !track ) {
+	  return 0;
+	}
+	else if ( pre.truth.deposit <= e_cut && !Geant4Hit::isGeantino(track) ) {
+	  clear();
+	  return 0;
+	}
+	double rho1 = pre.position.Rho();
+	double rho2 = post.position.Rho();
+	double rho = 0.5*(rho1+rho2);
+	Position pos = 0.5 * (pre.position + post.position);
+	double z = pos.z();
+	double r = sqrt(rho*rho+z*z);
+	Position path = post.position - pre.position;
+	double angle_O_pre_post = acos(-pre.position.Unit().Dot(path.Unit()));
+	double angle_O_post_pre = acos(post.position.Unit().Dot(path.Unit()));
+	double angle_O_P_pre = asin(pre.position.R()*sin(angle_O_pre_post)/r);
+	if(angle_O_pre_post>dd4hep::halfpi||angle_O_post_pre>dd4hep::halfpi){
+	  bool backward = angle_O_post_pre>angle_O_pre_post;
+	  double angle_O_P_pre = backward ? dd4hep::pi - asin(pre.position.R()*sin(angle_O_pre_post)/r) : asin(pre.position.R()*sin(angle_O_pre_post)/r);
+	  double pre2P = r/sin(angle_O_pre_post)*sin(angle_O_pre_post+angle_O_P_pre);
+	  pos = pre.position + pre2P*path.Unit();
+	}
+	Momentum mom = 0.5 * (pre.momentum + post.momentum);
+	Geant4TrackerHit* hit = new Geant4TrackerHit(pre.truth.trackID,
+						     pre.truth.pdgID,
+						     pre.truth.deposit,
+						     pre.truth.time);
+	hit->cellID   = cellID;
+	hit->position = pos;
+	hit->momentum = mom;
+	hit->length   = path.R();;
+	clear();
+	c->insert(hit);
+	return hit;
+      }
+    };
+  }
+}
+
+class TrackerCombineSensitiveDetector: public DDG4SensitiveDetector {
+ public:
+  TrackerCombineSensitiveDetector(const std::string& name, dd4hep::Detector& description);
+  
+  void Initialize(G4HCofThisEvent* HCE);
+  G4bool ProcessHits(G4Step* step, G4TouchableHistory* history);
+  void EndOfEvent(G4HCofThisEvent* HCE);
+  
+ protected:
+  HitCollection* m_hc = nullptr;
+
+ private:
+  dd4hep::sim::TrackerCombine userData;
+};
+#endif