diff --git a/Digitisers/SimpleDigi/CMakeLists.txt b/Digitisers/SimpleDigi/CMakeLists.txt
index 4e2687bc123eeb14f5b013af3cf0461ec48605db..28c8ea89d61d17c7d123b9d1366e94f49ac7aa29 100644
--- a/Digitisers/SimpleDigi/CMakeLists.txt
+++ b/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
diff --git a/Digitisers/SimpleDigi/src/CylinderDigiAlg.cpp b/Digitisers/SimpleDigi/src/CylinderDigiAlg.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d3c67e60c89c427ce12a7c4bc27561b1b971f120
--- /dev/null
+++ b/Digitisers/SimpleDigi/src/CylinderDigiAlg.cpp
@@ -0,0 +1,119 @@
+#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();
+}
diff --git a/Digitisers/SimpleDigi/src/CylinderDigiAlg.h b/Digitisers/SimpleDigi/src/CylinderDigiAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..af3ffe2100f8151e0076c7ff4643f369f9f836ab
--- /dev/null
+++ b/Digitisers/SimpleDigi/src/CylinderDigiAlg.h
@@ -0,0 +1,39 @@
+#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
diff --git a/Simulation/DetSimSD/CMakeLists.txt b/Simulation/DetSimSD/CMakeLists.txt
index fa8ba9ef3e13d27da31d8c25763278ba795ef47c..210e3639bd4ea751e65e55ba8fc87955dc56fb96 100644
--- a/Simulation/DetSimSD/CMakeLists.txt
+++ b/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}
diff --git a/Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp b/Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp
index 43dfa17828453c1545318eea45c0aec7e271976e..c911d47071bfd1410708fbebc1bd5d3d2501954d 100644
--- a/Simulation/DetSimSD/src/DriftChamberSensDetTool.cpp
+++ b/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;
}
diff --git a/Simulation/DetSimSD/src/GenericTrackerSensDetTool.cpp b/Simulation/DetSimSD/src/GenericTrackerSensDetTool.cpp
index 480f509aa57c563e90c032aac64abcaca4957a31..edb275a8d76e4d0710368611e916f4c026837e12 100644
--- a/Simulation/DetSimSD/src/GenericTrackerSensDetTool.cpp
+++ b/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;
}
diff --git a/Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.cpp b/Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0d53d3b99c915a4428be09f2dfceec72040f7195
--- /dev/null
+++ b/Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.cpp
@@ -0,0 +1,56 @@
+#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();
+}
diff --git a/Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.h b/Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.h
new file mode 100644
index 0000000000000000000000000000000000000000..cf0cb844eb64a076d8c39477c94a6fd62f677402
--- /dev/null
+++ b/Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.h
@@ -0,0 +1,95 @@
+// *********************************************************
+//
+// $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