From e537281d4d1ed1abafda1c518904882bcef5fdcf Mon Sep 17 00:00:00 2001
From: Chengdong Fu <fucd@ihep.ac.cn>
Date: Mon, 30 May 2022 09:39:31 +0800
Subject: [PATCH] new DC module for MarlinTrk&TOT
---
.../CRD_common_v01/DC_Straight_v01_01.xml | 41 +++
.../Tracker/StraightDriftChamber_v01_geo.cpp | 252 ++++++++++++++++++
Digitisers/SimpleDigi/CMakeLists.txt | 1 +
Digitisers/SimpleDigi/src/CylinderDigiAlg.cpp | 120 +++++++++
Digitisers/SimpleDigi/src/CylinderDigiAlg.h | 39 +++
Simulation/DetSimSD/CMakeLists.txt | 1 +
.../DetSimSD/src/DriftChamberSensDetTool.cpp | 18 +-
.../DetSimSD/src/DriftChamberSensDetTool.h | 1 +
.../src/TrackerCombineSensitiveDetector.cpp | 56 ++++
.../src/TrackerCombineSensitiveDetector.h | 95 +++++++
10 files changed, 618 insertions(+), 6 deletions(-)
create mode 100644 Detector/DetCRD/compact/CRD_common_v01/DC_Straight_v01_01.xml
create mode 100644 Detector/DetCRD/src/Tracker/StraightDriftChamber_v01_geo.cpp
create mode 100644 Digitisers/SimpleDigi/src/CylinderDigiAlg.cpp
create mode 100644 Digitisers/SimpleDigi/src/CylinderDigiAlg.h
create mode 100644 Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.cpp
create mode 100644 Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.h
diff --git a/Detector/DetCRD/compact/CRD_common_v01/DC_Straight_v01_01.xml b/Detector/DetCRD/compact/CRD_common_v01/DC_Straight_v01_01.xml
new file mode 100644
index 00000000..2502ab45
--- /dev/null
+++ b/Detector/DetCRD/compact/CRD_common_v01/DC_Straight_v01_01.xml
@@ -0,0 +1,41 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<lccdd>
+ <info name="DriftChamber" title="Drift Chamber only with straight wire" url="http://github.com/cepc/CEPCSW" author="FU Chengdong" status="development" version="v1">
+ <comment>Test with Drift Chamber</comment>
+ </info>
+
+ <define>
+ <constant name="DC_cell_size" value="15*mm"/>
+ <constant name="DC_layer_number" value="int((DC_chamber_layer_rend-DC_chamber_layer_rbegin)/DC_cell_size)"/>
+ <constant name="DC_sensitive_rmin" value="DC_chamber_layer_rbegin+0.5*(DC_chamber_layer_rend-DC_chamber_layer_rbegin-DC_cell_size*DC_layer_number)"/>
+ <constant name="DC_sensitive_rmax" value="DC_sensitive_rmin+DC_cell_size*DC_layer_number"/>
+ </define>
+
+ <detectors>
+ <detector id="DetID_DC" name="DriftChamber" type="StraightDriftChamber_v01" readout="DriftChamberHitsCollection" vis="DCVis" insideTrackingVolume="true" limits="dc_limits">
+ <envelope>
+ <shape type="BooleanShape" operation="Union" material="Air">
+ <shape type="Tube" rmin="DC_inner_radius" rmax="DC_outer_radius" dz="DC_half_length+DC_Endcap_dz" />
+ </shape>
+ </envelope>
+
+ <type_flags type="DetType_TRACKER + DetType_BARREL + DetType_GASEOUS + DetType_WIRE"/>
+
+ <shell>
+ <inner thickness="SDT_inner_wall_thickness" material="CarbonFiber" vis="LightGrayVis"/>
+ <outer thickness="SDT_outer_wall_thickness" material="CarbonFiber" vis="LightGrayVis"/>
+ <side thickness="0.1*mm" material="CarbonFiber" vis="LightGrayVis"/>
+ </shell>
+ <chamber rmin="DC_chamber_layer_rbegin" rmax="DC_chamber_layer_rend" zhalf="DC_half_length" material="GasHe_90Isob_10" vis="DCVis">
+ <layer id="0" name="MainChamber" rmin="DC_sensitive_rmin" rmax="DC_sensitive_rmax" number="DC_layer_number" material="GasHe_90Isob_10" vis="DCLayerVis"/>
+ </chamber>
+ </detector>
+ </detectors>
+
+ <readouts>
+ <readout name="DriftChamberHitsCollection">
+ <id>system:5,layer:7:9,chamber:8,cellID:32:16</id>
+ </readout>
+ </readouts>
+
+</lccdd>
diff --git a/Detector/DetCRD/src/Tracker/StraightDriftChamber_v01_geo.cpp b/Detector/DetCRD/src/Tracker/StraightDriftChamber_v01_geo.cpp
new file mode 100644
index 00000000..af878c09
--- /dev/null
+++ b/Detector/DetCRD/src/Tracker/StraightDriftChamber_v01_geo.cpp
@@ -0,0 +1,252 @@
+//====================================================================
+// Detector description implementation of the Drift Chamber with only straight
+//--------------------------------------------------------------------
+//
+// Author: Chengdong FU
+//
+//====================================================================
+
+#include <DD4hep/Detector.h>
+#include "DD4hep/DetFactoryHelper.h"
+#include "XML/Utilities.h"
+#include "DDRec/Surface.h"
+#include "DDRec/DetectorData.h"
+
+#include <map>
+
+using namespace std;
+
+using dd4hep::DetElement;
+using dd4hep::Detector;
+using dd4hep::Material;
+using dd4hep::PlacedVolume;
+using dd4hep::SensitiveDetector;
+using dd4hep::Volume;
+using dd4hep::Tube;
+using dd4hep::_toString;
+using dd4hep::Position;
+using dd4hep::rec::FixedPadSizeTPCData;
+using dd4hep::rec::Vector3D;
+using dd4hep::rec::VolCylinder;
+using dd4hep::rec::SurfaceType;
+using dd4hep::rec::volSurfaceList;
+using dd4hep::rec::VolPlane;
+
+static dd4hep::Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens) {
+ xml_det_t x_det = e;
+ Material air = description.air();
+ int det_id = x_det.id();
+ string name = x_det.nameStr();
+ DetElement tracker(name, det_id);
+
+ Volume envelope = dd4hep::xml::createPlacedEnvelope(description, e, tracker);
+ dd4hep::xml::setDetectorTypeFlag(e, tracker) ;
+ if(description.buildType()==dd4hep::BUILD_ENVELOPE) return tracker;
+ envelope.setVisAttributes(description.visAttributes("SeeThrough"));
+
+ sens.setType("tracker");
+ std::cout << " ** building StraightDriftChamber_v01 ..." << std::endl ;
+
+ xml_coll_t c_shell(x_det,_Unicode(shell));
+ xml_comp_t x_shell = c_shell;
+
+ xml_coll_t c_inner(x_shell,_U(inner));
+ xml_comp_t x_inner = c_inner;
+ double thickness_inner = x_inner.thickness();
+ Material mat_inner(description.material(x_inner.materialStr()));
+
+ xml_coll_t c_outer(x_shell,_U(outer));
+ xml_comp_t x_outer = c_outer;
+ double thickness_outer = x_outer.thickness();
+ Material mat_outer(description.material(x_outer.materialStr()));
+
+ xml_coll_t c_side(x_shell,_U(side));
+ xml_comp_t x_side = c_side;
+ double thickness_side = x_side.thickness();
+ Material mat_side(description.material(x_side.materialStr()));
+
+ xml_coll_t c_chamber(x_det,_U(chamber));
+ xml_comp_t x_chamber = c_chamber;
+ double rmin_chamber = x_chamber.rmin();
+ double rmax_chamber = x_chamber.rmax();
+ double halflength_chamber = x_chamber.zhalf();
+ Material mat_chamber(description.material(x_chamber.materialStr()));
+
+ xml_coll_t c_signal(x_det,_Unicode(signal_wire));
+ xml_comp_t x_signal = c_signal;
+ double rmin_signal = x_signal.rmin();
+ double rmax_signal = x_signal.rmax();
+ Material mat_signal = (description.material(x_signal.materialStr()));
+ Tube solid_signal(rmin_signal,rmax_signal,halflength_chamber);
+ Volume volume_signal(name+"_signal",solid_signal,mat_signal);
+ for(xml_coll_t c(x_signal, _U(tubs)); c; ++c) {
+ xml_comp_t x_tub = c;
+ string name_tub = x_tub.nameStr();
+ double rmin_tub = x_tub.rmin();
+ double rmax_tub = x_tub.rmax();
+ Material mat_tub(description.material(x_tub.materialStr()));
+ Tube solid_tub(rmin_tub,rmax_tub,halflength_chamber);
+ Volume volume_tub(name+"_signal"+name_tub,solid_tub,mat_tub);
+ volume_tub.setVisAttributes(description.visAttributes(x_tub.visStr()));
+ PlacedVolume phy_tub = volume_signal.placeVolume(volume_tub);
+ }
+
+ xml_coll_t c_field(x_det,_Unicode(field_wire));
+ xml_comp_t x_field = c_field;
+ double rmin_field = x_field.rmin();
+ double rmax_field = x_field.rmax();
+ Material mat_field = (description.material(x_field.materialStr()));
+ Tube solid_field(rmin_field,rmax_field,halflength_chamber);
+ Volume volume_field(name+"_field",solid_field,mat_field);
+ for(xml_coll_t c(x_field, _U(tubs)); c; ++c) {
+ xml_comp_t x_tub = c;
+ string name_tub = x_tub.nameStr();
+ double rmin_tub = x_tub.rmin();
+ double rmax_tub = x_tub.rmax();
+ Material mat_tub(description.material(x_tub.materialStr()));
+ Tube solid_tub(rmin_tub,rmax_tub,halflength_chamber);
+ Volume volume_tub(name+"_field"+name_tub,solid_tub,mat_tub);
+ volume_tub.setVisAttributes(description.visAttributes(x_tub.visStr()));
+ PlacedVolume phy_tub = volume_field.placeVolume(volume_tub);
+ }
+
+ Tube solid_chamber(rmin_chamber,rmax_chamber,halflength_chamber);
+ Volume volume_chamber(name+"_chamber",solid_chamber,mat_chamber);
+ volume_chamber.setVisAttributes(description.visAttributes(x_chamber.visStr()));
+ PlacedVolume phy_chamber = envelope.placeVolume(volume_chamber);
+ if(x_det.hasAttr(_U(id))){
+ phy_chamber.addPhysVolID("system",x_det.id());
+ }
+ DetElement det_chamber(tracker, name+"_chamber", 0);
+ det_chamber.setPlacement(phy_chamber);
+
+ Tube solid_inner(rmin_chamber-thickness_inner, rmin_chamber, halflength_chamber);
+ Volume volume_inner(name+"_inner_wall",solid_inner,mat_inner);
+ volume_inner.setVisAttributes(description.visAttributes(x_inner.visStr()));
+ PlacedVolume phy_inner = envelope.placeVolume(volume_inner);
+ DetElement det_inner(tracker, name+"_inner_wall", 0);
+ det_inner.setPlacement(phy_inner);
+ Vector3D ocyl_inner(rmin_chamber-0.5*thickness_inner, 0., 0.);
+ VolCylinder surfI(volume_inner, SurfaceType(SurfaceType::Helper), 0.5*thickness_inner, 0.5*thickness_inner, ocyl_inner);
+ volSurfaceList(tracker)->push_back(surfI);
+
+ Tube solid_outer(rmax_chamber, rmax_chamber+thickness_outer, halflength_chamber);
+ Volume volume_outer(name+"_outer_wall",solid_outer,mat_outer);
+ volume_outer.setVisAttributes(description.visAttributes(x_outer.visStr()));
+ PlacedVolume phy_outer = envelope.placeVolume(volume_outer);
+ DetElement det_outer(tracker, name+"_outer_wall", 0);
+ det_outer.setPlacement(phy_outer);
+ Vector3D ocyl_outer(rmax_chamber+0.5*thickness_inner, 0., 0.);
+ VolCylinder surfO(volume_outer, SurfaceType(SurfaceType::Helper), 0.5*thickness_outer, 0.5*thickness_outer, ocyl_outer);
+ volSurfaceList(tracker)->push_back(surfO);
+
+ Tube solid_side(rmin_chamber-thickness_inner, rmax_chamber+thickness_inner, 0.5*thickness_side);
+ Volume volume_side(name+"_side_wall",solid_side,mat_side);
+ volume_side.setVisAttributes(description.visAttributes(x_side.visStr()));
+ PlacedVolume phy_plus = envelope.placeVolume(volume_side,Position(0,0,halflength_chamber+0.5*thickness_side));
+ PlacedVolume phy_minus = envelope.placeVolume(volume_side,Position(0,0,-halflength_chamber-0.5*thickness_side));
+ DetElement det_plus(tracker, name+"_plusside_wall", 0);
+ det_plus.setPlacement(phy_plus);
+ DetElement det_minus(tracker, name+"_minusside_wall", 1);
+ det_minus.setPlacement(phy_minus);
+ Vector3D u(0., 1., 0.);
+ Vector3D v(1., 0., 0.);
+ Vector3D n(0., 0., 1.);
+ double mid_r = 0.5*(rmin_chamber-thickness_inner+rmax_chamber+thickness_inner);
+ Vector3D o(0., mid_r, 0.);
+ VolPlane surfS(volume_side, SurfaceType(SurfaceType::Helper), 0.5*thickness_side, 0.5*thickness_side, u, v, n, o);
+ volSurfaceList(det_plus)->push_back(surfS);
+ volSurfaceList(det_minus)->push_back(surfS);
+
+ int chamber_id=0, max_layer=0;
+ double rmin_sensitive=10000, rmax_sensitive=0, cell_height=0;
+ for(xml_coll_t c(x_chamber, _U(layer)); c; ++c) {
+ xml_comp_t x_layer = c;
+ string name_layer = x_layer.nameStr();
+ double rmin_layer = x_layer.rmin();
+ double rmax_layer = x_layer.rmax();
+ Material mat_layer(description.material(x_layer.materialStr()));
+ Tube solid_sub(rmin_layer, rmax_layer, halflength_chamber);
+ Volume volume_sub(name+name_layer,solid_sub,mat_layer);
+ volume_sub.setVisAttributes(description.visAttributes("SeeThrough"));
+ PlacedVolume phy_sub = volume_chamber.placeVolume(volume_sub);
+ DetElement det_sub(det_chamber, name+name_layer, x_layer.id());
+ det_sub.setPlacement(phy_sub);
+
+ int nlayer = x_layer.number();
+ double height = (rmax_layer-rmin_layer)/nlayer;
+ double radius = rmin_layer;
+ for(int layer_id=0; layer_id<nlayer; layer_id++){
+ Tube solid_layer(radius, radius+height, halflength_chamber);
+ Volume volume_layer(name+name_layer+_toString(layer_id, "_%d"),solid_layer,mat_layer);
+ volume_layer.setVisAttributes(description.visAttributes(x_layer.visStr()));
+
+ PlacedVolume phy = volume_sub.placeVolume(volume_layer);
+ if(x_layer.isSensitive()){
+ volume_layer.setSensitiveDetector(sens);
+ volume_layer.setLimitSet(description,x_det.limitsStr());
+ phy.addPhysVolID("chamber", chamber_id).addPhysVolID("layer", layer_id);
+ if(radius<rmin_sensitive){
+ rmin_sensitive = radius;
+ // TODO: more than one chamber, with different height; now only minimum radius chamber include
+ rmax_sensitive = radius + nlayer*height;
+ max_layer = nlayer;
+ cell_height = height;
+ }
+ double radius_center = radius+0.5*height;
+ double radius_edge = radius+rmax_field;
+ int ncell = floor(dd4hep::twopi*radius_center/height);
+ double dphi = dd4hep::twopi / ncell;
+ double offset = 0.;
+ if(layer_id%2!=0) offset = 0.5*dphi;
+ //std::cout << "debug: " << layer_id << " rmid=" << radius_center << " redge=" << radius_edge
+ // << " ncell=" << ncell << " dphi=" << dphi << " offset=" << offset << std::endl;
+ for(int icell=0;icell<ncell;icell++){
+ double phi = (icell+0.5)*dphi + offset;
+ volume_layer.placeVolume(volume_signal, Position(radius_center*cos(phi),radius_center*sin(phi),0));
+ volume_layer.placeVolume(volume_field, Position(radius_center*cos(phi+0.5*dphi),radius_center*sin(phi+0.5*dphi),0));
+ volume_layer.placeVolume(volume_field, Position(radius_edge*cos(phi+0.5*dphi),radius_edge*sin(phi+0.5*dphi),0));
+ volume_layer.placeVolume(volume_field, Position(radius_edge*cos(phi+0.25*dphi),radius_edge*sin(phi+0.25*dphi),0));
+ volume_layer.placeVolume(volume_field, Position(radius_edge*cos(phi),radius_edge*sin(phi),0));
+ volume_layer.placeVolume(volume_field, Position(radius_edge*cos(phi-0.25*dphi),radius_edge*sin(phi-0.25*dphi),0));
+
+ if(layer_id==nlayer){
+ double radius_max = radius+height-rmax_field;
+ volume_layer.placeVolume(volume_field, Position(radius_max*cos(phi+0.5*dphi),radius_max*sin(phi+0.5*dphi),0));
+ volume_layer.placeVolume(volume_field, Position(radius_max*cos(phi+0.25*dphi),radius_max*sin(phi+0.25*dphi),0));
+ volume_layer.placeVolume(volume_field, Position(radius_max*cos(phi),radius_max*sin(phi),0));
+ volume_layer.placeVolume(volume_field, Position(radius_max*cos(phi-0.25*dphi),radius_max*sin(phi-0.25*dphi),0));
+ }
+ }
+ }
+ DetElement det_layer(det_sub, name+name_layer+_toString(layer_id, "_%d"), layer_id);
+ det_layer.setPlacement(phy);
+ Vector3D ol(radius+0.5*height, 0., 0.);
+ SurfaceType type = x_layer.isSensitive()?SurfaceType(SurfaceType::Sensitive, SurfaceType::Invisible):SurfaceType(SurfaceType::Helper);
+ VolCylinder surf(volume_layer, type, 0.5*height, 0.5*height, ol);
+ volSurfaceList(det_layer)->push_back(surf);
+
+ radius += height;
+ }
+ if(x_layer.isSensitive()) chamber_id++;
+ }
+
+ FixedPadSizeTPCData* dcData = new FixedPadSizeTPCData;
+ dcData->zHalf = halflength_chamber+thickness_side;
+ dcData->rMin = rmin_chamber-thickness_inner;
+ dcData->rMax = rmax_chamber+thickness_outer;
+ dcData->innerWallThickness = thickness_inner;
+ dcData->outerWallThickness = thickness_outer;
+ dcData->rMinReadout = rmin_sensitive;
+ dcData->rMaxReadout = rmax_sensitive;
+ dcData->maxRow = max_layer;
+ dcData->padHeight = cell_height;
+ dcData->padWidth = cell_height;
+ dcData->driftLength = halflength_chamber;
+ dcData->zMinReadout = thickness_side/2.0;
+ tracker.addExtension<FixedPadSizeTPCData>(dcData);
+
+ return tracker;
+}
+
+DECLARE_DETELEMENT(StraightDriftChamber_v01, create_detector)
diff --git a/Digitisers/SimpleDigi/CMakeLists.txt b/Digitisers/SimpleDigi/CMakeLists.txt
index 4e2687bc..28c8ea89 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 00000000..8bc07eed
--- /dev/null
+++ b/Digitisers/SimpleDigi/src/CylinderDigiAlg.cpp
@@ -0,0 +1,120 @@
+#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.setEdx (simhit.getEDep()/simhit.getPathLength());
+ 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 00000000..af3ffe21
--- /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 fa8ba9ef..210e3639 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 43dfa178..f5ed24ed 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,18 @@ 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;
+ if(m_sdTypeOption==0){
+ DriftChamberSensitiveDetector* dcsd = new DriftChamberSensitiveDetector(name, *dd4hep_geo);
+ dcsd->setDedxSimTool(m_dedx_simtool);
+ sd = dcsd;
+ }
+ else if(m_sdTypeOption==1){
+ 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/DriftChamberSensDetTool.h b/Simulation/DetSimSD/src/DriftChamberSensDetTool.h
index e01445d5..537ff103 100644
--- a/Simulation/DetSimSD/src/DriftChamberSensDetTool.h
+++ b/Simulation/DetSimSD/src/DriftChamberSensDetTool.h
@@ -35,6 +35,7 @@ private:
SmartIF<IGeomSvc> m_geosvc;
ToolHandle<IDedxSimTool> m_dedx_simtool;
Gaudi::Property<std::string> m_dedx_sim_option{this, "DedxSimTool"};
+ Gaudi::Property<int> m_sdTypeOption{this, "TypeOption", 0};
};
diff --git a/Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.cpp b/Simulation/DetSimSD/src/TrackerCombineSensitiveDetector.cpp
new file mode 100644
index 00000000..0d53d3b9
--- /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 00000000..cf0cb844
--- /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
--
GitLab