diff --git a/Detector/DetCRD/CMakeLists.txt b/Detector/DetCRD/CMakeLists.txt index 6c242a58b96be073bfef767a8f5a56ad440ea048..1157f61a3abaf89748b76c1bebda4fe799a861a4 100644 --- a/Detector/DetCRD/CMakeLists.txt +++ b/Detector/DetCRD/CMakeLists.txt @@ -17,7 +17,10 @@ gaudi_add_module(DetCRD src/Muon/Muon_Barrel_v01_01.cpp src/Muon/Muon_Endcap_v01_02.cpp src/Tracker/SiTrackerSkewRing_v01_geo.cpp + src/Tracker/SiTrackerStitching_v01_geo.cpp src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp + src/Tracker/SiTrackerStaggeredLadder_v02_geo.cpp + src/Tracker/SiTrackerComposite_v01_geo.cpp src/Tracker/TPC_Simple_o1_v01.cpp src/Tracker/TPC_ModularEndcap_o1_v01.cpp src/Tracker/SiTracker_itkbarrel_v01_geo.cpp @@ -26,6 +29,7 @@ gaudi_add_module(DetCRD src/Tracker/SiTracker_otkendcap_v01_geo.cpp LINK ${DD4hep_COMPONENT_LIBRARIES} + Identifier ) set(LIBRARY_OUTPUT_PATH ${CMAKE_LIBRARY_OUTPUT_DIRECTORY}) diff --git a/Detector/DetCRD/compact/CRD_common_v02/VXD_Composite_v01_01.xml b/Detector/DetCRD/compact/CRD_common_v02/VXD_Composite_v01_01.xml new file mode 100644 index 0000000000000000000000000000000000000000..60ac647ebddd47153b65b0ed9aaf097c9750c6da --- /dev/null +++ b/Detector/DetCRD/compact/CRD_common_v02/VXD_Composite_v01_01.xml @@ -0,0 +1,165 @@ +<lccdd> + <info name="VXD_Stitching_v01_02" + title="CepC VXD with stitch module" + author="" + url="http://cepc.ihep.ac.cn" + status="developing" + version="v01"> + <comment>CepC vertex detector</comment> + </info> + <define> + <constant name="VXD_inner_radius" value="Vertex_inner_radius"/> + <constant name="VXD_outer_radius" value="Vertex_outer_radius"/> + <constant name="VXD_half_length" value="Vertex_half_length"/> + <constant name="VXDLayer1_length" value="260*mm" /> + <constant name="VXDLayer2_length" value="494*mm" /> + <constant name="VXDLayer3_length" value="749*mm" /> + <constant name="VXD_sensor_length" value="20*mm" /> + <constant name="VXD_sensor_xgap" value="0.01*mm"/> + <constant name="VXD_sensor_ygap" value="0.1*mm"/> + </define> + + <detectors> + <detector id="DetID_VXD" name="VXD" type="SiTrackerComposite_v01" vis="VXDVis" readout="VXDCollection" insideTrackingVolume="true" printLevel="INFO"> + <envelope> + <shape type="Assembly"/> + </envelope> + + <type_flags type="DetType_TRACKER + DetType_BARREL + DetType_PIXEL "/> + + <global sensitive_thickness="VXD_sensitive_thickness" support_thickness="VXD_support_thickness" sensor_length="VXD_sensor_length" + sensitive_mat="G4_Si" support_mat="G4_C" sensitive_threshold_KeV="64*keV" /> + <display ladder="SeeThrough" support="VXDSupportVis" flex="VXDFlexVis" sens_env="SeeThrough" sens="GrayVis" deadsensor="GreenVis" deadwire="RedVis"/> + + <shell rmin="70*mm" rmax="72.5*mm" zhalf="510*mm" material="CFRP_CMS" vis="LightGrayVis"/> + + <layer id="0" phi0="0" isBent="true"> + <module offset="0" phi="0" radius="12.2*mm" nx="10" ny="2" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="20*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="10.0*um" material="G4_Si"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module offset="0" phi="180*degree" radius="12.7*mm" nx="10" ny="2" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="20*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="10.0*um" material="G4_Si"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + </layer> + + <layer id="1" phi0="30*degree" isBent="true"> + <module offset="0" phi="0" radius="19.2*mm" nx="15" ny="3" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="21*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="10.0*um" material="G4_Si"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module offset="0" phi="180*degree" radius="19.7*mm" nx="15" ny="3" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="21*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="10.0*um" material="G4_Si"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + </layer> + + <layer id="2" phi0="60*degree" isBent="true"> + <module offset="0" phi="0" radius="25.9*mm" nx="20" ny="4" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="21*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="10.0*um" material="G4_Si"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module offset="0" phi="180*degree" radius="26.4*mm" nx="20" ny="4" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="21*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="10.0*um" material="G4_Si"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + </layer> + + <layer id="3" phi0="90*degree" isBent="true"> + <module offset="0" phi="0" radius="32.9*mm" nx="25" ny="5" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="22*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="10.0*um" material="G4_Si"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module offset="0" phi="180*degree" radius="33.4*mm" nx="25" ny="5" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="22*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="10.0*um" material="G4_Si"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + </layer> + + <layer id="4" ladder_radius="43.792*mm" ladder_offset="(8.7+11.7)*mm" n_ladders="25" n_sensors_per_side="0"> + <ladder isDoubleSided="true"> + <ladderSupport height="3.2*mm" length="VXDLayer3_length" thickness="370*um" width="17.4*mm" mat="CFRP_CMS"/> + <flex> + <slice length="VXDLayer3_length" thickness="12.5*um" width="17.4*mm" mat="Acrylicglue"/> <!--glue between flex and sensor/support--> + <slice length="VXDLayer3_length" thickness="12.5*um" width="17.4*mm" mat="Kapton"/> + <slice length="VXDLayer3_length" thickness="12.5*um" width="17.4*mm" mat="Acrylicglue"/> + <slice length="VXDLayer3_length" thickness=" 8.0*um" width="17.4*mm" mat="G4_Al"/> + <slice length="VXDLayer3_length" thickness="13.0*um" width="17.4*mm" mat="Kapton"/> + <slice length="VXDLayer3_length" thickness="12.5*um" width="17.4*mm" mat="Acrylicglue"/> + <slice length="VXDLayer3_length" thickness=" 8.0*um" width="17.4*mm" mat="G4_Al"/> + <slice length="VXDLayer3_length" thickness="13.0*um" width="17.4*mm" mat="Kapton"/> + <slice length="VXDLayer3_length" thickness="12.5*um" width="17.4*mm" mat="Acrylicglue"/> + <slice length="VXDLayer3_length" thickness="12.0*um" width="17.4*mm" mat="G4_Al"/> + <slice length="VXDLayer3_length" thickness="25.0*um" width="17.4*mm" mat="Kapton"/> + <slice length="VXDLayer3_length" thickness="12.0*um" width="17.4*mm" mat="G4_Al"/> + <slice length="VXDLayer3_length" thickness="12.5*um" width="17.4*mm" mat="Acrylicglue"/> + <slice length="VXDLayer3_length" thickness="13.0*um" width="17.4*mm" mat="Kapton"/> + <slice length="VXDLayer3_length" thickness=" 8.0*um" width="17.4*mm" mat="G4_Al"/> + <slice length="VXDLayer3_length" thickness="12.5*um" width="17.4*mm" mat="Acrylicglue"/> + <slice length="VXDLayer3_length" thickness="13.0*um" width="17.4*mm" mat="Kapton"/> + <slice length="VXDLayer3_length" thickness=" 8.0*um" width="17.4*mm" mat="G4_Al"/> + <slice length="VXDLayer3_length" thickness="12.5*um" width="17.4*mm" mat="Acrylicglue"/> + <slice length="VXDLayer3_length" thickness="12.5*um" width="17.4*mm" mat="Kapton"/> + <slice length="VXDLayer3_length" thickness="12.5*um" width="17.4*mm" mat="Acrylicglue"/> <!--glue between flex and sensor/support--> + </flex> + <sensor n_sensors="29" gap="0.1*mm" thickness="50*um" active_length="25.6*mm" active_width="12.8*mm" dead_width="2*mm" sensor_mat="G4_Si" + deadwire_length="(29*(25.6+0.1)-0.1)*mm" deadwire_width="2.6*mm" deadwire_thickness="(50/10)*um" deadwire_mat="G4_Al"/> + </ladder> + </layer> + </detector> + </detectors> + + <readouts> + <readout name="VXDCollection"> + <!--segmentation type="CartesianGridYZ" grid_size_y="0.016*mm" grid_size_z="0.016*mm"/> + <id>system:5,side:-2,layer:9,module:8,sensor:32:8,y:-12,z:-12</id--> + <!-- old tracking not use senor id: 24-31 bit--> + <id>system:5,side:-2,layer:9,module:8,sensor:32:8</id> + </readout> + </readouts> +</lccdd> diff --git a/Detector/DetCRD/compact/CRD_common_v02/VXD_Stitching_v01_01.xml b/Detector/DetCRD/compact/CRD_common_v02/VXD_Stitching_v01_01.xml new file mode 100644 index 0000000000000000000000000000000000000000..a84be6994fd023424f87ba75b8a3084c6c68c9ee --- /dev/null +++ b/Detector/DetCRD/compact/CRD_common_v02/VXD_Stitching_v01_01.xml @@ -0,0 +1,207 @@ +<lccdd> + <info name="VXD_Stitching_v01_01" + title="CepC VXD with stitch module" + author="" + url="http://cepc.ihep.ac.cn" + status="developing" + version="v01"> + <comment>CepC vertex detector</comment> + </info> + <define> + <constant name="VXD_inner_radius" value="Vertex_inner_radius"/> + <constant name="VXD_outer_radius" value="Vertex_outer_radius"/> + <constant name="VXD_half_length" value="Vertex_half_length"/> + <!--constant name="VXDLayer1_length" value="260*mm" /> + <constant name="VXDLayer2_length" value="494*mm" /> + <constant name="VXDLayer3_length" value="749*mm" /--> + <constant name="VXD_sensor_length" value="20*mm" /> + <constant name="VXD_sensor_xgap" value="0.01*mm"/> + <constant name="VXD_sensor_ygap" value="0.1*mm"/> + </define> + + <detectors> + <detector id="DetID_VXD" name="VXD" type="SiTrackerStitching_v01" vis="VXDVis" readout="VXDCollection" insideTrackingVolume="true"> + <envelope> + <shape type="Assembly"/> + </envelope> + + <type_flags type="DetType_TRACKER + DetType_BARREL + DetType_PIXEL "/> + + <shell rmin="70*mm" rmax="VXD_outer_radius" zhalf="VXD_half_length" material="CFRP_CMS" vis="LightGrayVis"/> + + <modules> + <module type="Module1A" radius="12.2*mm" nx="10" ny="2" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="20*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="12.5*um" material="Kapton"/> + <slice thickness="10.0*um" material="G4_Al"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module type="Module1B" radius="12.7*mm" nx="10" ny="2" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="20*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="12.5*um" material="Kapton"/> + <slice thickness="10.0*um" material="G4_Al"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module type="Module2A" radius="13.2*mm" nx="10" ny="2" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="21*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="12.5*um" material="Kapton"/> + <slice thickness="10.0*um" material="G4_Al"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module type="Module2B" radius="13.7*mm" nx="10" ny="2" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="21*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="12.5*um" material="Kapton"/> + <slice thickness="10.0*um" material="G4_Al"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module type="Module3A" radius="25.9*mm" nx="20" ny="4" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="21*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="12.5*um" material="Kapton"/> + <slice thickness="10.0*um" material="G4_Al"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module type="Module3B" radius="26.4*mm" nx="20" ny="4" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="21*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="12.5*um" material="Kapton"/> + <slice thickness="10.0*um" material="G4_Al"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module type="Module4A" radius="26.9*mm" nx="20" ny="4" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="22*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="12.5*um" material="Kapton"/> + <slice thickness="10.0*um" material="G4_Al"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module type="Module4B" radius="27.4*mm" nx="20" ny="4" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="22*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="12.5*um" material="Kapton"/> + <slice thickness="10.0*um" material="G4_Al"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module type="Module5A" radius="39.6*mm" nx="30" ny="6" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="22*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="12.5*um" material="Kapton"/> + <slice thickness="10.0*um" material="G4_Al"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module type="Module5B" radius="40.1*mm" nx="30" ny="6" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="22*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="12.5*um" material="Kapton"/> + <slice thickness="10.0*um" material="G4_Al"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module type="Module6A" radius="40.6*mm" nx="30" ny="6" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="23*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="12.5*um" material="Kapton"/> + <slice thickness="10.0*um" material="G4_Al"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + + <module type="Module6B" radius="41.1*mm" nx="30" ny="6" xdead="VXD_sensor_xgap" ydead="VXD_sensor_ygap" material="G4_Si"> + <sensor thickness="50*um" width="23*mm" length="20*mm" material="G4_Si" vis="VXDVis"/> + <flex vis="GrayVis"> + <slice thickness="12.5*um" material="Kapton"/> + <slice thickness="10.0*um" material="G4_Al"/> + </flex> + <electronics thickness="100*um" width="5*mm" material="Kapton" vis="RedVis"/> + <readout thickness="100*um" width="5*mm" material="Kapton" vis="GreenVis"/> + <driver thickness="100*um" width="8*mm" material="Kapton" vis="BlueVis"/> + </module> + </modules> + + <layer id="0" phi0="0"> + <module type="Module1A" offset="0" phi="0" z="0"/> + <module type="Module1B" offset="0" phi="180*degree" z="0"/> + </layer> + + <layer id="1" phi0="30*degree" length="0"> + <module type="Module2A" offset="0" phi="0" z="0"/> + <module type="Module2B" offset="0" phi="180*degree" z="0"/> + </layer> + + <layer id="2" phi0="60*degree" length="0"> + <module type="Module3A" offset="0" phi="0" z="0"/> + <module type="Module3B" offset="0" phi="180*degree" z="0"/> + </layer> + + <layer id="3" phi0="90*degree" length="0"> + <module type="Module4A" offset="0" phi="0" z="0"/> + <module type="Module4B" offset="0" phi="180*degree" z="0"/> + </layer> + + <layer id="4" phi0="120*degree" length="0"> + <module type="Module5A" offset="0" phi="0" z="0"/> + <module type="Module5B" offset="0" phi="180*degree" z="0"/> + </layer> + + <layer id="5" phi0="150*degree" length="0"> + <module type="Module6A" offset="0" phi="0" z="0"/> + <module type="Module6B" offset="0" phi="180*degree" z="0"/> + </layer> + + </detector> + + </detectors> + + <readouts> + <readout name="VXDCollection"> + <!--segmentation type="CartesianGridYZ" grid_size_y="0.016*mm" grid_size_z="0.016*mm"/> + <id>system:5,side:-2,layer:9,module:8,active:8,sensor:8,y:-12,z:-12</id--> + <id>system:5,side:-2,layer:9,module:8,active:8,sensor:8</id> + </readout> + </readouts> +</lccdd> diff --git a/Detector/DetCRD/src/Tracker/SiTrackerComposite_v01_geo.cpp b/Detector/DetCRD/src/Tracker/SiTrackerComposite_v01_geo.cpp new file mode 100644 index 0000000000000000000000000000000000000000..5bd1172b8c33390259988875282323440f584465 --- /dev/null +++ b/Detector/DetCRD/src/Tracker/SiTrackerComposite_v01_geo.cpp @@ -0,0 +1,554 @@ +//==================================================================== +// cepcgeo - CEPC silicon detector models in DD4hep +//-------------------------------------------------------------------- +// Chengdong FU, IHEP +// email: fucd@ihep.ac.cn +// $Id$ +//==================================================================== +#include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/DD4hepUnits.h" +#include "DD4hep/DetType.h" +#include "DD4hep/Printout.h" +#include "DDRec/Surface.h" +#include "DDRec/DetectorData.h" +#include "XML/Utilities.h" +#include <cmath> + +#include "Identifier/CEPCDetectorData.h" + +using namespace std; + +using dd4hep::Box; +using dd4hep::Tube; +using dd4hep::DetElement; +using dd4hep::Material; +using dd4hep::Position; +using dd4hep::RotationY; +using dd4hep::RotationZ; +using dd4hep::RotationZYX; +using dd4hep::Transform3D; +using dd4hep::Rotation3D; +using dd4hep::Volume; +using dd4hep::_toString; +using dd4hep::rec::volSurfaceList; +using dd4hep::rec::ZPlanarData; +using dd4hep::mm; + +/** Construction of the VXD detector, ported from Mokka driver SIT_Simple_Pixel.cc + * + * Mokka History: + * Feb 7th 2011, Steve Aplin - original version + * F.Gaede, DESY, Jan 2014 - dd4hep SIT pixel + * Hao Zeng, IHEP, July 2021 + * Chengdong FU, IHEP, Sep 2024 - composite from SiTrackerStaggeredLadder_v02 and SiTrackerStitching_v01 + */ + +static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4hep::SensitiveDetector sens) { + + xml_det_t x_det = e; + Material air = theDetector.air(); + int det_id = x_det.id(); + string name = x_det.nameStr(); + DetElement det(name, det_id); + + Volume envelope = dd4hep::xml::createPlacedEnvelope(theDetector, e, det); + dd4hep::xml::setDetectorTypeFlag(e, det) ; + if(theDetector.buildType()==dd4hep::BUILD_ENVELOPE) return det; + envelope.setVisAttributes(theDetector.visAttributes("SeeThrough")); + + dd4hep::PrintLevel oldLevel = dd4hep::printLevel(); + if (x_det.hasAttr(_Unicode(printLevel))) { + dd4hep::PrintLevel printLevel = dd4hep::printLevel(x_det.attr<string>(_Unicode(printLevel))); + // global level alway larger than local + if (printLevel<oldLevel) dd4hep::setPrintLevel(printLevel); + } + + if (x_det.hasChild(_U(sensitive))) { + xml_dim_t sd_typ = x_det.child(_U(sensitive)); + sens.setType(sd_typ.typeStr()); + } + else { + sens.setType("tracker"); + } + + dd4hep::printout(dd4hep::INFO, "Construct", "** building SiTrackerComposite_v01 ..."); + + dd4hep::rec::CompositeData* compositeData = new dd4hep::rec::CompositeData; + + //fetch the display parameters + xml_comp_t x_display(x_det.child(_Unicode(display))); + std::string ladderVis = x_display.attr<string>(_Unicode(ladder)); + std::string supportVis = x_display.attr<string>(_Unicode(support)); + std::string flexVis = x_display.attr<string>(_Unicode(flex)); + std::string sensEnvVis = x_display.attr<string>(_Unicode(sens_env)); + std::string sensVis = x_display.attr<string>(_Unicode(sens)); + std::string deadsensVis = x_display.attr<string>(_Unicode(deadsensor)); + std::string deadwireVis = x_display.attr<string>(_Unicode(deadwire)); + + //fetch the shell parameters + if (x_det.hasChild(_Unicode(shell))) { + xml_comp_t x_shell(x_det.child(_Unicode(shell))); + double rmin_shell = x_shell.rmin(); + double rmax_shell = x_shell.rmax(); + double zhalf_shell = x_shell.zhalf(); + Tube shellSolid(rmin_shell, rmax_shell, zhalf_shell); + Volume shellLogical(name + "_ShellLogical", shellSolid, theDetector.material(x_shell.materialStr())); + shellLogical.setVisAttributes(theDetector.visAttributes(x_shell.visStr())); + envelope.placeVolume(shellLogical); + + compositeData->zHalfShell = zhalf_shell; + compositeData->rInnerShell = rmin_shell; + compositeData->rOuterShell = rmax_shell; + } + + for(xml_coll_t layer_i(x_det,_U(layer)); layer_i; ++layer_i){ + xml_comp_t x_layer(layer_i); + + dd4hep::PlacedVolume pv; + int layer_id = x_layer.id(); + bool isBent = false; + if (x_layer.hasAttr(_Unicode(isBent))) isBent = x_layer.attr<bool>(_Unicode(isBent)); + + dd4hep::printout(dd4hep::INFO, "Construct", "layer_id: %02d --- %s", layer_id, isBent ? "Stitching" : "Planar"); + + double z = x_layer.hasAttr(_U(z)) ? x_layer.z() : 0; + + if (isBent) { + double phi0 = x_layer.phi0(); + + dd4hep::rec::CylindricalData::LayerLayout thisLayer; + + int module_id = 0; + for (xml_coll_t module_i(x_layer, _U(module)); module_i; ++module_i) { + if (module_id>1) dd4hep::printout(dd4hep::ERROR, "Construct", "Not support more than two modules, possible wrong in reconstruction"); + + string name_module = name + dd4hep::_toString(layer_id, "_Layer%02d") + dd4hep::_toString(module_id, "_Stave%02d"); + DetElement moduleDE(det, name_module, x_det.id()); + + xml_comp_t x_module(module_i); + + double offset = x_module.offset(); + double phi = x_module.phi(); + + double radius = x_module.radius(); + double xdead = x_module.attr<double>(_Unicode(xdead)); + double ydead = x_module.attr<double>(_Unicode(ydead)); + int nx = x_module.attr<int>(_Unicode(nx)); + int ny = x_module.attr<int>(_Unicode(ny)); + Material mat = theDetector.material(x_module.materialStr()); + + xml_comp_t x_sensor(x_module.child(_U(sensor))); + double sensor_thickness = x_sensor.thickness(); + double sensor_width = x_sensor.width(); + double sensor_length = x_sensor.length(); + double sensor_dphi = sensor_width/radius; + Material sensor_mat = theDetector.material(x_sensor.materialStr()); + + Tube sensor_solid(radius, radius+sensor_thickness, sensor_length/2, -sensor_width/radius/2, sensor_width/radius/2); + Volume sensor_volume(name_module + "Sensor", sensor_solid, sensor_mat); + sensor_volume.setSensitiveDetector(sens); + if (x_det.hasAttr(_U(limits))) sensor_volume.setLimitSet(theDetector, x_det.limitsStr()); + sensor_volume.setVisAttributes(theDetector.visAttributes(x_sensor.visStr())); + + xml_comp_t x_flex(x_module.child(_Unicode(flex))); + std::vector<std::pair<double, Material> > flexs; + double flex_thickness = 0; + for (xml_coll_t slice_i(x_flex, _U(slice)); slice_i; ++slice_i) { + xml_comp_t x_slice(slice_i); + double thickness = x_slice.thickness(); + Material mat_slice = theDetector.material(x_slice.materialStr()); + flexs.push_back(std::make_pair(thickness, mat_slice)); + flex_thickness += thickness; + dd4hep::printout(dd4hep::DEBUG, "Construct", "flex %s: %f mm", mat_slice.name(), thickness/dd4hep::mm); + } + double flex_width = sensor_width*ny + ydead*(ny-1); + double flex_length = sensor_length*nx + xdead*(nx-1); + double flex_radius = radius+sensor_thickness; + double flex_dphi = flex_width/radius; + Tube flex_solid(flex_radius, flex_radius+flex_thickness, flex_length/2, -flex_dphi/2, flex_dphi/2); + Volume flex_volume(name_module + "Flex", flex_solid, air); + flex_volume.setVisAttributes(theDetector.visAttributes(x_flex.visStr())); + double start_radius = flex_radius; + for (unsigned islice=0; islice<flexs.size(); islice++) { + dd4hep::printout(dd4hep::DEBUG, "Construct", "flex start radius: %f mm", start_radius/dd4hep::mm); + Tube slice_solid(start_radius, start_radius+flexs[islice].first, flex_length/2, -flex_dphi/2, flex_dphi/2); + Volume slice_volume(name_module + dd4hep::_toString(int(islice), "Flex_%d"), slice_solid, flexs[islice].second); + flex_volume.placeVolume(slice_volume); + start_radius += flexs[islice].first; + } + + double service_radius = radius + sensor_thickness; + + xml_comp_t x_electronics(x_module.child(_Unicode(electronics))); + double electronics_thickness = x_electronics.thickness(); + double electronics_width = x_electronics.width(); + double electronics_dphi = electronics_width/radius; + Material electronics_mat = theDetector.material(x_electronics.materialStr()); + Tube electronics_solid(service_radius, service_radius+electronics_thickness, flex_length/2, -electronics_dphi, 0); + Volume electronics_volume(name_module + "Electronics", electronics_solid, electronics_mat); + electronics_volume.setVisAttributes(theDetector.visAttributes(x_electronics.visStr())); + + xml_comp_t x_readout(x_module.child(_U(readout))); + double readout_thickness = x_readout.thickness(); + double readout_width = x_readout.width(); + double readout_dphi = readout_width/radius; + Material readout_mat = theDetector.material(x_readout.materialStr()); + Tube readout_solid(service_radius, service_radius+readout_thickness, flex_length/2, 0, readout_dphi); + Volume readout_volume(name_module + "Electronics", readout_solid, readout_mat); + readout_volume.setVisAttributes(theDetector.visAttributes(x_readout.visStr())); + + xml_comp_t x_driver(x_module.child(_Unicode(driver))); + double driver_thickness = x_driver.thickness(); + // width -> length, length->width + double driver_length = x_driver.width(); + double driver_width = flex_width + electronics_width + readout_width; + double driver_dphi = driver_width/radius; + Material driver_mat = theDetector.material(x_driver.materialStr()); + Tube driver_solid(service_radius, service_radius+driver_thickness, driver_length/2, -flex_dphi/2-electronics_dphi, flex_dphi/2+readout_dphi); + Volume driver_volume(name_module + "Driver", driver_solid, driver_mat); + driver_volume.setVisAttributes(theDetector.visAttributes(x_driver.visStr())); + + double module_length = flex_length + 2*driver_length; + double module_width = driver_width; + double module_dphi = module_width/radius; + double module_thickness = sensor_thickness + std::max(std::max(std::max(flex_thickness, driver_thickness), readout_thickness), electronics_thickness); + Tube module_solid(radius, radius+module_thickness, module_length/2, -flex_dphi/2-electronics_dphi, flex_dphi/2+readout_dphi); + Volume module_volume(name_module, module_solid, air); + module_volume.setVisAttributes(theDetector.visAttributes("SeeThrough")); + + Tube board_solid(radius, radius+sensor_thickness, module_length/2, -flex_dphi/2-electronics_dphi, flex_dphi/2+readout_dphi); + Volume board_volume(name_module + "Board", board_solid, mat); + board_volume.setVisAttributes(theDetector.visAttributes(x_module.visStr())); + module_volume.placeVolume(board_volume); + + for (int ix=0; ix<nx; ix++) { + double z = -flex_length/2+sensor_length/2+ix*(sensor_length+xdead); + for (int iy=0; iy<ny; iy++) { + double delta = -flex_dphi/2+sensor_dphi/2+iy*(sensor_dphi+ydead/radius); + Transform3D tran(RotationZ(delta), Position(0, 0, z)); + dd4hep::PlacedVolume pv = board_volume.placeVolume(sensor_volume, tran); + + int sensor_id = ix + nx*iy; + pv.addPhysVolID("sensor", sensor_id); + + dd4hep::rec::Vector3D ocyl(radius + 0.5*sensor_thickness, 0., 0.); + dd4hep::rec::VolCylinder surf(sensor_volume, dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Sensitive), + 0.5*sensor_thickness, module_thickness-0.5*sensor_thickness, ocyl); + DetElement sensorDE(moduleDE, name_module + dd4hep::_toString(sensor_id, "_Sensor%02d"), x_det.id()); + sensorDE.setPlacement(pv); + volSurfaceList(sensorDE)->push_back(surf); + } + } + + module_volume.placeVolume(flex_volume); + module_volume.placeVolume(electronics_volume, RotationZYX(-flex_dphi/2,0,0)); + module_volume.placeVolume(readout_volume, RotationZYX(flex_dphi/2,0,0)); + module_volume.placeVolume(driver_volume, Position(0, 0, -flex_length/2-driver_length/2)); + module_volume.placeVolume(driver_volume, Position(0, 0, flex_length/2+driver_length/2)); + + if (module_id == 0) { + thisLayer.zHalf = flex_length/2.0; + thisLayer.radiusSensitive = radius; + thisLayer.thicknessSensitive = sensor_thickness; + thisLayer.radiusSupport = flex_radius; + thisLayer.thicknessSupport = flex_thickness; + //thisLayer.width = flex_width; + thisLayer.phi0 = phi0 + phi; + thisLayer.rgap = radius; + } + else { + thisLayer.rgap = radius - thisLayer.rgap; + thisLayer.dphi = phi0 + phi - thisLayer.phi0; + } + + Transform3D tran(RotationZ(phi0+phi), Position(offset*cos(phi), offset*sin(phi), z)); + pv = envelope.placeVolume(module_volume, tran); + pv.addPhysVolID("layer", layer_id).addPhysVolID("module", module_id); + + moduleDE.setPlacement(pv); + + module_id++; + } + + compositeData->layersBent.push_back(thisLayer); + } + else { + double sensitive_radius = x_layer.attr<double>(_Unicode(ladder_radius)); + int n_sensors_per_ladder = x_layer.attr<int>(_Unicode(n_sensors_per_side)); + int n_ladders = x_layer.attr<int>(_Unicode(n_ladders)) ; + double ladder_offset = x_layer.attr<double>(_Unicode(ladder_offset)); + double ladder_radius = sqrt(ladder_offset*ladder_offset + sensitive_radius*sensitive_radius); + double ladder_phi0 = -atan(ladder_offset/sensitive_radius); + dd4hep::printout(dd4hep::INFO, "Construct", "ladder_radius = %f mm, sensitive_radius = %f mm, n_sensors_per_ladder = %d", ladder_radius/mm, sensitive_radius/mm, n_sensors_per_ladder); + + std::string layerName = dd4hep::_toString(layer_id , "layer_%02d"); + dd4hep::Assembly layer_assembly(layerName); + pv = envelope.placeVolume(layer_assembly); + dd4hep::DetElement layerDE(det, layerName, x_det.id()); + layerDE.setPlacement(pv); + + const double ladder_dphi = ( dd4hep::twopi / n_ladders ) ; + dd4hep::printout(dd4hep::DEBUG, "Construct", "ladder_dphi: %f (%f degree)", ladder_dphi, ladder_dphi/dd4hep::degree); + + //fetch the ladder parameters + xml_comp_t x_ladder(x_layer.child(_Unicode(ladder))); + + //fetch the ladder support parameters + xml_comp_t x_ladder_support(x_ladder.child(_Unicode(ladderSupport))); + double support_length = x_ladder_support.attr<double>(_Unicode(length)); + double support_thickness = x_ladder_support.attr<double>(_Unicode(thickness)); + double support_height = x_ladder_support.attr<double>(_Unicode(height)); + double support_width = x_ladder_support.attr<double>(_Unicode(width)); + Material support_mat; + if(x_ladder_support.hasAttr(_Unicode(mat))) { + support_mat = theDetector.material(x_ladder_support.attr<string>(_Unicode(mat))); + } + else { + support_mat = theDetector.material(x_ladder_support.materialStr()); + } + dd4hep::printout(dd4hep::INFO, "Construct", "support_length = %f mm, support_thickness = %f mm, support_width = %f mm", support_length/mm, support_thickness/mm, support_width/mm); + + //fetch the flex parameters + double flex_thickness(0); + double flex_width(0); + double flex_length(0); + xml_comp_t x_flex(x_ladder.child(_Unicode(flex))); + for(xml_coll_t flex_i(x_flex,_U(slice)); flex_i; ++flex_i) { + xml_comp_t x_flex_slice(flex_i); + double x_flex_slice_thickness = x_flex_slice.attr<double>(_Unicode(thickness)); + double x_flex_slice_width = x_flex_slice.attr<double>(_Unicode(width)); + double x_flex_slice_length = x_flex_slice.attr<double>(_Unicode(length)); + flex_thickness += x_flex_slice_thickness; + if (x_flex_slice_width > flex_width) flex_width = x_flex_slice_width; + if (x_flex_slice_length > flex_length) flex_length = x_flex_slice_length; + dd4hep::printout(dd4hep::DEBUG, "Construct", "x_flex_slice_thickness: %f mm", x_flex_slice_thickness/mm); + } + dd4hep::printout(dd4hep::INFO, "Construct", "flex_length = %f mm, flex_thickness = %f mm, flex_width = %f mm", flex_length/mm, flex_thickness/mm, flex_width/mm); + + //fetch the sensor parameters + xml_comp_t x_sensor(x_ladder.child(_Unicode(sensor))); + int n_sensors_per_side = x_sensor.attr<int>(_Unicode(n_sensors)); + double dead_gap = x_sensor.attr<double>(_Unicode(gap)); + double sensor_thickness = x_sensor.attr<double>(_Unicode(thickness)); + double sensor_active_len = x_sensor.attr<double>(_Unicode(active_length)); + double sensor_active_width = x_sensor.attr<double>(_Unicode(active_width)); + double sensor_dead_width = x_sensor.attr<double>(_Unicode(dead_width)); + double sensor_deadwire_length = x_sensor.attr<double>(_Unicode(deadwire_length)); + double sensor_deadwire_width = x_sensor.attr<double>(_Unicode(deadwire_width)); + double sensor_deadwire_thickness = x_sensor.attr<double>(_Unicode(deadwire_thickness)); + Material sensor_mat = theDetector.material(x_sensor.attr<string>(_Unicode(sensor_mat))); + Material sensor_deadwire_mat = theDetector.material(x_sensor.attr<string>(_Unicode(deadwire_mat))); + + dd4hep::printout(dd4hep::INFO, "Construct", "sensor_active_length = %f mm, sensor_thickness = %f mm, sensor_active_width = %f mm", sensor_active_len/mm, sensor_thickness/mm, sensor_active_width/mm); + dd4hep::printout(dd4hep::INFO, "Construct", "sensor_dead_width = %f mm, dead_gap = %f mm, n_sensors_per_side = %d", sensor_dead_width/mm, dead_gap/mm, n_sensors_per_side); + + //create ladder logical volume + Box LadderSolid((support_height+2*sensor_thickness+2*flex_thickness)/2.0, support_width/2.0, support_length/2.0); + Volume LadderLogical(name + dd4hep::_toString(layer_id, "_LadderLogical_%02d"), LadderSolid, air); + // create flex envelope logical volume + Box FlexEnvelopeSolid(flex_thickness/2.0, flex_width/2.0, flex_length/2.0); + Volume FlexEnvelopeLogical(name + dd4hep::_toString( layer_id, "_FlexEnvelopeLogical_%02d"), FlexEnvelopeSolid, air); + FlexEnvelopeLogical.setVisAttributes(theDetector.visAttributes("SeeThrough")); + + //create the flex layers inside the flex envelope + double flex_start_height(-flex_thickness/2.); + int index = 0; + for(xml_coll_t flex_i(x_flex,_U(slice)); flex_i; ++flex_i){ + xml_comp_t x_flex_slice(flex_i); + double x_flex_slice_thickness = x_flex_slice.attr<double>(_Unicode(thickness)); + double x_flex_slice_width = x_flex_slice.attr<double>(_Unicode(width)); + double x_flex_slice_length = x_flex_slice.attr<double>(_Unicode(length)); + Material x_flex_slice_mat; + if(x_flex_slice.hasAttr(_Unicode(mat))) { + x_flex_slice_mat = theDetector.material(x_flex_slice.attr<string>(_Unicode(mat))); + } + else { + x_flex_slice_mat = theDetector.material(x_flex_slice.materialStr()); + } + // Material x_flex_slice_mat = theDetector.material(x_flex_slice.attr<string>(_Unicode(mat))); + Box FlexLayerSolid(x_flex_slice_thickness/2.0, x_flex_slice_width/2.0, x_flex_slice_length/2.0); + Volume FlexLayerLogical(name + dd4hep::_toString( layer_id, "_FlexLayerLogical_%02d") + dd4hep::_toString( index, "index_%02d"), FlexLayerSolid, x_flex_slice_mat); + FlexLayerLogical.setVisAttributes(theDetector.visAttributes(flexVis)); + double flex_slice_height = flex_start_height + x_flex_slice_thickness/2.; + pv = FlexEnvelopeLogical.placeVolume(FlexLayerLogical, Position(flex_slice_height, 0., 0.)); + flex_start_height += x_flex_slice_thickness; + index++; + } + + //place the flex envelope inside the ladder envelope + pv = LadderLogical.placeVolume(FlexEnvelopeLogical, Position((support_height + flex_thickness)/2.0, 0., 0.)); //top side + //define the transformation3D(only need a combination of translation and rotation) + Transform3D tran_mirro(RotationZYX(0., dd4hep::twopi/2.0, 0.), Position(-(support_height + flex_thickness)/2.0, 0., 0.)); + pv = LadderLogical.placeVolume(FlexEnvelopeLogical, tran_mirro); //bottom side + + //create sensor envelope logical volume + Box SensorTopEnvelopeSolid(sensor_thickness/2.0, support_width/2.0, support_length/2.0); + Volume SensorTopEnvelopeLogical(name + dd4hep::_toString( layer_id, "_SensorEnvelopeLogical_%02d"), SensorTopEnvelopeSolid, air); + Box SensorBottomEnvelopeSolid(sensor_thickness/2.0, support_width/2.0, support_length/2.0); + Volume SensorBottomEnvelopeLogical(name + dd4hep::_toString(layer_id, "_SensorEnvelopeLogical_%02d"), SensorBottomEnvelopeSolid, air); + SensorTopEnvelopeLogical.setVisAttributes(theDetector.visAttributes(sensEnvVis)); + + //create sensor logical volume + Box SensorSolid(sensor_thickness/2.0, sensor_active_width/2.0, sensor_active_len/2.0); + Volume SensorTopLogical(name + dd4hep::_toString(layer_id+1, "_SensorLogical_%02d"), SensorSolid, sensor_mat); + Volume SensorBottomLogical(name + dd4hep::_toString(layer_id, "_SensorLogical_%02d"), SensorSolid, sensor_mat); + SensorTopLogical.setSensitiveDetector(sens); + SensorBottomLogical.setSensitiveDetector(sens); + if (x_det.hasAttr(_U(limits))) { + SensorTopLogical.setLimitSet(theDetector, x_det.limitsStr()); + SensorBottomLogical.setLimitSet(theDetector, x_det.limitsStr()); + } + SensorTopLogical.setVisAttributes(theDetector.visAttributes(sensVis)); + SensorBottomLogical.setVisAttributes(theDetector.visAttributes(sensVis)); + + //create dead sensor logical volume + Box SensorDeadSolid(sensor_thickness / 2.0, sensor_dead_width / 2.0, sensor_active_len / 2.0); + Volume SensorDeadLogical(name + dd4hep::_toString( layer_id, "_SensorDeadLogical_%02d"), SensorDeadSolid, sensor_mat); + SensorDeadLogical.setVisAttributes(theDetector.visAttributes(deadsensVis)); + + //create dead wire logical volume + Box SensorDeadWireSolid(sensor_deadwire_thickness / 2.0, sensor_deadwire_width / 2.0, sensor_deadwire_length / 2.0); + Volume SensorDeadWireLogical(name + dd4hep::_toString( layer_id, "_SensorDeadWireLogical_%02d"), SensorDeadWireSolid, sensor_deadwire_mat); + SensorDeadWireLogical.setVisAttributes(theDetector.visAttributes(deadwireVis)); + + //place the dead wire in the sensor envelope + // pv = SensorTopEnvelopeLogical.placeVolume(SensorDeadWireLogical, Position(0.0, (sensor_active_width-support_width/2.0) + sensor_dead_width/2.0 + sensor_deadwire_width/2.0, 0.0)); + // pv = SensorBottomEnvelopeLogical.placeVolume(SensorDeadWireLogical, Position(0.0, (sensor_active_width-support_width/2.0) + sensor_dead_width/2.0 + sensor_deadwire_width/2.0, 0.0)); + pv = SensorTopEnvelopeLogical.placeVolume(SensorDeadWireLogical, Position(0.0, (-support_width/2.0) + (sensor_deadwire_width/2.0), 0.0)); + pv = SensorBottomEnvelopeLogical.placeVolume(SensorDeadWireLogical, Position(0.0, (-support_width/2.0) + (sensor_deadwire_width/2.0), 0.0)); + + // place the active sensor and dead sensor inside the sensor envelope + std::vector<dd4hep::PlacedVolume> TopSensor_pv; + std::vector<dd4hep::PlacedVolume> BottomSensor_pv; + for(int isensor=0; isensor < n_sensors_per_side; ++isensor){ + double sensor_total_z = n_sensors_per_side*sensor_active_len + dead_gap*(n_sensors_per_side-1); + double xpos = 0.0; + double ypos_active = (support_width/2.0) - (sensor_active_width/2.0); + double ypos_dead = (-support_width/2.0) + sensor_deadwire_width + (sensor_dead_width/2.0); + double zpos = -sensor_total_z/2.0 + sensor_active_len/2.0 + isensor*(sensor_active_len + dead_gap); + pv = SensorTopEnvelopeLogical.placeVolume(SensorTopLogical, Position(xpos,ypos_active,zpos)); + //pv.addPhysVolID("topsensor", isensor ) ; + pv.addPhysVolID("layer", layer_id+1).addPhysVolID("sensor", isensor); + TopSensor_pv.push_back(pv); + pv = SensorBottomEnvelopeLogical.placeVolume(SensorBottomLogical, Position(xpos,ypos_active,zpos)); + //pv.addPhysVolID("bottomsensor", isensor ) ; + pv.addPhysVolID("layer", layer_id ).addPhysVolID("sensor", isensor); + BottomSensor_pv.push_back(pv); + pv = SensorTopEnvelopeLogical.placeVolume(SensorDeadLogical, Position(xpos,ypos_dead,zpos)); + pv = SensorBottomEnvelopeLogical.placeVolume(SensorDeadLogical, Position(xpos,ypos_dead,zpos)); + } + //place the sensor envelope inside the ladder envelope + pv = LadderLogical.placeVolume(SensorTopEnvelopeLogical, + Position(support_height/2.0 + flex_thickness + sensor_thickness/2.0, 0., 0.));//top-side sensors + Position pos(-(support_height/2.0 + flex_thickness + sensor_thickness/2.0), 0., 0.); + pv = LadderLogical.placeVolume(SensorBottomEnvelopeLogical, pos);//bottom-side sensors + + //create the ladder support + Box LadderSupportSolid(support_height/2.0, support_width/2.0, support_length/2.0); + Volume LadderSupportLogical(name + _toString( layer_id,"_SupLogical_%02d"), LadderSupportSolid, support_mat); + LadderSupportLogical.setVisAttributes(theDetector.visAttributes(supportVis)); + + //create ladder support cavity + Box LadderSupportCavitySolid(support_height/2.0-support_thickness/2.0, support_width/2.0-support_thickness/2.0, support_length/2.0); + Volume LadderSupportCavityLogical(name + _toString( layer_id,"_SupCavityLogical_%02d"), LadderSupportCavitySolid, air); + LadderSupportCavityLogical.setVisAttributes(theDetector.visAttributes("SeeThrough")); + + pv = LadderSupportLogical.placeVolume(LadderSupportCavityLogical); + pv = LadderLogical.placeVolume(LadderSupportLogical); + + for(int i = 0; i < n_ladders; i++){ + std::stringstream ladder_enum; + ladder_enum << "vxt_ladder_" << layer_id << "_" << i; + DetElement ladderDE(layerDE, ladder_enum.str(), x_det.id()); + dd4hep::printout(dd4hep::DEBUG, "Construct", "start building %s", ladder_enum.str().c_str()); + + //====== create the meassurement surface =================== + dd4hep::rec::Vector3D o(0,0,0); + dd4hep::rec::Vector3D u( 0., 1., 0.); + dd4hep::rec::Vector3D v( 0., 0., 1.); + dd4hep::rec::Vector3D n( 1., 0., 0.); + // fucd: SensorLogical only sensor_thickness, support need another surface, todo + double inner_thick_top = sensor_thickness/2.0; + double outer_thick_top = sensor_thickness/2.0;//support_height/2.0 + flex_thickness + sensor_thickness/2.0; + double inner_thick_bottom = sensor_thickness/2.0;//support_height/2.0 + flex_thickness + sensor_thickness/2.0; + double outer_thick_bottom = sensor_thickness/2.0; + dd4hep::rec::VolPlane surfTop(SensorTopLogical, + dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Sensitive), + inner_thick_top, outer_thick_top, u, v, n, o); + dd4hep::rec::VolPlane surfBottom(SensorBottomLogical, + dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Sensitive), + inner_thick_bottom, outer_thick_bottom, u, v, n, o); + + for(int isensor=0; isensor < n_sensors_per_side; ++isensor){ + std::stringstream topsensor_str; + std::stringstream bottomsensor_str; + topsensor_str << ladder_enum.str() << "_top_" << isensor; + // std::cout << "\tstart building " << topsensor_str.str() << ":" << endl; + bottomsensor_str << ladder_enum.str() << "_bottom_" << isensor; + // std::cout << "\tstart building " << bottomsensor_str.str() << ":" << endl; + DetElement topsensorDE(ladderDE, topsensor_str.str(), x_det.id()); + DetElement bottomsensorDE(ladderDE, bottomsensor_str.str(), x_det.id()); + topsensorDE.setPlacement(TopSensor_pv[isensor]); + volSurfaceList(topsensorDE)->push_back(surfTop); + // std::cout << "\t" << topsensor_str.str() << " done." << endl; + bottomsensorDE.setPlacement(BottomSensor_pv[isensor]); + // std::cout << "\t" << bottomsensor_str.str() << " done." << endl; + volSurfaceList(bottomsensorDE)->push_back(surfBottom); + } + Transform3D tr (RotationZYX(ladder_dphi*i,0.,0.),Position(ladder_radius*cos(ladder_phi0+ladder_dphi*i), ladder_radius*sin(ladder_phi0+ladder_dphi*i), 0.)); + pv = layer_assembly.placeVolume(LadderLogical,tr); + pv.addPhysVolID("module", i ) ; + ladderDE.setPlacement(pv); + dd4hep::printout(dd4hep::DEBUG, "Construct", "%s done.", ladder_enum.str().c_str()); + if(i==0) dd4hep::printout(dd4hep::DEBUG, "Construct", "(x,y) = (%f, %f)mm", ladder_radius*cos(ladder_phi0)/mm, ladder_radius*sin(ladder_phi0)/mm); + } + + // package the reconstruction data + dd4hep::rec::ZPlanarData::LayerLayout topLayer; + dd4hep::rec::ZPlanarData::LayerLayout bottomLayer; + + topLayer.ladderNumber = n_ladders; + topLayer.phi0 = 0.; + topLayer.sensorsPerLadder = n_sensors_per_side; + topLayer.lengthSensor = sensor_active_len; + topLayer.distanceSupport = sensitive_radius; + topLayer.thicknessSupport = support_thickness / 2.0; + topLayer.offsetSupport = -ladder_offset; + topLayer.widthSupport = support_width; + topLayer.zHalfSupport = support_length / 2.0; + topLayer.distanceSensitive = sensitive_radius + support_height / 2.0 + flex_thickness; + topLayer.thicknessSensitive = sensor_thickness; + topLayer.offsetSensitive = -ladder_offset + (support_width/2.0 - sensor_active_width/2.0); + topLayer.widthSensitive = sensor_active_width; + topLayer.zHalfSensitive = (n_sensors_per_side*(sensor_active_len + dead_gap) - dead_gap) / 2.0; + + bottomLayer.ladderNumber = n_ladders; + bottomLayer.phi0 = 0.; + bottomLayer.sensorsPerLadder = n_sensors_per_side; + bottomLayer.lengthSensor = sensor_active_len; + bottomLayer.distanceSupport = sensitive_radius - support_height / 2.0 - flex_thickness; + bottomLayer.thicknessSupport = support_thickness / 2.0; + bottomLayer.offsetSupport = -ladder_offset; + bottomLayer.widthSupport = support_width; + bottomLayer.zHalfSupport = support_length / 2.0; + bottomLayer.distanceSensitive = sensitive_radius - support_height / 2.0 - sensor_thickness - flex_thickness; + bottomLayer.thicknessSensitive = sensor_thickness; + bottomLayer.offsetSensitive = -ladder_offset + (support_width/2.0 - sensor_active_width/2.0); + bottomLayer.widthSensitive = sensor_active_width; + bottomLayer.zHalfSensitive = (n_sensors_per_side*(sensor_active_len + dead_gap) - dead_gap) / 2.0; + + compositeData->layersPlanar.push_back(bottomLayer); + compositeData->layersPlanar.push_back(topLayer); + } + } + if (dd4hep::printLevel()<=dd4hep::INFO) std::cout << (*compositeData) << endl; + det.addExtension<dd4hep::rec::CompositeData>(compositeData); + + if (x_det.hasAttr(_U(combineHits))) det.setCombineHits(x_det.attr<bool>(_U(combineHits)),sens); + + dd4hep::printout(dd4hep::INFO, "Construct", "SiTrackerComposite_v01 done."); + dd4hep::setPrintLevel(oldLevel); + + return det; +} +DECLARE_DETELEMENT(SiTrackerComposite_v01,create_element) diff --git a/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp b/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp index 64fe00bbf274765aac5d35ed2dba86edc40151a6..b00c38536b373830ff0a172170b1815e54912a7f 100644 --- a/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp +++ b/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp @@ -84,7 +84,7 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h else { sens.setType("tracker"); } - std::cout << " ** building SiTrackerSkewBarrel_v01 ... " << sens.type() << std::endl ; + std::cout << " ** building SiTrackerStaggeredLadder_v01 ... " << sens.type() << std::endl ; dd4hep::rec::ZPlanarData* zPlanarData = new dd4hep::rec::ZPlanarData; @@ -264,10 +264,16 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h //create sensor logical volume Box SensorSolid(sensor_thickness / 2.0, sensor_active_width / 2.0, sensor_active_len / 2.0); - Volume SensorLogical(name + dd4hep::_toString( layer_id, "_SensorLogical_%02d"), SensorSolid, sensor_mat); - SensorLogical.setSensitiveDetector(sens); - if (x_det.hasAttr(_U(limits))) SensorLogical.setLimitSet(theDetector, x_det.limitsStr()); - SensorLogical.setVisAttributes(theDetector.visAttributes(sensVis)); + Volume SensorTopLogical(name + dd4hep::_toString(layer_id*2+1, "_SensorLogical_%02d"), SensorSolid, sensor_mat); + Volume SensorBottomLogical(name + dd4hep::_toString(layer_id*2, "_SensorLogical_%02d"), SensorSolid, sensor_mat); + SensorTopLogical.setSensitiveDetector(sens); + SensorBottomLogical.setSensitiveDetector(sens); + if (x_det.hasAttr(_U(limits))) { + SensorTopLogical.setLimitSet(theDetector, x_det.limitsStr()); + SensorBottomLogical.setLimitSet(theDetector, x_det.limitsStr()); + } + SensorTopLogical.setVisAttributes(theDetector.visAttributes(sensVis)); + SensorBottomLogical.setVisAttributes(theDetector.visAttributes(sensVis)); //create dead sensor logical volume Box SensorDeadSolid(sensor_thickness / 2.0, sensor_dead_width / 2.0, sensor_active_len / 2.0); @@ -294,11 +300,11 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h double ypos_active = (support_width/2.0) - (sensor_active_width/2.0); double ypos_dead = (-support_width/2.0) + sensor_deadwire_width + (sensor_dead_width/2.0); double zpos = -sensor_total_z/2.0 + sensor_active_len/2.0 + isensor*(sensor_active_len + dead_gap); - pv = SensorTopEnvelopeLogical.placeVolume(SensorLogical, Position(xpos,ypos_active,zpos)); + pv = SensorTopEnvelopeLogical.placeVolume(SensorTopLogical, Position(xpos,ypos_active,zpos)); //pv.addPhysVolID("topsensor", isensor ) ; pv.addPhysVolID("layer", layer_id*2+1).addPhysVolID("active", 0).addPhysVolID("sensor", isensor) ; TopSensor_pv.push_back(pv); - pv = SensorBottomEnvelopeLogical.placeVolume(SensorLogical, Position(xpos,ypos_active,zpos)); + pv = SensorBottomEnvelopeLogical.placeVolume(SensorBottomLogical, Position(xpos,ypos_active,zpos)); //pv.addPhysVolID("bottomsensor", isensor ) ; pv.addPhysVolID("layer", layer_id*2 ).addPhysVolID("active", 0).addPhysVolID("sensor", isensor) ; BottomSensor_pv.push_back(pv); @@ -309,7 +315,8 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h //place the sensor envelope inside the ladder envelope pv = LadderLogical.placeVolume(SensorTopEnvelopeLogical, Position(support_height/2.0 + flex_thickness + sensor_thickness/2.0, 0., 0.));//top-side sensors - Transform3D tran_sen(RotationZYX(0., dd4hep::twopi/2.0, 0.), Position(-(support_height/2.0 + flex_thickness + sensor_thickness/2.0), 0., 0.)); + //Transform3D tran_sen(RotationZYX(0., dd4hep::twopi/2.0, 0.), Position(-(support_height/2.0 + flex_thickness + sensor_thickness/2.0), 0., 0.)); + Transform3D tran_sen(RotationZYX(0., 0, 0.), Position(-(support_height/2.0 + flex_thickness + sensor_thickness/2.0), 0., 0.)); pv = LadderLogical.placeVolume(SensorBottomEnvelopeLogical,tran_sen);//bottom-side sensors //create the ladder support envelope @@ -321,11 +328,7 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h Box LadderSupportSolid(support_thickness / 2.0 , support_width / 2.0 , support_length / 2.0); Volume LadderSupportLogical(name + _toString( layer_id,"_SupLogical_%02d"), LadderSupportSolid, support_mat); LadderSupportLogical.setVisAttributes(theDetector.visAttributes(supportVis)); - - //vxd.setVisAttributes(theDetector, sensVis, SensorLogical); - // vxd.setVisAttributes(theDetector, sensEnvVis, SensorBottomEnvelopeLogical); - // vxd.setVisAttributes(theDetector, ladderVis, LadderLogical); - + pv = LadderSupportEnvelopeLogical.placeVolume(LadderSupportLogical); pv = LadderLogical.placeVolume(LadderSupportEnvelopeLogical); @@ -335,21 +338,28 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h DetElement ladderDE(layerDE, ladder_enum.str(), x_det.id()); std::cout << "start building " << ladder_enum.str() << ":" << endl; + //Transform3D tr (RotationZYX(ladder_dphi*i,0.,0.),Position(ladder_radius*cos(ladder_phi0+ladder_dphi*i), ladder_radius*sin(ladder_phi0+ladder_dphi*i), 0.)); + //pv = layer_assembly.placeVolume(LadderLogical,tr); + //pv.addPhysVolID("module", i); + //ladderDE.setPlacement(pv); + //std::cout << ladder_enum.str() << " done." << endl; + //if(i==0) std::cout << "xy=" << ladder_radius*cos(ladder_phi0) << " " << ladder_radius*sin(ladder_phi0) << std::endl; + //====== create the meassurement surface =================== dd4hep::rec::Vector3D o(0,0,0); - dd4hep::rec::Vector3D u( 0., 0., 1.); - dd4hep::rec::Vector3D v( 0., 1., 0.); - dd4hep::rec::Vector3D n( 1., 0., 0.); + dd4hep::rec::Vector3D u(0., 1., 0.); + dd4hep::rec::Vector3D v(0., 0., 1.); + dd4hep::rec::Vector3D n(1., 0., 0.); double inner_thick_top = sensor_thickness/2.0; double outer_thick_top = support_height/2.0 + flex_thickness + sensor_thickness/2.0; double inner_thick_bottom = support_height/2.0 + flex_thickness + sensor_thickness/2.0; double outer_thick_bottom = sensor_thickness/2.0; - dd4hep::rec::VolPlane surfTop( SensorLogical , - dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Sensitive), - inner_thick_top, outer_thick_top , u,v,n,o ) ; - dd4hep::rec::VolPlane surfBottom( SensorLogical , - dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Sensitive), - inner_thick_bottom, outer_thick_bottom, u,v,n,o ) ; + dd4hep::rec::VolPlane surfTop(SensorTopLogical, + dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Sensitive), + inner_thick_top, outer_thick_top, u, v, n, o); + dd4hep::rec::VolPlane surfBottom(SensorBottomLogical, + dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Sensitive), + inner_thick_bottom, outer_thick_bottom, u, v, n, o); for(int isensor=0; isensor < n_sensors_per_side; ++isensor){ diff --git a/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v02_geo.cpp b/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v02_geo.cpp new file mode 100644 index 0000000000000000000000000000000000000000..12352383e4a6909124f5180b53161b41bcf8581a --- /dev/null +++ b/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v02_geo.cpp @@ -0,0 +1,423 @@ +//==================================================================== +// cepcvxdgeo - CEPC vertex detector models in DD4hep +//-------------------------------------------------------------------- +// Hao Zeng, IHEP +// email: zenghao@ihep.ac.cn +// $Id$ +//==================================================================== +#include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/DD4hepUnits.h" +#include "DD4hep/DetType.h" +#include "DDRec/Surface.h" +#include "DDRec/DetectorData.h" +#include "XML/Utilities.h" +#include <cmath> + +using namespace std; + +using dd4hep::Box; +using dd4hep::Tube; +using dd4hep::DetElement; +using dd4hep::Material; +using dd4hep::Position; +using dd4hep::RotationY; +using dd4hep::RotationZYX; +using dd4hep::Transform3D; +using dd4hep::Rotation3D; +using dd4hep::Volume; +using dd4hep::_toString; +using dd4hep::rec::volSurfaceList; +using dd4hep::rec::ZPlanarData; +using dd4hep::mm; + +/** helper struct */ +struct VXD_Layer { + int n_ladders; + int n_sensors_per_ladder; + double sensor_length; + double half_z; + double sensitive_inner_radius ; + double support_inner_radius ; + double ladder_width ; + double ladder_dphi ; +}; + +//std::vector<VXD_Layer> _VXD_Layers; + +// /** helper struct */ +// struct extended_reconstruction_parameters { +// double sensor_length_mm; +// double strip_width_mm; +// double strip_length_mm; +// double strip_pitch_mm; +// double strip_angle_deg; +// }; + +//extended_reconstruction_parameters _e_r_p; + + +/** Construction of the VXD detector, ported from Mokka driver SIT_Simple_Pixel.cc + * + * Mokka History: + * Feb 7th 2011, Steve Aplin - original version + * F.Gaede, DESY, Jan 2014 - dd4hep SIT pixel + + * @author Hao Zeng, IHEP, July 2021 + */ +static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4hep::SensitiveDetector sens) { + + xml_det_t x_det = e; + Material air = theDetector.air(); + int det_id = x_det.id(); + string name = x_det.nameStr(); + DetElement vxd(name, det_id); + + Volume envelope = dd4hep::xml::createPlacedEnvelope(theDetector, e, vxd); + dd4hep::xml::setDetectorTypeFlag(e, vxd) ; + if(theDetector.buildType()==dd4hep::BUILD_ENVELOPE) return vxd; + envelope.setVisAttributes(theDetector.visAttributes("SeeThrough")); + + if (x_det.hasChild(_U(sensitive))) { + xml_dim_t sd_typ = x_det.child(_U(sensitive)); + sens.setType(sd_typ.typeStr()); + } + else { + sens.setType("tracker"); + } + std::cout << " ** building SiTrackerStaggeredLadder_v02 ... " << sens.type() << std::endl ; + + dd4hep::rec::ZPlanarData* zPlanarData = new dd4hep::rec::ZPlanarData; + + // fetch the global parameters + + //fetch the display parameters + xml_comp_t x_display(x_det.child(_Unicode(display))); + std::string ladderVis = x_display.attr<string>(_Unicode(ladder)); + std::string supportVis = x_display.attr<string>(_Unicode(support)); + std::string flexVis = x_display.attr<string>(_Unicode(flex)); + std::string sensEnvVis = x_display.attr<string>(_Unicode(sens_env)); + std::string sensVis = x_display.attr<string>(_Unicode(sens)); + std::string deadsensVis = x_display.attr<string>(_Unicode(deadsensor)); + std::string deadwireVis = x_display.attr<string>(_Unicode(deadwire)); + + //fetch the shell parameters + if (x_det.hasChild(_Unicode(shell))) { + xml_comp_t x_shell(x_det.child(_Unicode(shell))); + double rmin_shell = x_shell.rmin(); + double rmax_shell = x_shell.rmax(); + double zhalf_shell = x_shell.zhalf(); + Tube shellSolid(rmin_shell, rmax_shell, zhalf_shell); + Volume shellLogical(name + "_ShellLogical", shellSolid, theDetector.material(x_shell.materialStr())); + shellLogical.setVisAttributes(theDetector.visAttributes(x_shell.visStr())); + envelope.placeVolume(shellLogical); + + zPlanarData->zHalfShell = zhalf_shell; + zPlanarData->rInnerShell = rmin_shell; + zPlanarData->rOuterShell = rmax_shell; + } + + for(xml_coll_t layer_i(x_det,_U(layer)); layer_i; ++layer_i){ + xml_comp_t x_layer(layer_i); + + dd4hep::PlacedVolume pv; + int layer_id = x_layer.attr<int>(_Unicode(layer_id)); + + std::cout << "layer_id: " << layer_id << endl; + + double sensitive_radius = x_layer.attr<double>(_Unicode(ladder_radius)); + int n_sensors_per_ladder = x_layer.attr<int>(_Unicode(n_sensors_per_side)); + int n_ladders = x_layer.attr<int>(_Unicode(n_ladders)) ; + double ladder_offset = x_layer.attr<double>(_Unicode(ladder_offset)); + double ladder_radius = sqrt(ladder_offset*ladder_offset + sensitive_radius*sensitive_radius); + double ladder_phi0 = -atan(ladder_offset/sensitive_radius); + std::cout << "ladder_radius: " << ladder_radius/mm <<" mm" << endl; + + std::cout << "sensitive_radius: " << sensitive_radius/mm << " mm" << endl; + std::cout << "n_sensors_per_ladder: " << n_sensors_per_ladder << endl; + + std::string layerName = dd4hep::_toString( layer_id , "layer_%d" ); + dd4hep::Assembly layer_assembly( layerName ) ; + pv = envelope.placeVolume( layer_assembly ) ; + dd4hep::DetElement layerDE( vxd , layerName , x_det.id() ); + layerDE.setPlacement( pv ) ; + + const double ladder_dphi = ( dd4hep::twopi / n_ladders ) ; + std::cout << "ladder_dphi: " << ladder_dphi << endl; + + //fetch the ladder parameters + xml_comp_t x_ladder(x_layer.child(_Unicode(ladder))); + // db = XMLHandlerDB(x_ladder); + + //fetch the ladder support parameters + xml_comp_t x_ladder_support(x_ladder.child(_Unicode(ladderSupport))); + double support_length = x_ladder_support.attr<double>(_Unicode(length)); + double support_thickness = x_ladder_support.attr<double>(_Unicode(thickness)); + double support_height = x_ladder_support.attr<double>(_Unicode(height)); + double support_width = x_ladder_support.attr<double>(_Unicode(width)); + Material support_mat; + if(x_ladder_support.hasAttr(_Unicode(mat))) + { + support_mat = theDetector.material(x_ladder_support.attr<string>(_Unicode(mat))); + } + else + { + support_mat = theDetector.material(x_ladder_support.materialStr()); + } + std::cout << "support_length: " << support_length/mm << " mm" << endl; + std::cout << "support_thickness: " << support_thickness/mm << " mm" << endl; + std::cout << "support_width: " << support_width/mm << " mm" << endl; + + //fetch the flex parameters + double flex_thickness(0); + double flex_width(0); + double flex_length(0); + xml_comp_t x_flex(x_ladder.child(_Unicode(flex))); + for(xml_coll_t flex_i(x_flex,_U(slice)); flex_i; ++flex_i){ + xml_comp_t x_flex_slice(flex_i); + double x_flex_slice_thickness = x_flex_slice.attr<double>(_Unicode(thickness)); + double x_flex_slice_width = x_flex_slice.attr<double>(_Unicode(width)); + double x_flex_slice_length = x_flex_slice.attr<double>(_Unicode(length)); + flex_thickness += x_flex_slice_thickness; + if (x_flex_slice_width > flex_width) flex_width = x_flex_slice_width; + if (x_flex_slice_length > flex_length) flex_length = x_flex_slice_length; + //std::cout << "x_flex_slice_thickness: " << x_flex_slice_thickness/mm << " mm" << endl; + } + std::cout << "flex_thickness: " << flex_thickness/mm << " mm" << endl; + std::cout << "flex_width: " << flex_width/mm << " mm" << endl; + std::cout << "flex_length: " << flex_length/mm << " mm" << endl; + + //fetch the sensor parameters + xml_comp_t x_sensor(x_ladder.child(_Unicode(sensor))); + int n_sensors_per_side = x_sensor.attr<int>(_Unicode(n_sensors)); + double dead_gap = x_sensor.attr<double>(_Unicode(gap)); + double sensor_thickness = x_sensor.attr<double>(_Unicode(thickness)); + double sensor_active_len = x_sensor.attr<double>(_Unicode(active_length)); + double sensor_active_width = x_sensor.attr<double>(_Unicode(active_width)); + double sensor_dead_width = x_sensor.attr<double>(_Unicode(dead_width)); + double sensor_deadwire_length = x_sensor.attr<double>(_Unicode(deadwire_length)); + double sensor_deadwire_width = x_sensor.attr<double>(_Unicode(deadwire_width)); + double sensor_deadwire_thickness = x_sensor.attr<double>(_Unicode(deadwire_thickness)); + Material sensor_mat = theDetector.material(x_sensor.attr<string>(_Unicode(sensor_mat))); + Material sensor_deadwire_mat = theDetector.material(x_sensor.attr<string>(_Unicode(deadwire_mat))); + + std::cout << "n_sensors_per_side: " << n_sensors_per_side << endl; + std::cout << "dead_gap: " << dead_gap/mm << " mm" << endl; + std::cout << "sensor_thickness: " << sensor_thickness/mm << " mm" << endl; + std::cout << "sensor_active_len: " << sensor_active_len/mm << " mm" << endl; + std::cout << "sensor_active_width: " << sensor_active_width/mm << " mm" << endl; + std::cout << "sensor_dead_width: " << sensor_dead_width/mm << " mm" << endl; + + //create ladder logical volume + Box LadderSolid((support_height+2*sensor_thickness+2*flex_thickness)/2.0, + support_width / 2.0, support_length / 2.0); + Volume LadderLogical(name + dd4hep::_toString( layer_id, "_LadderLogical_%02d"), + LadderSolid, air); + // create flex envelope logical volume + Box FlexEnvelopeSolid(flex_thickness / 2.0, flex_width / 2.0, flex_length / 2.0); + Volume FlexEnvelopeLogical(name + dd4hep::_toString( layer_id, "_FlexEnvelopeLogical_%02d"), FlexEnvelopeSolid, air); + FlexEnvelopeLogical.setVisAttributes(theDetector.visAttributes("SeeThrough")); + + //create the flex layers inside the flex envelope + double flex_start_height(-flex_thickness/2.); + int index = 0; + for(xml_coll_t flex_i(x_flex,_U(slice)); flex_i; ++flex_i){ + xml_comp_t x_flex_slice(flex_i); + double x_flex_slice_thickness = x_flex_slice.attr<double>(_Unicode(thickness)); + double x_flex_slice_width = x_flex_slice.attr<double>(_Unicode(width)); + double x_flex_slice_length = x_flex_slice.attr<double>(_Unicode(length)); + Material x_flex_slice_mat; + if(x_flex_slice.hasAttr(_Unicode(mat))) { + x_flex_slice_mat = theDetector.material(x_flex_slice.attr<string>(_Unicode(mat))); + } + else { + x_flex_slice_mat = theDetector.material(x_flex_slice.materialStr()); + } + // Material x_flex_slice_mat = theDetector.material(x_flex_slice.attr<string>(_Unicode(mat))); + Box FlexLayerSolid(x_flex_slice_thickness/2.0, x_flex_slice_width/2.0, x_flex_slice_length/2.0); + Volume FlexLayerLogical(name + dd4hep::_toString(layer_id, "_FlexLayerLogical_%02d") + dd4hep::_toString(index, "index_%02d"), FlexLayerSolid, x_flex_slice_mat); + FlexLayerLogical.setVisAttributes(theDetector.visAttributes(flexVis)); + double flex_slice_height = flex_start_height + x_flex_slice_thickness/2.; + pv = FlexEnvelopeLogical.placeVolume(FlexLayerLogical, Position(flex_slice_height, 0., 0.)); + //std::cout << "flex thickness = " << x_flex_slice_thickness << std::endl; + //std::cout << "flex width = " << x_flex_slice_width << std::endl; + //std::cout << "flex length = " << x_flex_slice_length << std::endl; + //std::cout << "flex material: " << x_flex_slice_mat << std::endl; + flex_start_height += x_flex_slice_thickness; + index++; + } + + //place the flex envelope inside the ladder envelope + pv = LadderLogical.placeVolume(FlexEnvelopeLogical, Position((support_height + flex_thickness) / 2.0, 0., 0.)); //top side + //define the transformation3D(only need a combination of translation and rotation) + Transform3D tran_mirro(RotationZYX(0., dd4hep::twopi/2.0, 0.), Position(-(support_height + flex_thickness) / 2.0, 0., 0.)); + pv = LadderLogical.placeVolume(FlexEnvelopeLogical, tran_mirro); //bottom side + + //create sensor envelope logical volume + Box SensorTopEnvelopeSolid(sensor_thickness / 2.0, support_width / 2.0, support_length / 2.0); + Volume SensorTopEnvelopeLogical(name + dd4hep::_toString( layer_id, "_SensorEnvelopeLogical_%02d"), SensorTopEnvelopeSolid, air); + Box SensorBottomEnvelopeSolid(sensor_thickness / 2.0, support_width / 2.0, support_length / 2.0); + Volume SensorBottomEnvelopeLogical(name + dd4hep::_toString( layer_id, "_SensorEnvelopeLogical_%02d"), SensorBottomEnvelopeSolid, air); + SensorTopEnvelopeLogical.setVisAttributes(theDetector.visAttributes(sensEnvVis)); + + //create sensor logical volume + Box SensorSolid(sensor_thickness / 2.0, sensor_active_width / 2.0, sensor_active_len / 2.0); + Volume SensorTopLogical(name + dd4hep::_toString(layer_id*2+1, "_SensorLogical_%02d"), SensorSolid, sensor_mat); + Volume SensorBottomLogical(name + dd4hep::_toString(layer_id*2, "_SensorLogical_%02d"), SensorSolid, sensor_mat); + SensorTopLogical.setSensitiveDetector(sens); + SensorBottomLogical.setSensitiveDetector(sens); + if (x_det.hasAttr(_U(limits))) { + SensorTopLogical.setLimitSet(theDetector, x_det.limitsStr()); + SensorBottomLogical.setLimitSet(theDetector, x_det.limitsStr()); + } + SensorTopLogical.setVisAttributes(theDetector.visAttributes(sensVis)); + SensorBottomLogical.setVisAttributes(theDetector.visAttributes(sensVis)); + + //create dead sensor logical volume + Box SensorDeadSolid(sensor_thickness / 2.0, sensor_dead_width / 2.0, sensor_active_len / 2.0); + Volume SensorDeadLogical(name + dd4hep::_toString( layer_id, "_SensorDeadLogical_%02d"), SensorDeadSolid, sensor_mat); + SensorDeadLogical.setVisAttributes(theDetector.visAttributes(deadsensVis)); + + //create dead wire logical volume + Box SensorDeadWireSolid(sensor_deadwire_thickness / 2.0, sensor_deadwire_width / 2.0, sensor_deadwire_length / 2.0); + Volume SensorDeadWireLogical(name + dd4hep::_toString( layer_id, "_SensorDeadWireLogical_%02d"), SensorDeadWireSolid, sensor_deadwire_mat); + SensorDeadWireLogical.setVisAttributes(theDetector.visAttributes(deadwireVis)); + + //place the dead wire in the sensor envelope + // pv = SensorTopEnvelopeLogical.placeVolume(SensorDeadWireLogical, Position(0.0, (sensor_active_width-support_width/2.0) + sensor_dead_width/2.0 + sensor_deadwire_width/2.0, 0.0)); + // pv = SensorBottomEnvelopeLogical.placeVolume(SensorDeadWireLogical, Position(0.0, (sensor_active_width-support_width/2.0) + sensor_dead_width/2.0 + sensor_deadwire_width/2.0, 0.0)); + pv = SensorTopEnvelopeLogical.placeVolume(SensorDeadWireLogical, Position(0.0, (-support_width/2.0) + (sensor_deadwire_width/2.0), 0.0)); + pv = SensorBottomEnvelopeLogical.placeVolume(SensorDeadWireLogical, Position(0.0, (-support_width/2.0) + (sensor_deadwire_width/2.0), 0.0)); + + // place the active sensor and dead sensor inside the sensor envelope + std::vector<dd4hep::PlacedVolume> TopSensor_pv; + std::vector<dd4hep::PlacedVolume> BottomSensor_pv; + for(int isensor=0; isensor < n_sensors_per_side; ++isensor){ + double sensor_total_z = n_sensors_per_side*sensor_active_len + dead_gap*(n_sensors_per_side-1); + double xpos = 0.0; + double ypos_active = (support_width/2.0) - (sensor_active_width/2.0); + double ypos_dead = (-support_width/2.0) + sensor_deadwire_width + (sensor_dead_width/2.0); + double zpos = -sensor_total_z/2.0 + sensor_active_len/2.0 + isensor*(sensor_active_len + dead_gap); + pv = SensorTopEnvelopeLogical.placeVolume(SensorTopLogical, Position(xpos,ypos_active,zpos)); + //pv.addPhysVolID("topsensor", isensor ) ; + pv.addPhysVolID("layer", layer_id*2+1).addPhysVolID("active", 0).addPhysVolID("sensor", isensor) ; + TopSensor_pv.push_back(pv); + pv = SensorBottomEnvelopeLogical.placeVolume(SensorBottomLogical, Position(xpos,ypos_active,zpos)); + //pv.addPhysVolID("bottomsensor", isensor ) ; + pv.addPhysVolID("layer", layer_id*2 ).addPhysVolID("active", 0).addPhysVolID("sensor", isensor) ; + BottomSensor_pv.push_back(pv); + pv = SensorTopEnvelopeLogical.placeVolume(SensorDeadLogical, Position(xpos,ypos_dead,zpos)); + pv = SensorBottomEnvelopeLogical.placeVolume(SensorDeadLogical, Position(xpos,ypos_dead,zpos)); + + } + //place the sensor envelope inside the ladder envelope + pv = LadderLogical.placeVolume(SensorTopEnvelopeLogical, + Position(support_height/2.0 + flex_thickness + sensor_thickness/2.0, 0., 0.));//top-side sensors + Transform3D tran_sen(RotationZYX(0., 0., 0.), Position(-(support_height/2.0 + flex_thickness + sensor_thickness/2.0), 0., 0.)); + pv = LadderLogical.placeVolume(SensorBottomEnvelopeLogical,tran_sen);//bottom-side sensors + + //create the ladder support + Box LadderSupportSolid(support_height/2.0, support_width/2.0, support_length/2.0); + Volume LadderSupportLogical(name + _toString( layer_id,"_SupLogical_%02d"), LadderSupportSolid, support_mat); + LadderSupportLogical.setVisAttributes(theDetector.visAttributes(supportVis)); + + //create ladder support cavity + Box LadderSupportCavitySolid(support_height/2.0-support_thickness/2.0, support_width/2.0-support_thickness/2.0, support_length/2.0); + Volume LadderSupportCavityLogical(name + _toString( layer_id,"_SupCavityLogical_%02d"), LadderSupportCavitySolid, air); + LadderSupportCavityLogical.setVisAttributes(theDetector.visAttributes("SeeThrough")); + + pv = LadderSupportLogical.placeVolume(LadderSupportCavityLogical); + pv = LadderLogical.placeVolume(LadderSupportLogical); + + for(int i = 0; i < n_ladders; i++){ + std::stringstream ladder_enum; + ladder_enum << "vxt_ladder_" << layer_id << "_" << i; + DetElement ladderDE(layerDE, ladder_enum.str(), x_det.id()); + //std::cout << "start building " << ladder_enum.str() << ":" << endl; + + //====== create the meassurement surface =================== + dd4hep::rec::Vector3D o(0,0,0); + dd4hep::rec::Vector3D u(0., 1., 0.); + dd4hep::rec::Vector3D v(0., 0., 1.); + dd4hep::rec::Vector3D n(1., 0., 0.); + double inner_thick_top = sensor_thickness/2.0; + double outer_thick_top = support_height/2.0 + flex_thickness + sensor_thickness/2.0; + double inner_thick_bottom = support_height/2.0 + flex_thickness + sensor_thickness/2.0; + double outer_thick_bottom = sensor_thickness/2.0; + dd4hep::rec::VolPlane surfTop(SensorTopLogical, + dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Sensitive), + inner_thick_top, outer_thick_top , u, v, n, o); + dd4hep::rec::VolPlane surfBottom(SensorBottomLogical, + dd4hep::rec::SurfaceType(dd4hep::rec::SurfaceType::Sensitive), + inner_thick_bottom, outer_thick_bottom, u, v, n, o); + + + for(int isensor=0; isensor < n_sensors_per_side; ++isensor){ + std::stringstream topsensor_str; + std::stringstream bottomsensor_str; + topsensor_str << ladder_enum.str() << "_top_" << isensor; + // std::cout << "\tstart building " << topsensor_str.str() << ":" << endl; + bottomsensor_str << ladder_enum.str() << "_bottom_" << isensor; + // std::cout << "\tstart building " << bottomsensor_str.str() << ":" << endl; + DetElement topsensorDE(ladderDE, topsensor_str.str(), x_det.id()); + DetElement bottomsensorDE(ladderDE, bottomsensor_str.str(), x_det.id()); + topsensorDE.setPlacement(TopSensor_pv[isensor]); + volSurfaceList(topsensorDE)->push_back(surfTop); + // std::cout << "\t" << topsensor_str.str() << " done." << endl; + bottomsensorDE.setPlacement(BottomSensor_pv[isensor]); + // std::cout << "\t" << bottomsensor_str.str() << " done." << endl; + volSurfaceList(bottomsensorDE)->push_back(surfBottom); + } + Transform3D tr (RotationZYX(ladder_dphi*i,0.,0.),Position(ladder_radius*cos(ladder_phi0+ladder_dphi*i), ladder_radius*sin(ladder_phi0+ladder_dphi*i), 0.)); + pv = layer_assembly.placeVolume(LadderLogical,tr); + pv.addPhysVolID("module", i ) ; + ladderDE.setPlacement(pv); + //std::cout << ladder_enum.str() << " done." << endl; + if(i==0) std::cout << "xy=" << ladder_radius*cos(ladder_phi0) << " " << ladder_radius*sin(ladder_phi0) << std::endl; + } + + // package the reconstruction data + dd4hep::rec::ZPlanarData::LayerLayout topLayer; + dd4hep::rec::ZPlanarData::LayerLayout bottomLayer; + + topLayer.ladderNumber = n_ladders; + topLayer.phi0 = 0.; + topLayer.sensorsPerLadder = n_sensors_per_side; + topLayer.lengthSensor = sensor_active_len; + topLayer.distanceSupport = sensitive_radius; + topLayer.thicknessSupport = support_thickness / 2.0; + topLayer.offsetSupport = -ladder_offset; + topLayer.widthSupport = support_width; + topLayer.zHalfSupport = support_length / 2.0; + topLayer.distanceSensitive = sensitive_radius + support_height / 2.0 + flex_thickness; + topLayer.thicknessSensitive = sensor_thickness; + topLayer.offsetSensitive = -ladder_offset + (support_width/2.0 - sensor_active_width/2.0); + topLayer.widthSensitive = sensor_active_width; + topLayer.zHalfSensitive = (n_sensors_per_side*(sensor_active_len + dead_gap) - dead_gap) / 2.0; + + bottomLayer.ladderNumber = n_ladders; + bottomLayer.phi0 = 0.; + bottomLayer.sensorsPerLadder = n_sensors_per_side; + bottomLayer.lengthSensor = sensor_active_len; + bottomLayer.distanceSupport = sensitive_radius - support_height / 2.0 - flex_thickness; + bottomLayer.thicknessSupport = support_thickness / 2.0; + bottomLayer.offsetSupport = -ladder_offset; + bottomLayer.widthSupport = support_width; + bottomLayer.zHalfSupport = support_length / 2.0; + bottomLayer.distanceSensitive = sensitive_radius - support_height / 2.0 - sensor_thickness - flex_thickness; + bottomLayer.thicknessSensitive = sensor_thickness; + bottomLayer.offsetSensitive = -ladder_offset + (support_width/2.0 - sensor_active_width/2.0); + bottomLayer.widthSensitive = sensor_active_width; + bottomLayer.zHalfSensitive = (n_sensors_per_side*(sensor_active_len + dead_gap) - dead_gap) / 2.0; + + zPlanarData->layers.push_back(bottomLayer); + zPlanarData->layers.push_back(topLayer); + } + std::cout << (*zPlanarData) << endl; + vxd.addExtension< ZPlanarData >(zPlanarData); + if ( x_det.hasAttr(_U(combineHits)) ) { + vxd.setCombineHits(x_det.attr<bool>(_U(combineHits)),sens); + } + std::cout << "vxd done." << endl; + return vxd; +} +DECLARE_DETELEMENT(SiTrackerStaggeredLadder_v02,create_element) diff --git a/Detector/DetCRD/src/Tracker/SiTrackerStitching_v01_geo.cpp b/Detector/DetCRD/src/Tracker/SiTrackerStitching_v01_geo.cpp new file mode 100644 index 0000000000000000000000000000000000000000..1ac541b8517bfb9265f42f2192d1bfa51e89ad4c --- /dev/null +++ b/Detector/DetCRD/src/Tracker/SiTrackerStitching_v01_geo.cpp @@ -0,0 +1,276 @@ +//==================================================================== +// cepcvxdgeo - CEPC vertex detector models in DD4hep +//-------------------------------------------------------------------- +// Hao Zeng, IHEP +// email: zenghao@ihep.ac.cn +// $Id$ +//==================================================================== +#include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/DD4hepUnits.h" +#include "DD4hep/DetType.h" +#include "DDRec/Surface.h" +#include "DDRec/DetectorData.h" +#include "XML/Utilities.h" +#include <cmath> + +#include "Identifier/CEPCDetectorData.h" + +using namespace std; + +using dd4hep::Box; +using dd4hep::Tube; +using dd4hep::DetElement; +using dd4hep::Material; +using dd4hep::Position; +using dd4hep::RotationZ; +using dd4hep::RotationZYX; +using dd4hep::Transform3D; +using dd4hep::Rotation3D; +using dd4hep::Volume; +using dd4hep::_toString; +using dd4hep::rec::volSurfaceList; +using dd4hep::rec::CylindricalData; +using dd4hep::mm; + +/** helper struct */ +struct VXD_Module { + // int n_ladders; + //int n_sensors_per_ladder; + //double sensor_length; + double half_z; + double sensitive_inner_radius; + double support_inner_radius; + double support_thickness; + double width; + double ladder_dphi; +}; + +/** Construction of the VXD detector, stitching sensor + * @author + */ +static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4hep::SensitiveDetector sens) { + + xml_det_t x_det = e; + Material air = theDetector.air(); + int det_id = x_det.id(); + string name = x_det.nameStr(); + DetElement det(name, det_id); + + Volume envelope = dd4hep::xml::createPlacedEnvelope(theDetector, e, det); + dd4hep::xml::setDetectorTypeFlag(e, det); + if (theDetector.buildType()==dd4hep::BUILD_ENVELOPE) return det; + envelope.setVisAttributes(theDetector.visAttributes("SeeThrough")); + + if (x_det.hasChild(_U(sensitive))) { + xml_dim_t sd_typ = x_det.child(_U(sensitive)); + sens.setType(sd_typ.typeStr()); + } + else { + sens.setType("tracker"); + } + std::cout << " ** building SiTrackerStitchingLadder_v01 ... " << sens.type() << std::endl ; + + dd4hep::rec::CylindricalData* cylindricalData = new dd4hep::rec::CylindricalData; + //fetch the shell parameters + if (x_det.hasChild(_Unicode(shell))) { + xml_comp_t x_shell(x_det.child(_Unicode(shell))); + double rmin_shell = x_shell.rmin(); + double rmax_shell = x_shell.rmax(); + double zhalf_shell = x_shell.zhalf(); + Tube shellSolid(rmin_shell, rmax_shell, zhalf_shell); + Volume shellLogical(name + "_ShellLogical", shellSolid, theDetector.material(x_shell.materialStr())); + shellLogical.setVisAttributes(theDetector.visAttributes(x_shell.visStr())); + dd4hep::PlacedVolume pv = envelope.placeVolume(shellLogical); + + DetElement shellDE(det, name + "_Shell", x_det.id()); + shellDE.setPlacement(pv); + + cylindricalData->zHalfShell = zhalf_shell; + cylindricalData->rInnerShell = rmin_shell; + cylindricalData->rOuterShell = rmax_shell; + } + + std::map<std::string, Volume> module_volumes; + std::map<std::string, VXD_Module> module_data; + + xml_comp_t x_modules(x_det.child(_U(modules))); + for (xml_coll_t module_i(x_modules, _U(module)); module_i; ++module_i) { + xml_comp_t x_module(module_i); + std::string type = x_module.typeStr(); + double radius = x_module.radius(); + double xdead = x_module.attr<double>(_Unicode(xdead)); + double ydead = x_module.attr<double>(_Unicode(ydead)); + int nx = x_module.attr<int>(_Unicode(nx)); + int ny = x_module.attr<int>(_Unicode(ny)); + Material mat = theDetector.material(x_module.materialStr()); + + xml_comp_t x_sensor(x_module.child(_U(sensor))); + double sensor_thickness = x_sensor.thickness(); + double sensor_width = x_sensor.width(); + double sensor_length = x_sensor.length(); + double sensor_dphi = sensor_width/radius; + Material sensor_mat = theDetector.material(x_sensor.materialStr()); + + Tube sensor_solid(radius, radius+sensor_thickness, sensor_length/2, -sensor_width/radius/2, sensor_width/radius/2); + Volume sensor_volume(name + type + "Sensor", sensor_solid, sensor_mat); + sensor_volume.setSensitiveDetector(sens); + if (x_det.hasAttr(_U(limits))) sensor_volume.setLimitSet(theDetector, x_det.limitsStr()); + sensor_volume.setVisAttributes(theDetector.visAttributes(x_sensor.visStr())); + + xml_comp_t x_flex(x_module.child(_Unicode(flex))); + std::vector<std::pair<double, Material> > flexs; + double flex_thickness = 0; + for (xml_coll_t slice_i(x_flex, _U(slice)); slice_i; ++slice_i) { + xml_comp_t x_slice(slice_i); + double thickness = x_slice.thickness(); + Material mat = theDetector.material(x_slice.materialStr()); + flexs.push_back(std::make_pair(thickness, mat)); + flex_thickness += thickness; + std::cout << "flex thickness: " << thickness << std::endl; + } + double flex_width = sensor_width*ny + ydead*(ny-1); + double flex_length = sensor_length*nx + xdead*(nx-1); + double flex_radius = radius+sensor_thickness; + double flex_dphi = flex_width/radius; + Tube flex_solid(flex_radius, flex_radius+flex_thickness, flex_length/2, -flex_dphi/2, flex_dphi/2); + Volume flex_volume(name + type + "Flex", flex_solid, air); + flex_volume.setVisAttributes(theDetector.visAttributes(x_flex.visStr())); + double start_radius = flex_radius; + for (unsigned islice=0; islice<flexs.size(); islice++) { + //std::cout << "flex start radius " << start_radius << endl; + Tube slice_solid(start_radius, start_radius+flexs[islice].first, flex_length/2, -flex_dphi/2, flex_dphi/2); + Volume slice_volume(name + type + dd4hep::_toString(int(islice), "Flex_%d"), slice_solid, flexs[islice].second); + flex_volume.placeVolume(slice_volume); + start_radius += flexs[islice].first; + } + + double service_radius = radius + sensor_thickness; + + xml_comp_t x_electronics(x_module.child(_Unicode(electronics))); + double electronics_thickness = x_electronics.thickness(); + double electronics_width = x_electronics.width(); + double electronics_dphi = electronics_width/radius; + Material electronics_mat = theDetector.material(x_electronics.materialStr()); + Tube electronics_solid(service_radius, service_radius+electronics_thickness, flex_length/2, -electronics_dphi, 0); + Volume electronics_volume(name + type + "Electronics", electronics_solid, electronics_mat); + electronics_volume.setVisAttributes(theDetector.visAttributes(x_electronics.visStr())); + + xml_comp_t x_readout(x_module.child(_U(readout))); + double readout_thickness = x_readout.thickness(); + double readout_width = x_readout.width(); + double readout_dphi = readout_width/radius; + Material readout_mat = theDetector.material(x_readout.materialStr()); + Tube readout_solid(service_radius, service_radius+readout_thickness, flex_length/2, 0, readout_dphi); + Volume readout_volume(name + type + "Electronics", readout_solid, readout_mat); + readout_volume.setVisAttributes(theDetector.visAttributes(x_readout.visStr())); + + xml_comp_t x_driver(x_module.child(_Unicode(driver))); + double driver_thickness = x_driver.thickness(); + // width -> length, length->width + double driver_length = x_driver.width(); + double driver_width = flex_width + electronics_width + readout_width; + double driver_dphi = driver_width/radius; + Material driver_mat = theDetector.material(x_driver.materialStr()); + Tube driver_solid(service_radius, service_radius+driver_thickness, driver_length/2, -flex_dphi/2-electronics_dphi, flex_dphi/2+readout_dphi); + Volume driver_volume(name + type + "Driver", driver_solid, driver_mat); + driver_volume.setVisAttributes(theDetector.visAttributes(x_driver.visStr())); + + double module_length = flex_length + 2*driver_length; + double module_width = driver_width; + double module_dphi = module_width/radius; + double module_thickness = sensor_thickness + std::max(std::max(std::max(flex_thickness, driver_thickness), readout_thickness), electronics_thickness); + Tube module_solid(radius, radius+module_thickness, module_length/2, -flex_dphi/2-electronics_dphi, flex_dphi/2+readout_dphi); + Volume module_volume(name + type, module_solid, air); + module_volume.setVisAttributes(theDetector.visAttributes("SeeThrough")); + + Tube board_solid(radius, radius+sensor_thickness, module_length/2, -flex_dphi/2-electronics_dphi, flex_dphi/2+readout_dphi); + Volume board_volume(name + type + "Board", board_solid, mat); + board_volume.setVisAttributes(theDetector.visAttributes(x_module.visStr())); + module_volume.placeVolume(board_volume); + + for (int ix=0; ix<nx; ix++) { + double z = -flex_length/2+sensor_length/2+ix*(sensor_length+xdead); + for (int iy=0; iy<ny; iy++) { + double delta = -flex_dphi/2+sensor_dphi/2+iy*(sensor_dphi+ydead/radius); + Transform3D tran(RotationZ(delta), Position(0, 0, z)); + dd4hep::PlacedVolume pv = board_volume.placeVolume(sensor_volume, tran); + + int sensor_id = ix + nx*iy; + pv.addPhysVolID("sensor", sensor_id); + } + } + + module_volume.placeVolume(flex_volume); + module_volume.placeVolume(electronics_volume, RotationZYX(-flex_dphi/2,0,0)); + module_volume.placeVolume(readout_volume, RotationZYX(flex_dphi/2,0,0)); + module_volume.placeVolume(driver_volume, Position(0, 0, -flex_length/2-driver_length/2)); + module_volume.placeVolume(driver_volume, Position(0, 0, flex_length/2+driver_length/2)); + module_volumes[type] = module_volume; + + VXD_Module thisModule; + thisModule.half_z = flex_length/2; + thisModule.sensitive_inner_radius = radius; + thisModule.support_inner_radius = flex_radius; + thisModule.support_thickness = flex_thickness; + thisModule.width = flex_width; + module_data[type] = thisModule; + } + + for(xml_coll_t layer_i(x_det,_U(layer)); layer_i; ++layer_i){ + xml_comp_t x_layer(layer_i); + + int layer_id = x_layer.id(); + std::cout << "layer_id: " << layer_id << endl; + + dd4hep::rec::CylindricalData::LayerLayout thisLayer; + + double phi0 = x_layer.phi0(); + + //thisLayer.phi0 = phi0; + thisLayer.id = layer_id; + + int module_id = 0; + for (xml_coll_t module_i(x_layer, _U(module)); module_i; ++module_i) { + if (module_id>1) std::cerr << "Not support more than two modules, possible wrong in reconstruction" << std::endl; + xml_comp_t x_module(module_i); + std::string type = x_module.typeStr(); + std::cout << ">>Module " << type << endl; + double offset = x_module.offset(); + double phi = x_module.phi(); + double z = x_module.z(); + Transform3D tran(RotationZ(phi0+phi), Position(offset*cos(phi), offset*sin(phi), z)); + dd4hep::PlacedVolume pv = envelope.placeVolume(module_volumes[type], tran); + pv.addPhysVolID("layer", layer_id).addPhysVolID("module", module_id); + + DetElement moduleDE(det, name + dd4hep::_toString(layer_id, "_Layer%02d") + dd4hep::_toString(module_id, "_Stave%02d"), x_det.id()); + moduleDE.setPlacement(pv); + + VXD_Module thisModule = module_data[type]; + if (module_id==0) { + thisLayer.phi0 = phi0+phi; + thisLayer.zHalf = thisModule.half_z; + thisLayer.radiusSensitive = thisModule.sensitive_inner_radius; + thisLayer.thicknessSensitive = thisModule.support_inner_radius - thisModule.sensitive_inner_radius; + thisLayer.radiusSupport = thisModule.support_inner_radius; + thisLayer.thicknessSupport = thisModule.support_thickness; + } + else if (module_id==1) { + thisLayer.rgap = thisModule.sensitive_inner_radius - thisLayer.radiusSensitive; + thisLayer.dphi = phi0+phi-thisLayer.phi0; + } + + module_id++; + } + cylindricalData->layers.push_back(thisLayer); + } + + std::cout << (*cylindricalData) << std::endl; + det.addExtension<dd4hep::rec::CylindricalData>( cylindricalData ); + + if ( x_det.hasAttr(_U(combineHits)) ) { + det.setCombineHits(x_det.attr<bool>(_U(combineHits)),sens); + } + std::cout << name << " done." << endl; + return det; +} +DECLARE_DETELEMENT(SiTrackerStitching_v01,create_element) diff --git a/Detector/DetInterface/include/DetInterface/IGeomSvc.h b/Detector/DetInterface/include/DetInterface/IGeomSvc.h index 51c26406a307d3e66724d207a730bbeca823c216..5c63e5a64dd2e3683e14cf4e0b5c1c0b8db7f48a 100644 --- a/Detector/DetInterface/include/DetInterface/IGeomSvc.h +++ b/Detector/DetInterface/include/DetInterface/IGeomSvc.h @@ -12,6 +12,7 @@ #include "GaudiKernel/IService.h" #include "DDRec/DetectorData.h" +#include "DDRec/SurfaceManager.h" #include <map> namespace dd4hep { @@ -43,6 +44,7 @@ public: // short cut to retrieve the Decoder according to the Readout name virtual Decoder* getDecoder(const std::string& readout_name) = 0; + virtual const dd4hep::rec::SurfaceMap* getSurfaceMap(const std::string& det_name) = 0; virtual ~IGeomSvc() {} }; diff --git a/Detector/GeomSvc/src/GeomSvc.cpp b/Detector/GeomSvc/src/GeomSvc.cpp index 31443757ef16cb43bd2115d1de4c28929b7bb836..a35dc650e06695e258b0fa5bcf03b00c1d591a9f 100644 --- a/Detector/GeomSvc/src/GeomSvc.cpp +++ b/Detector/GeomSvc/src/GeomSvc.cpp @@ -35,9 +35,13 @@ GeomSvc::initialize() { StatusCode sc = Service::initialize(); m_dd4hep_geo = &(dd4hep::Detector::getInstance()); + // if DEBUG, set + //dd4hep::PrintLevel level = msgLevel(MSG::INFO) ? dd4hep::printLevel() : dd4hep::setPrintLevel(dd4hep::DEBUG); // if failed to load the compact, a runtime error will be thrown. m_dd4hep_geo->fromCompact(m_dd4hep_xmls.value()); - + // recover to old level, if not, too many DD4hep print + //dd4hep::setPrintLevel(level); + return sc; } @@ -45,6 +49,8 @@ StatusCode GeomSvc::finalize() { StatusCode sc; + // m_surface_manager has added as extension of Detector, so not delete? + //if (m_surface_manager != nullptr) delete m_surface_manager; dd4hep::Detector::destroyInstance(); return sc; @@ -95,3 +101,26 @@ GeomSvc::getDecoder(const std::string& readout_name) { return decoder; } + +const dd4hep::rec::SurfaceMap* +GeomSvc::getSurfaceMap(const std::string& det_name) { + if (m_surface_manager == nullptr) { + dd4hep::rec::SurfaceManager* surfaceMgr = nullptr; + // first check whether exist + try { + surfaceMgr = m_dd4hep_geo->extension<dd4hep::rec::SurfaceManager>(); + } + catch (std::runtime_error& e) { + info() << e.what() << " " << surfaceMgr << endmsg; + surfaceMgr = nullptr; + } + + if (surfaceMgr) { + m_surface_manager = surfaceMgr; + } + else { + m_dd4hep_geo->addExtension<dd4hep::rec::SurfaceManager>(m_surface_manager = new dd4hep::rec::SurfaceManager(*m_dd4hep_geo)); + } + } + return m_surface_manager->map(det_name); +} diff --git a/Detector/GeomSvc/src/GeomSvc.h b/Detector/GeomSvc/src/GeomSvc.h index 119e2ea5870db041ad48bf4bf0175af290af51cd..fb253bb189e5f00ff28b9cfd51fed1fec460d4d7 100644 --- a/Detector/GeomSvc/src/GeomSvc.h +++ b/Detector/GeomSvc/src/GeomSvc.h @@ -36,6 +36,7 @@ class GeomSvc: public extends<Service, IGeomSvc> { private: Decoder* getDecoder(const std::string& readout_name) override; + const dd4hep::rec::SurfaceMap* getSurfaceMap(const std::string& det_name) override; private: // DD4hep XML compact file path @@ -43,6 +44,7 @@ private: // dd4hep::Detector* m_dd4hep_geo; + dd4hep::rec::SurfaceManager* m_surface_manager = nullptr; }; #endif // GeomSvc_h diff --git a/Detector/Identifier/include/Identifier/CEPCConf.h b/Detector/Identifier/include/Identifier/CEPCConf.h index 95544ffcfecc175327192a43c41482eb62ed4829..a6e47b6b987d3fd81836735bbb3ed76f28d30709 100644 --- a/Detector/Identifier/include/Identifier/CEPCConf.h +++ b/Detector/Identifier/include/Identifier/CEPCConf.h @@ -37,6 +37,7 @@ namespace CEPCConf{ static const int ONE_DIMENSIONAL = 29; static const int COMPOSITE_SPACEPOINT = 30; static const int PLANAR = 3; // 3 is compatible with old tracking codes, 31 or 28 is better in future to modify uniformly + static const int CYLINDER = 28; }; struct TrkHitQualityBit{ diff --git a/Detector/Identifier/include/Identifier/CEPCDetectorData.h b/Detector/Identifier/include/Identifier/CEPCDetectorData.h new file mode 100644 index 0000000000000000000000000000000000000000..2cf98d3d2045807cea2f5edfc92acc5216cbf523 --- /dev/null +++ b/Detector/Identifier/include/Identifier/CEPCDetectorData.h @@ -0,0 +1,192 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// Author : +// +//========================================================================== +#ifndef CEPC_DETECTORDATA_H +#define CEPC_DETECTORDATA_H + +#include "DDRec/DetectorData.h" +#include <boost/io/ios_state.hpp> + +namespace dd4hep { + namespace rec { + /** Simple data structure with key parameters for + * reconstruction of a cylindrical detector + * + * @author + * @date Sep, 02 2024 + * @version $Id: $ + */ + struct CylindricalStruct { + /// The half length (z) of the support shell (w/o gap) - 0. if no shell exists. + double zHalfShell; + /// The length of the gap in mm (gap position at z=0). + double gapShell; + /// The inner radius of the support shell. + double rInnerShell; + /// The outer radius of the support shell. + double rOuterShell; + + /**Internal helper struct for defining the layer layout. Layers are defined + * with a sensitive part and a support part. + */ + struct LayerLayout { + // default c'tor with zero initialization + LayerLayout() : + id(0), + zHalf(0), + radiusSensitive(0), + thicknessSensitive(0), + radiusSupport(0), + thicknessSupport(0), + phi0(0), + rgap(0), + dphi(0), + offset(0) { + } + int id; + /// half length + double zHalf; + /// inner radius + double radiusSensitive; + /// thickness + double thicknessSensitive; + /// + double radiusSupport; + /// + double thicknessSupport; + /// + double phi0; + /// + double rgap; + /// + double dphi; + /// + double offset; + }; + std::vector<LayerLayout> layers; + }; + typedef StructExtension<CylindricalStruct> CylindricalData; + + static std::ostream& operator<<( std::ostream& io , const CylindricalData& d ) { + boost::io::ios_base_all_saver ifs(io); + + io << " -- CylindricalData: " << std::scientific << std::endl; + io << " zHalfShell : " << d.zHalfShell << std::endl; + io << " gapShell : " << d.gapShell << std::endl; + io << " rInnerShell : " << d.rInnerShell << std::endl; + io << " rOuterShell : " << d.rOuterShell << std::endl; + + std::vector<CylindricalData::LayerLayout> layers = d.layers; + + io << " Layers : " << std::endl + << " phi0 zHalf radiusSensitive thicknessSensitive radiusSupport thicknessSupport radialgap deltaphi" << std::endl; + + for(unsigned i=0,N=layers.size(); i<N; ++i){ + + CylindricalData::LayerLayout l = layers[i]; + + io << " " << l.phi0 + << " " << l.zHalf + << " " << l.radiusSensitive + << " " << l.thicknessSensitive + << " " << l.radiusSupport + << " " << l.thicknessSupport + << " " << l.rgap + << " " << l.dphi + << std::endl; + } + + return io; + }; + + /** Simple data structure with key parameters for + * reconstruction of a composite of cylindrical and planar detector + * + * @author + * @date Sep, 02 2024 + * @version $Id: $ + */ + struct CompositeStruct { + /// The half length (z) of the support shell (w/o gap) - 0. if no shell exists. + double zHalfShell; + /// The length of the gap in mm (gap position at z=0). + double gapShell; + /// The inner radius of the support shell. + double rInnerShell; + /// The outer radius of the support shell. + double rOuterShell; + + std::vector<CylindricalStruct::LayerLayout> layersBent; + std::vector<ZPlanarStruct::LayerLayout> layersPlanar; + }; + typedef StructExtension<CompositeStruct> CompositeData; + + static std::ostream& operator<<( std::ostream& io , const CompositeData& d ) { + boost::io::ios_base_all_saver ifs(io); + + io << " -- CompositeData: " << std::scientific << std::endl; + io << " zHalfShell : " << d.zHalfShell << std::endl; + io << " gapShell : " << d.gapShell << std::endl; + io << " rInnerShell : " << d.rInnerShell << std::endl; + io << " rOuterShell : " << d.rOuterShell << std::endl; + + std::vector<CylindricalData::LayerLayout> layersBent = d.layersBent; + + io << " Bent Layers : " << std::endl + << " phi0 zHalf radiusSensitive thicknessSensitive radiusSupport thicknessSupport radialgap deltaphi" << std::endl; + + for(unsigned i=0,N=layersBent.size(); i<N; ++i){ + CylindricalData::LayerLayout l = layersBent[i]; + + io << " " << l.phi0 + << " " << l.zHalf + << " " << l.radiusSensitive + << " " << l.thicknessSensitive + << " " << l.radiusSupport + << " " << l.thicknessSupport + << " " << l.rgap + << " " << l.dphi + << std::endl; + } + + std::vector<ZPlanarData::LayerLayout> layersPlanar = d.layersPlanar; + + io << " Planar Layers : " << std::endl + << " nLadder phi0 nSensors lengthSensor distSupport thickSupport offsetSupport widthSupport zHalfSupport distSense " + << " thickSense offsetSense widthSense zHalfSense" << std::endl; + + for (unsigned i=0,N=layersPlanar.size(); i<N; ++i) { + ZPlanarData::LayerLayout l = layersPlanar[i]; + + io << " " << l.ladderNumber + << " " << l.phi0 + << " " << l.sensorsPerLadder + << " " << l.lengthSensor + << " " << l.distanceSupport + << " " << l.thicknessSupport + << " " << l.offsetSupport + << " " << l.widthSupport + << " " << l.zHalfSupport + << " " << l.distanceSensitive + << " " << l.thicknessSensitive + << " " << l.offsetSensitive + << " " << l.widthSensitive + << " " << l.zHalfSensitive + << std::endl ; + } + + return io; + }; + } +} + +#endif diff --git a/Digitisers/CMakeLists.txt b/Digitisers/CMakeLists.txt index 27da8e2e5b28a073b6a6b8b9ce9d8ddf475e6058..77de90cc5fc9b7100ca2ba66c491765ee90d8a6e 100644 --- a/Digitisers/CMakeLists.txt +++ b/Digitisers/CMakeLists.txt @@ -1,3 +1,4 @@ +add_subdirectory(DigiTool) add_subdirectory(DCHDigi) add_subdirectory(G2CDArbor) add_subdirectory(SimHitMerge) diff --git a/Digitisers/DigiTool/CMakeLists.txt b/Digitisers/DigiTool/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..0628cb818f87434e086f5279000f84772d4dbe05 --- /dev/null +++ b/Digitisers/DigiTool/CMakeLists.txt @@ -0,0 +1,11 @@ +################################################################################ +# Package: DigiTool +################################################################################ + +gaudi_add_header_only_library(DigiTool) + +install(TARGETS DigiTool + EXPORT CEPCSWTargets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + COMPONENT dev) diff --git a/Digitisers/DigiTool/include/DigiTool/IDigiTool.h b/Digitisers/DigiTool/include/DigiTool/IDigiTool.h new file mode 100644 index 0000000000000000000000000000000000000000..84057edc23bc6e9022ecb8e6deae930150ee5428 --- /dev/null +++ b/Digitisers/DigiTool/include/DigiTool/IDigiTool.h @@ -0,0 +1,33 @@ +#ifndef IDigiTool_h +#define IDigiTool_h + +/* + * Description: + * IDigiTool is used to perform digitizer + * + * The interface: + * * Call: peform on digitizer + * + */ + +#include "GaudiKernel/IAlgTool.h" + +namespace edm4hep{ + class SimTrackerHit; + class SimTrackerHitCollection; + class TrackerHitCollection; + class MCRecoTrackerAssociationCollection; +} + +class IDigiTool: virtual public IAlgTool { + public: + + DeclareInterfaceID(IDigiTool, 0, 1); + virtual ~IDigiTool() {} + + virtual StatusCode Call(const edm4hep::SimTrackerHitCollection* simCol, edm4hep::TrackerHitCollection* hitCol, + edm4hep::MCRecoTrackerAssociationCollection* assCol) = 0; + virtual StatusCode Call(edm4hep::SimTrackerHit simhit, edm4hep::TrackerHitCollection* hitCol, + edm4hep::MCRecoTrackerAssociationCollection* assCol) = 0; +}; +#endif diff --git a/Digitisers/SimpleDigi/CMakeLists.txt b/Digitisers/SimpleDigi/CMakeLists.txt index 28c8ea89d61d17c7d123b9d1366e94f49ac7aa29..dfebaaea1daebd39f3176e85452cf74166ba9623 100644 --- a/Digitisers/SimpleDigi/CMakeLists.txt +++ b/Digitisers/SimpleDigi/CMakeLists.txt @@ -4,7 +4,10 @@ gaudi_add_module(SimpleDigi src/TPCDigiAlg.cpp src/voxel.cpp src/CylinderDigiAlg.cpp + src/SiTrackerDigiAlg.cpp + src/SmearDigiTool.cpp LINK GearSvc + DigiTool EventSeeder TrackSystemSvcLib DataHelperLib diff --git a/Digitisers/SimpleDigi/src/CylinderDigiAlg.cpp b/Digitisers/SimpleDigi/src/CylinderDigiAlg.cpp index d3c67e60c89c427ce12a7c4bc27561b1b971f120..7cf254961a68af65350e2ac376d7d7c02468b637 100644 --- a/Digitisers/SimpleDigi/src/CylinderDigiAlg.cpp +++ b/Digitisers/SimpleDigi/src/CylinderDigiAlg.cpp @@ -1,5 +1,7 @@ #include "CylinderDigiAlg.h" +#include "Identifier/CEPCConf.h" + #include "edm4hep/Vector3f.h" #include "DD4hep/Detector.h" @@ -18,10 +20,10 @@ 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"); + declareProperty("SimTrackHitCollection", m_inputColHdls, "Handle of the Input SimTrackerHit collection"); // Output collections - declareProperty("OutputTrackerHitCollection", m_outputColHdls, "Handle of the output TrackerHit collection"); + declareProperty("TrackerHitCollection", m_outputColHdls, "Handle of the output TrackerHit collection"); declareProperty("TrackerHitAssociationCollection", m_assColHdls, "Handle of the Association collection between SimTrackerHit and TrackerHit"); } @@ -74,8 +76,9 @@ StatusCode CylinderDigiAlg::execute(){ 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" ); + int module = m_decoder->get(cellId, "module"); + int sensor = m_decoder->get(cellId, "sensor" ); auto& pos = simhit.getPosition(); auto& mom = simhit.getMomentum(); @@ -93,9 +96,9 @@ StatusCode CylinderDigiAlg::execute(){ 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.setType(1<<CEPCConf::TrkHitTypeBit::CYLINDER); trkHit.addToRawHits(simhit.getObjectID()); - debug() << "Hit " << simhit.id() << ": " << pos << " -> " << trkHit.getPosition() << "s:" << system << " c:" << chamber << " l:" << layer + debug() << "Hit " << simhit.id() << ": " << pos << " -> " << trkHit.getPosition() << "s:" << system << " l:" << layer << " m:" << module << " s:" << sensor << " pt = " << pt << " " << mom.x << " " << mom.y << " " << mom.z << endmsg; auto ass = assVec->create(); diff --git a/Digitisers/SimpleDigi/src/SiTrackerDigiAlg.cpp b/Digitisers/SimpleDigi/src/SiTrackerDigiAlg.cpp new file mode 100644 index 0000000000000000000000000000000000000000..e24a65293aec26e57d5fb1daa57b8c96b81333b4 --- /dev/null +++ b/Digitisers/SimpleDigi/src/SiTrackerDigiAlg.cpp @@ -0,0 +1,70 @@ +#include "SiTrackerDigiAlg.h" + +#include "Identifier/CEPCConf.h" + +#include "edm4hep/Vector3f.h" + +#include "DD4hep/Detector.h" +#include <DD4hep/Objects.h> +#include "DD4hep/DD4hepUnits.h" +#include "DDRec/Vector3D.h" +#include "DDRec/ISurface.h" + +#include "GaudiKernel/INTupleSvc.h" +#include "GaudiKernel/MsgStream.h" +#include "GaudiKernel/IRndmGen.h" +#include "GaudiKernel/IRndmGenSvc.h" +#include "GaudiKernel/RndmGenerators.h" + +DECLARE_COMPONENT(SiTrackerDigiAlg) + +SiTrackerDigiAlg::SiTrackerDigiAlg(const std::string& name, ISvcLocator* svcLoc) +: GaudiAlgorithm(name, svcLoc) { + // Input collections + declareProperty("SimTrackHitCollection", m_inputColHdls, "Handle of the Input SimTrackerHit collection"); + + // Output collections + declareProperty("TrackerHitCollection", m_outputColHdls, "Handle of the output TrackerHit collection"); + declareProperty("TrackerHitAssociationCollection", m_assColHdls, "Handle of the Association collection between SimTrackerHit and TrackerHit"); +} + +StatusCode SiTrackerDigiAlg::initialize() { + m_digiTool = ToolHandle<IDigiTool>(m_digiToolName.value()); + + info() << "DigiTool " << m_digiTool.typeAndName() << " found" << endmsg; + + info() << "SiTrackerDigiAlg::initialized" << endmsg; + return GaudiAlgorithm::initialize(); +} + + +StatusCode SiTrackerDigiAlg::execute(){ + auto trkCol = m_outputColHdls.createAndPut(); + auto assCol = m_assColHdls.createAndPut(); + + const edm4hep::SimTrackerHitCollection* simCol = nullptr; + try { + simCol = m_inputColHdls.get(); + } + catch(GaudiException &e){ + debug() << "Collection " << m_inputColHdls.fullKey() << " is unavailable in event " << m_nEvt << endmsg; + return StatusCode::SUCCESS; + } + if (simCol->size() == 0) return StatusCode::SUCCESS; + debug() << m_inputColHdls.fullKey() << " has SimTrackerHit "<< simCol->size() << endmsg; + + m_digiTool->Call(simCol, trkCol, assCol); + + debug() << "Created " << trkCol->size() << " hits, " + << simCol->size() - trkCol->size() << " hits got dismissed for being out of boundary" + << endmsg; + + m_nEvt++; + + return StatusCode::SUCCESS; +} + +StatusCode SiTrackerDigiAlg::finalize(){ + info() << "Processed " << m_nEvt << " events " << endmsg; + return GaudiAlgorithm::finalize(); +} diff --git a/Digitisers/SimpleDigi/src/SiTrackerDigiAlg.h b/Digitisers/SimpleDigi/src/SiTrackerDigiAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..e25427715f7ea5ea6dcc87ca573c84020f465d62 --- /dev/null +++ b/Digitisers/SimpleDigi/src/SiTrackerDigiAlg.h @@ -0,0 +1,34 @@ +#ifndef SiTrackerDigiAlg_H +#define SiTrackerDigiAlg_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 "DigiTool/IDigiTool.h" + +class SiTrackerDigiAlg : public GaudiAlgorithm{ + public: + + SiTrackerDigiAlg(const std::string& name, ISvcLocator* svcLoc); + + virtual StatusCode initialize(); + virtual StatusCode execute(); + virtual StatusCode finalize(); + +protected: + // Input collections + DataHandle<edm4hep::SimTrackerHitCollection> m_inputColHdls{"VXDCollection", Gaudi::DataHandle::Reader, this}; + // Output collections + DataHandle<edm4hep::TrackerHitCollection> m_outputColHdls{"VXDTrackerHits", Gaudi::DataHandle::Writer, this}; + DataHandle<edm4hep::MCRecoTrackerAssociationCollection> m_assColHdls{"VXDTrackerHitAssociationCollection", Gaudi::DataHandle::Writer, this}; + + Gaudi::Property<std::string> m_digiToolName{this, "DigiTool", "SmearDigiTool/VXD"}; + + ToolHandle<IDigiTool> m_digiTool; + int m_nEvt=0; +}; +#endif diff --git a/Digitisers/SimpleDigi/src/SmearDigiTool.cpp b/Digitisers/SimpleDigi/src/SmearDigiTool.cpp new file mode 100644 index 0000000000000000000000000000000000000000..ea9091807d8838871b938f655afe1e4dc8f7ec13 --- /dev/null +++ b/Digitisers/SimpleDigi/src/SmearDigiTool.cpp @@ -0,0 +1,301 @@ +#include "SmearDigiTool.h" + +#include "DataHelper/TrackerHitHelper.h" +#include "Identifier/CEPCConf.h" +#include "DetInterface/IGeomSvc.h" + +#include "edm4hep/Vector3f.h" + +#include "DD4hep/Detector.h" +#include <DD4hep/Objects.h> +#include "DD4hep/DD4hepUnits.h" +#include "DDRec/ISurface.h" + +#include "GaudiKernel/INTupleSvc.h" +#include "GaudiKernel/MsgStream.h" +#include "GaudiKernel/IRndmGen.h" +#include "GaudiKernel/RndmGenerators.h" + +#include "CLHEP/Vector/ThreeVector.h" +#include "CLHEP/Units/SystemOfUnits.h" + +DECLARE_COMPONENT(SmearDigiTool) + +StatusCode SmearDigiTool::initialize() { + if (m_parameterize) { + if (m_parU.size()!=10 || m_parV.size()!=10) { + fatal() << "parameters number must be 10! now " << m_parU.size() << " for U and " << m_parV.size() << " for V" << endmsg; + return StatusCode::FAILURE; + } + } + else { + if (m_resU.size() != m_resV.size()) { + fatal() << "Inconsistent number of resolutions given for U and V coordinate: " + << "ResolutionU :" << m_resU.size() << " != ResolutionV : " << m_resV.size() + << endmsg; + return StatusCode::FAILURE; + } + } + + m_geosvc = service<IGeomSvc>("GeomSvc"); + if(!m_geosvc){ + error() << "Failed to get the GeomSvc" << endmsg; + return StatusCode::FAILURE; + } + + m_surfaces = m_geosvc->getSurfaceMap(m_detName); + debug() << "Surface map size: " << m_surfaces->size() << endmsg; + if (msgLevel(MSG::VERBOSE)) { + unsigned long old = 0; + for (const auto& pair : *m_surfaces) { + verbose() << pair.first << " : " << pair.second << endmsg; + if (old==pair.first) error() << old << " repeat!" << endmsg; + old = pair.first; + } + } + + debug() << "Readout name: " << m_readoutName << endmsg; + m_decoder = m_geosvc->getDecoder(m_readoutName); + if(!m_decoder){ + error() << "Failed to get the decoder. " << endmsg; + return StatusCode::FAILURE; + } + + m_randSvc = service<IRndmGenSvc>("RndmGenSvc"); + + info() << "initialized" << endmsg; + return StatusCode::SUCCESS; +} + +StatusCode SmearDigiTool::Call(const edm4hep::SimTrackerHitCollection* simCol, edm4hep::TrackerHitCollection* hitCol, + edm4hep::MCRecoTrackerAssociationCollection* assCol) { + for (auto simhit : *simCol) { + StatusCode sc = Call(simhit, hitCol, assCol); + if (sc.isFailure()) return sc; + } + return StatusCode::SUCCESS; +} + +StatusCode SmearDigiTool::Call(edm4hep::SimTrackerHit simhit, edm4hep::TrackerHitCollection* hitCol, edm4hep::MCRecoTrackerAssociationCollection* assCol) { + if (!simhit.isAvailable()) { + error() << "input SimTrackerHit not available!" << endmsg; + return StatusCode::SUCCESS; + } + + auto e = simhit.getEDep(); + if (e <= m_eThreshold) return StatusCode::SUCCESS; + if (m_randSvc->generator(Rndm::Flat(0, 1))->shoot() > m_efficiency) return StatusCode::SUCCESS; + auto t = simhit.getTime(); + + auto cellId = simhit.getCellID(); + int system = m_decoder->get(cellId, "system"); + int side = m_decoder->get(cellId, "side" ); + int layer = m_decoder->get(cellId, "layer" ); + int module = m_decoder->get(cellId, "module"); + int sensor = m_decoder->get(cellId, "sensor"); + + auto& pos = simhit.getPosition(); + auto& mom = simhit.getMomentum(); + + debug() << "Hit " << simhit.id() << " cell: " << cellId << " d:" << system << " l:" << layer << " m:" << module << " s:" << sensor + << " p = " << mom.x << ", " << mom.y << ", " << mom.z << " e:" << e << " t:" << t << endmsg; + + dd4hep::rec::ISurface* surface = nullptr; + auto it = m_surfaces->find(cellId); + if (it != m_surfaces->end()) { + surface = it->second; + if (!surface) { + fatal() << "found surface for cell id " << cellId << ", but NULL" << endmsg; + return StatusCode::FAILURE; + } + } + else { + fatal() << "not found surface for cell id " << cellId << endmsg; + return StatusCode::FAILURE; + } + + dd4hep::rec::Vector3D oldPos(pos.x*dd4hep::mm/CLHEP::mm, pos.y*dd4hep::mm/CLHEP::mm, pos.z*dd4hep::mm/CLHEP::mm); + dd4hep::rec::Vector3D uVec = surface->u(oldPos); + dd4hep::rec::Vector3D vVec = surface->v(oldPos); + float u_direction[2]; + u_direction[0] = uVec.theta(); + u_direction[1] = uVec.phi(); + + float v_direction[2]; + v_direction[0] = vVec.theta(); + v_direction[1] = vVec.phi(); + + debug() << " U: " << uVec << endmsg; + debug() << " V: " << vVec << endmsg; + debug() << " N: " << surface->normal() << endmsg; + debug() << " O: " << surface->origin() << endmsg; + + float resU(0), resV(0), resT(0); + if (!m_parameterize) { + if ((m_resU.size() > 1 && layer >= (int)m_resU.size()) || (m_resV.size() > 1 && layer >= (int)m_resV.size())) { + fatal() << "layer exceeds resolution vector, please check input parameters ResolutionU and ResolutionV" << endmsg; + return StatusCode::FAILURE; + } + + resU = (m_resU.size() > 1 ? m_resU.value().at(layer) : m_resU.value().at(0)); + resV = (m_resV.size() > 1 ? m_resV.value().at(layer) : m_resV.value().at(0)); + } + else { // Riccardo's parameterized model + CLHEP::Hep3Vector momVec(mom[0], mom[1], mom[2]); + CLHEP::Hep3Vector uVecCLHEP(uVec[0], uVec[1], uVec[2]); + CLHEP::Hep3Vector vVecCLHEP(vVec[0], vVec[1], vVec[2]); + const double alpha = uVecCLHEP.azimAngle(momVec, vVecCLHEP); + const double cotanAlpha = 1./tan(alpha); + // TODO: title angle (PI/2), magnetic field (3) + const double tanLorentzAngle = (side==0) ? 0. : 0.053 * 3 * cos(M_PI/2.); + const double x = fabs(-cotanAlpha - tanLorentzAngle); + resU = m_parU[0] + m_parU[1] * x + m_parU[2] * exp(-m_parU[9] * x) * cos(m_parU[3] * x + m_parU[4]) + + m_parU[5] * exp(-0.5 * pow(((x - m_parU[6]) / m_parU[7]), 2)) + m_parU[8] * pow(x, 0.5); + + const double beta = vVecCLHEP.azimAngle(momVec, uVecCLHEP); + const double cotanBeta = 1./tan(beta); + const double y = fabs(-cotanBeta); + resV = m_parV[0] + m_parV[1] * y + m_parV[2] * exp(-m_parV[9] * y) * cos(m_parV[3] * y + m_parV[4]) + + m_parV[5] * exp(-0.5 * pow(((y - m_parV[6]) / m_parV[7]), 2)) + m_parV[8] * pow(y, 0.5); + } + resU *= dd4hep::mm/CLHEP::mm; + resV *= dd4hep::mm/CLHEP::mm; + // parameterize only for position now, todo + resT = (m_resT.size() > 1 ? m_resT.value().at(layer) : m_resT.value().at(0)); + // from ps (input unit) to ns (record unit, Geant4) + resT *= CLHEP::ps/CLHEP::ns; + debug() << " --- will smear hit with resU = " << resU/dd4hep::mm << " resV = " << resV/dd4hep::mm << " resT = " << resT << endmsg; + + auto& typeSurface = surface->type(); + debug() << " Surface id: " << surface->id() << " type:" << endmsg; + debug() << " " << typeSurface << endmsg; + verbose() << " count: " << m_surfaces->count(cellId) << " inner: " << surface->innerMaterial().name() + << " outer: " << surface->outerMaterial().name() << endmsg; + + if (typeSurface.isPlane() || typeSurface.isCylinder()) { + dd4hep::rec::Vector2D localPoint = surface->globalToLocal(oldPos); + dd4hep::rec::Vector3D local3D(localPoint.u(), localPoint.v(), 0); + // A small check, if the hit is in the boundaries: + if (!surface->insideBounds(oldPos, 1e100)) { + error() << " global: (" << oldPos.x() << " " << oldPos.y() << " " << oldPos.z() + << ") local: (" << localPoint.u() << ", " << localPoint.v() << " )" + << " is not within boundaries. Hit is skipped." << endmsg; + return StatusCode::SUCCESS; + } + dd4hep::rec::Vector3D globalPointSmeared;//CLHEP::Hep3Vector globalPoint(pos[0],pos[1],pos[2]); + dd4hep::rec::Vector2D localPointSmeared; + + debug() << std::setprecision(8) << " Before smearing global: (" << pos[0]/CLHEP::mm << " " << pos[1]/CLHEP::mm << " " << pos[2]/CLHEP::mm << ") " + << "local: (" << localPoint.u()/dd4hep::mm << " " << localPoint.v()/dd4hep::mm << ")" << endmsg; + + unsigned tries = 0; + + bool accept_hit = false; + + // Now try to smear the hit and make sure it is still on the surface + while (tries < 100) { + if (tries > 0) { + debug() << "retry smearing for side" << side << " layer"<< layer<< " module" << module + << " sensor" << sensor << " : retries " << tries << endmsg; + } + + double du = m_randSvc->generator(Rndm::Gauss(0, resU))->shoot(); + double dv = m_randSvc->generator(Rndm::Gauss(0, resV))->shoot(); + localPointSmeared.u() = localPoint.u() + du; + localPointSmeared.v() = localPoint.v() + dv; + + dd4hep::rec::Vector3D local3DSmeared(localPointSmeared.u(), localPointSmeared.v(), 0); + globalPointSmeared = surface->localToGlobal(localPointSmeared); + + //check if hit is in boundaries + if (surface->insideBounds(globalPointSmeared) && fabs(du) <= m_maxPull*resU && fabs(dv) <= m_maxPull*resV) { + accept_hit = true; + break; + } + tries++; + } + + if (accept_hit == false) { + debug() << "hit could not be smeared within ladder after 100 tries: hit dropped" << endmsg; + return StatusCode::SUCCESS; + } + + // for 1D strip measurements: set v to 0! Only the measurement in u counts! + if(m_isStrip || (resU != 0 && resV == 0)) localPointSmeared.v() = 0. ; + // convert back to global position for TrackerHit + globalPointSmeared = surface->localToGlobal(localPointSmeared); + + debug() << " After smearing global: (" + << globalPointSmeared.x()/dd4hep::mm <<" "<< globalPointSmeared.y()/dd4hep::mm <<" "<< globalPointSmeared.z()/dd4hep::mm << ") " + << "local: (" + << localPointSmeared.u()/dd4hep::mm << " " << localPointSmeared.v()/dd4hep::mm << ")" << endmsg; + + auto outhit = hitCol->create(); + outhit.setCellID(cellId); + + edm4hep::Vector3d smearedPos(globalPointSmeared.x()*CLHEP::mm/dd4hep::mm, + globalPointSmeared.y()*CLHEP::mm/dd4hep::mm, + globalPointSmeared.z()*CLHEP::mm/dd4hep::mm); + outhit.setPosition(smearedPos); + // recover CLHEP/Geant4 unit + resU /= dd4hep::mm/CLHEP::mm; + resV /= dd4hep::mm/CLHEP::mm; + + std::bitset<32> type; + if (typeSurface.isPlane() && m_usePlanarTag) { + std::array<float, 6> cov; + cov[0] = u_direction[0]; + cov[1] = u_direction[1]; + cov[2] = resU; + cov[3] = v_direction[0]; + cov[4] = v_direction[1]; + cov[5] = resV; + outhit.setCovMatrix(cov); + + type.set(CEPCConf::TrkHitTypeBit::PLANAR); + } + else if (typeSurface.isPlane()) { + outhit.setCovMatrix(CEPC::ConvertToCovXYZ(resU, u_direction[0], u_direction[1], resV, v_direction[0], v_direction[1])); + } + else if (typeSurface.isCylinder()) { + outhit.setCovMatrix(std::array<float, 6>{resU*resU/2., 0, resU*resU/2, 0, 0, resV*resV}); + type.set(CEPCConf::TrkHitTypeBit::CYLINDER); + } + + if(m_isStrip || (resU != 0 && resV == 0)){ + type.set(CEPCConf::TrkHitTypeBit::ONE_DIMENSIONAL); + } + outhit.setType((int)type.to_ulong()); + outhit.setEDep(e); + float dt = m_randSvc->generator(Rndm::Gauss(0, resT))->shoot(); + outhit.setTime(simhit.getTime() + dt); + // make the relation + auto ass = assCol->create(); + + float weight = 1.0; + + debug() <<" Set relation between " + << " sim hit " << simhit.id() + << " to tracker hit " << outhit.id() + << " with a weight of " << weight + << endmsg; + + outhit.addToRawHits(simhit.getObjectID()); + ass.setSim(simhit); + ass.setRec(outhit); + ass.setWeight(weight); + + debug() << "-------------------------------------------------------" << endmsg; + } + else { + fatal() << "Not plane and cylinder: " << typeSurface << endmsg; + return StatusCode::FAILURE; + } + + return StatusCode::SUCCESS; +} + +StatusCode SmearDigiTool::finalize(){ + StatusCode sc; + return sc; +} diff --git a/Digitisers/SimpleDigi/src/SmearDigiTool.h b/Digitisers/SimpleDigi/src/SmearDigiTool.h new file mode 100644 index 0000000000000000000000000000000000000000..2f8f9fa112f7f1bfdab8f900f51487896fd7ce20 --- /dev/null +++ b/Digitisers/SimpleDigi/src/SmearDigiTool.h @@ -0,0 +1,61 @@ +#ifndef SmearDigiTool_h +#define SmearDigiTool_h + +#include "GaudiKernel/AlgTool.h" +#include "GaudiKernel/IRndmGenSvc.h" + +#include "DigiTool/IDigiTool.h" +#include "DetInterface/IGeomSvc.h" + +#include "edm4hep/SimTrackerHitCollection.h" +#include "edm4hep/TrackerHitCollection.h" +#include "edm4hep/MCRecoTrackerAssociationCollection.h" + +#include "DDRec/DetectorData.h" +#include "DDRec/SurfaceManager.h" + +#include <vector> + +class SmearDigiTool : public extends<AlgTool, IDigiTool> { + public: + using extends::extends; + //SmearDigiTool(void* p) { m_pAlgUsing=p; }; + + virtual StatusCode Call(const edm4hep::SimTrackerHitCollection* simCol, edm4hep::TrackerHitCollection* hitCol, + edm4hep::MCRecoTrackerAssociationCollection* assCol) override; + virtual StatusCode Call(edm4hep::SimTrackerHit simhit, edm4hep::TrackerHitCollection* hitCol, + edm4hep::MCRecoTrackerAssociationCollection* assCol) override; + + StatusCode initialize() override; + StatusCode finalize() override; + + private: + Gaudi::Property<std::string> m_detName{this, "DetName", "VXD"}; + Gaudi::Property<std::string> m_readoutName{this, "Readout", "VXDCollection"}; + + Gaudi::Property<float> m_eThreshold{this, "EnergyThreshold", 0}; + Gaudi::Property<float> m_efficiency{this, "Efficiency", 1}; + + // resolution in direction of u - either one per layer or one for all layers + Gaudi::Property<std::vector<float> > m_resU{this, "ResolutionU", {0.004}}; + // resolution in direction of v - either one per layer or one for all layers + Gaudi::Property<std::vector<float> > m_resV{this, "ResolutionV", {0.004}}; + // resolution of time - either one per layer or one for all layers, ps as unit + Gaudi::Property<std::vector<float> > m_resT{this, "ResolutionT", {0.0}}; + // whether hits are 1D strip hits + Gaudi::Property<bool> m_isStrip{this, "IsStrip", false}; + // whether use Planar tag for type and cov, if true, CEPCConf::TrkHitTypeBit::PLANAR bit is set as true + // cov[0]=thetaU, cov[1]=phiU, cov[2]=resU, cov[0]=thetaV, cov[1]=phiV, cov[2]=resV + Gaudi::Property<bool> m_usePlanarTag{this, "UsePlanarTag", true}; + Gaudi::Property<float> m_maxPull{this, "PullCutToRetry", 1000.}; + Gaudi::Property<bool> m_parameterize{this, "ParameterizeResolution", false}; + Gaudi::Property<std::vector<float> > m_parU{this, "ParametersU", {0}}; + Gaudi::Property<std::vector<float> > m_parV{this, "ParametersV", {0}}; + + SmartIF<IRndmGenSvc> m_randSvc; + SmartIF<IGeomSvc> m_geosvc; + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + const dd4hep::rec::SurfaceMap* m_surfaces; + //void* m_pAlgUsing = nullptr; +}; +#endif diff --git a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp index 9515f3f3274a9097439db791871c5b102294c145..7921fd3e1fa5d27ae41d235c90298b58c1ad221f 100644 --- a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp +++ b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp @@ -1,4 +1,6 @@ #include "SiliconTrackingAlg.h" + +#include "Identifier/CEPCConf.h" #include "GearSvc/IGearSvc.h" #include "EventSeeder/IEventSeeder.h" #include "TrackSystemSvc/ITrackSystemSvc.h" @@ -506,12 +508,12 @@ int SiliconTrackingAlg::InitialiseFTD() { gear::Vector3D Z(0.0,0.0,1.0); const float eps = 1.0e-07; - // V must be the global z axis + // V must be the global z axis if( fabs(V.dot(Z)) > eps ) { error() << "SiliconTrackingAlg: FTD Hit measurment vectors V is not in the global X-Y plane. \n\n exit(1) called from file " << __FILE__ << " and line " << __LINE__ << endmsg; exit(1); } - + if( fabs(U.dot(Z)) > eps ) { error() << "SiliconTrackingAlg: FTD Hit measurment vectors U is not in the global X-Y plane. \n\n exit(1) called from file " << __FILE__ << " and line " << __LINE__ << endmsg; exit(1); @@ -720,29 +722,41 @@ int SiliconTrackingAlg::InitialiseVTX() { for (int ielem=0; ielem<nelem; ++ielem) { //for(auto hit : *hitVTXCol){ auto hit = hitVTXCol->at(ielem); - //gear::Vector3D U(1.0,hit->getU()[1],hit->getU()[0],gear::Vector3D::spherical); - //gear::Vector3D V(1.0,hit->getV()[1],hit->getV()[0],gear::Vector3D::spherical); - gear::Vector3D U(1.0,hit.getCovMatrix()[1],hit.getCovMatrix()[0],gear::Vector3D::spherical); - gear::Vector3D V(1.0,hit.getCovMatrix()[4],hit.getCovMatrix()[3],gear::Vector3D::spherical); - gear::Vector3D Z(0.0,0.0,1.0); - //debug() << "covMatrix : " << hit->getCovMatrix()[0] << " " << hit->getCovMatrix()[1] << endmsg; - const float eps = 1.0e-07; - // V must be the global z axis - if( fabs(1.0 - V.dot(Z)) > eps ) { - error() << "SiliconTrackingAlg: VXD Hit measurment vectors V is not equal to the global Z axis. \n\n exit(1) called from file " << __FILE__ << " and line " << __LINE__ << endmsg; - exit(1); + + int type = hit.getType(); + debug() << "type = " << type << endmsg; + double resRPhi, resZ; + if (UTIL::BitSet32(type)[CEPCConf::TrkHitTypeBit::PLANAR]) { + //gear::Vector3D U(1.0,hit->getU()[1],hit->getU()[0],gear::Vector3D::spherical); + //gear::Vector3D V(1.0,hit->getV()[1],hit->getV()[0],gear::Vector3D::spherical); + gear::Vector3D U(1.0,hit.getCovMatrix()[1],hit.getCovMatrix()[0],gear::Vector3D::spherical); + gear::Vector3D V(1.0,hit.getCovMatrix()[4],hit.getCovMatrix()[3],gear::Vector3D::spherical); + gear::Vector3D Z(0.0,0.0,1.0); + //debug() << "covMatrix : " << hit->getCovMatrix()[0] << " " << hit->getCovMatrix()[1] << endmsg; + const float eps = 1.0e-07; + // V must be the global z axis + if( fabs(1.0 - V.dot(Z)) > eps ) { + error() << "SiliconTrackingAlg: VXD Hit measurment vectors V is not equal to the global Z axis. \n\n exit(1) called from file " << __FILE__ << " and line " << __LINE__ << endmsg; + exit(1); + } + + if( fabs(U.dot(Z)) > eps ) { + error() << "SiliconTrackingAlg: VXD Hit measurment vectors U is not in the global X-Y plane. \n\n exit(1) called from file " << __FILE__ << " and line " << __LINE__ << endmsg; + exit(1); + } + // SJA:FIXME: just use planar res for now + resRPhi = hit.getCovMatrix()[2]; + resZ = hit.getCovMatrix()[5]; } - - if( fabs(U.dot(Z)) > eps ) { - error() << "SiliconTrackingAlg: VXD Hit measurment vectors U is not in the global X-Y plane. \n\n exit(1) called from file " << __FILE__ << " and line " << __LINE__ << endmsg; - exit(1); + else { + resRPhi = sqrt(hit.getCovMatrix()[0]+hit.getCovMatrix()[2]); + resZ = sqrt(hit.getCovMatrix()[5]); } TrackerHitExtended * hitExt = new TrackerHitExtended(hit); //debug() << "Saved TrackerHit id in TrackerHitExtended " << ielem << ": " << hitExt->getTrackerHit().id() << std::endl; - - // SJA:FIXME: just use planar res for now - hitExt->setResolutionRPhi(hit.getCovMatrix()[2]); - hitExt->setResolutionZ(hit.getCovMatrix()[5]); + + hitExt->setResolutionRPhi(resRPhi); + hitExt->setResolutionZ(resZ); // set type is now only used in one place where it is set to 0 to reject hits from a fit, set to INT_MAX to try and catch any missuse hitExt->setType(int(INT_MAX)); @@ -3004,6 +3018,9 @@ StatusCode SiliconTrackingAlg::setupGearGeom(){ pVXDDetMain = &gearMgr->getVXDParameters(); pVXDLayerLayout = &(pVXDDetMain->getVXDLayerLayout()); _nLayersVTX = pVXDLayerLayout->getNLayers(); + + const std::vector<int> ids = pVXDDetMain->getIntVals("VTXLayerIds"); + _nLayersVTX += ids.size(); } catch( gear::UnknownParameterException& e){ diff --git a/Service/GearSvc/CMakeLists.txt b/Service/GearSvc/CMakeLists.txt index 686ad1cd1570e3aa9ebca4af781aebb4c1e07c9a..13b502dd4e5ae5544afc5dd27733089e616cb820 100644 --- a/Service/GearSvc/CMakeLists.txt +++ b/Service/GearSvc/CMakeLists.txt @@ -8,7 +8,8 @@ gaudi_add_module(GearSvcPlugins ${GEAR_LIBRARIES} ${DD4hep_COMPONENT_LIBRARIES} DetInterface - DetSegmentation + DetSegmentation + Identifier ) target_include_directories(GearSvcPlugins diff --git a/Service/GearSvc/src/GearSvc.cpp b/Service/GearSvc/src/GearSvc.cpp index 3facee7a602188035f4d96abe6d7dc5d9b27fe1f..e4f4976e187cbf9c3792a5817cccebb193c52bd0 100644 --- a/Service/GearSvc/src/GearSvc.cpp +++ b/Service/GearSvc/src/GearSvc.cpp @@ -1,6 +1,7 @@ #include "GearSvc.h" #include "DetInterface/IGeomSvc.h" #include "DetSegmentation/GridDriftChamber.h" +#include "Identifier/CEPCDetectorData.h" #include "gearxml/GearXML.h" #include "gearimpl/GearMgrImpl.h" @@ -22,17 +23,6 @@ #include "DD4hep/DD4hepUnits.h" #include "CLHEP/Units/SystemOfUnits.h" -struct helpLayer { - double distance =0; - double offset =0; - double thickness =0; - double length =0; - double width =0; - double radLength =0; - double z =0; - double foam_spacer_radLength =0; -}; - static const double deg_to_rad = dd4hep::degree/CLHEP::rad; static const double rad_to_deg = dd4hep::rad/CLHEP::degree; @@ -91,13 +81,18 @@ StatusCode GearSvc::initialize() if(it->first=="Tube"||it->first=="BeamPipe"){ sc = convertBeamPipe(sub); } - else if(it->first=="VXD"){ + else if(it->first=="VXD" || it->first=="VTX"){ sc = convertVXD(sub); + if (sc.isRecoverable()) sc = convertStitching(sub); + if (sc.isRecoverable()) sc = convertComposite(sub); + if (!sc.isSuccess()) { + error() << it->first << " extension not read" << endmsg; + } } - else if(it->first=="FTD"){ + else if(it->first=="FTD" || it->first=="ITKEndcap"){ sc = convertFTD(sub); } - else if(it->first=="SIT"){ + else if(it->first=="SIT" || it->first=="ITKBarrel"){ sc = convertSIT(sub); } else if(it->first=="TPC"){ @@ -106,7 +101,7 @@ StatusCode GearSvc::initialize() else if(it->first=="DriftChamber"){ sc = convertDC(sub); } - else if(it->first=="SET"){ + else if(it->first=="SET" || it->first=="OTKBarrel"){ sc = convertSET(sub); } else if(it->first=="EcalBarrel"||it->first=="EcalEndcap"||it->first=="EcalPlug"|| @@ -197,67 +192,41 @@ StatusCode GearSvc::convertBeamPipe(dd4hep::DetElement& pipe){ StatusCode GearSvc::convertVXD(dd4hep::DetElement& vxd){ StatusCode sc; - //fucd: another method to obtain parameters, but not fully for KalDet + // always use extension data dd4hep::rec::ZPlanarData* vxdData = nullptr; - bool extensionDataValid = true; try{ vxdData = vxd.extension<dd4hep::rec::ZPlanarData>(); } catch(std::runtime_error& e){ - extensionDataValid = false; info() << e.what() << " " << vxdData << endmsg; + return StatusCode::RECOVERABLE; } if(vxdData){ int vxdType = gear::ZPlanarParameters::CMOS; - gear::ZPlanarParametersImpl* gearVXD = new gear::ZPlanarParametersImpl( vxdType, vxdData->rInnerShell/dd4hep::mm, vxdData->rOuterShell/dd4hep::mm, - vxdData->zHalfShell/dd4hep::mm , vxdData->gapShell/dd4hep::mm , 0. ) ; + gear::ZPlanarParametersImpl* gearVXD = new gear::ZPlanarParametersImpl(vxdType, vxdData->rInnerShell/dd4hep::mm, vxdData->rOuterShell/dd4hep::mm, + vxdData->zHalfShell/dd4hep::mm, vxdData->gapShell/dd4hep::mm, 0.); for(unsigned i=0,n=vxdData->layers.size() ; i<n; ++i){ const dd4hep::rec::ZPlanarData::LayerLayout& l = vxdData->layers[i]; // FIXME set rad lengths to 0 -> need to get from dd4hep .... gearVXD->addLayer(l.ladderNumber, l.phi0, - l.distanceSupport/dd4hep::mm, l.offsetSupport/dd4hep::mm, l.thicknessSupport/dd4hep::mm, l.zHalfSupport/dd4hep::mm, l.widthSupport/dd4hep::mm, 0., - l.distanceSensitive/dd4hep::mm, l.offsetSensitive/dd4hep::mm, l.thicknessSensitive/dd4hep::mm, l.zHalfSensitive/dd4hep::mm, l.widthSensitive/dd4hep::mm, 0.); + l.distanceSupport/dd4hep::mm, l.offsetSupport/dd4hep::mm, l.thicknessSupport/dd4hep::mm, + l.zHalfSupport/dd4hep::mm, l.widthSupport/dd4hep::mm, 0., + l.distanceSensitive/dd4hep::mm, l.offsetSensitive/dd4hep::mm, l.thicknessSensitive/dd4hep::mm, + l.zHalfSensitive/dd4hep::mm, l.widthSensitive/dd4hep::mm, 0.); } m_gearMgr->setVXDParameters(gearVXD); - dd4hep::rec::MaterialManager matMgr( dd4hep::Detector::getInstance().world().volume() ) ; const dd4hep::rec::ZPlanarData::LayerLayout& l = vxdData->layers[0] ; double offset = l.offsetSupport; - //dd4hep::rec::Vector3D a( l.distanceSensitive + l.thicknessSensitive, l.phi0 , 0. , dd4hep::rec::Vector3D::cylindrical ) ; - //dd4hep::rec::Vector3D b( l.distanceSupport + l.thicknessSupport, l.phi0 , 0. , dd4hep::rec::Vector3D::cylindrical ) ; - dd4hep::rec::Vector3D a( l.distanceSensitive + l.thicknessSensitive, l.offsetSupport, 2.*dd4hep::mm); - dd4hep::rec::Vector3D b( l.distanceSupport + l.thicknessSupport, l.offsetSupport, 2.*dd4hep::mm); - const dd4hep::rec::MaterialVec& materials = matMgr.materialsBetween( a , b ) ; - dd4hep::rec::MaterialData mat = ( materials.size() > 1 ? matMgr.createAveragedMaterial( materials ) : materials[0].first ) ; - - std::cout << " ####### found materials between points : " << a << " and " << b << " : " ; - for( unsigned i=0,n=materials.size();i<n;++i){ - std::cout << materials[i].first.name() << "[" << materials[i].second << "], " ; - } - std::cout << std::endl ; - std::cout << " averaged material : " << mat << std::endl ; - gear::SimpleMaterialImpl* VXDSupportMaterial = new gear::SimpleMaterialImpl("VXDSupportMaterial", mat.A(), mat.Z(), - mat.density()/(dd4hep::kg/(dd4hep::g*dd4hep::m3)), - mat.radiationLength()/dd4hep::mm, - mat.interactionLength()/dd4hep::mm); + dd4hep::rec::Vector3D a( l.distanceSensitive + l.thicknessSensitive, offset, 2.*dd4hep::mm); + dd4hep::rec::Vector3D b( l.distanceSupport + l.thicknessSupport, offset, 2.*dd4hep::mm); + gear::SimpleMaterialImpl* VXDSupportMaterial = CreateGearMaterial(a, b, "VXDSupportMaterial"); m_gearMgr->registerSimpleMaterial(VXDSupportMaterial); if (vxdData->rOuterShell>vxdData->rInnerShell) { dd4hep::rec::Vector3D a1( vxdData->rInnerShell, 0, 2.*dd4hep::mm); dd4hep::rec::Vector3D b1( vxdData->rOuterShell, 0, 2.*dd4hep::mm); - const dd4hep::rec::MaterialVec& materials1 = matMgr.materialsBetween( a1 , b1 ) ; - dd4hep::rec::MaterialData mat1 = ( materials1.size() > 1 ? matMgr.createAveragedMaterial( materials1 ) : materials1[0].first ) ; - - std::cout << " ####### found materials between points : " << a1 << " and " << b1 << " : " ; - for( unsigned i=0,n=materials1.size();i<n;++i){ - std::cout << materials1[i].first.name() << "[" << materials1[i].second << "], " ; - } - std::cout << std::endl ; - std::cout << " averaged material : " << mat1 << std::endl ; - gear::SimpleMaterialImpl* VXDShellMaterial = new gear::SimpleMaterialImpl("VXDShellMaterial", mat1.A(), mat1.Z(), - mat1.density()/(dd4hep::kg/(dd4hep::g*dd4hep::m3)), - mat1.radiationLength()/dd4hep::mm, - mat1.interactionLength()/dd4hep::mm); + gear::SimpleMaterialImpl* VXDShellMaterial = CreateGearMaterial(a1, b1, "VXDShellMaterial"); m_gearMgr->registerSimpleMaterial(VXDShellMaterial); } @@ -270,418 +239,133 @@ StatusCode GearSvc::convertVXD(dd4hep::DetElement& vxd){ << thisLayer.distanceSensitive/dd4hep::mm << "," << thisLayer.offsetSensitive/dd4hep::mm << "," << thisLayer.thicknessSensitive/dd4hep::mm << "," << thisLayer.zHalfSensitive/dd4hep::mm << "," << thisLayer.widthSensitive/dd4hep::mm << ",NULL" << endmsg; } - return sc; } + return StatusCode::SUCCESS; +} - std::vector<helpLayer> helpSensitives; - std::vector<helpLayer> helpLadders; - std::vector<int> helpNumberLadders; - std::vector<double> helpPhi0; - int helpCount=0; - int type=0; - double shellInnerRadius, shellOuterRadius, shellHalfLength, gap, shellRadLength; - int nLadders=0; - double phi0=0; - helpLayer thisLadder; - double beryllium_ladder_block_length=0,end_electronics_half_z=0,side_band_electronics_width=0; - double rAlu=0, drAlu=0, rSty=0, drSty=0, dzSty=0, rInner=0, aluEndcapZ=0, aluHalfZ=0, alu_RadLen=0, Cryostat_dEdx=0; - double VXDSupportDensity=0, VXDSupportZeff=0, VXDSupportAeff=0, VXDSupportRadLen=0, VXDSupportIntLen=0; - double styDensity=0, styZeff=0, styAeff=0, styRadLen=0, styIntLen=0; - dd4hep::Volume vxd_vol = vxd.volume(); - for(int i=0;i<vxd_vol->GetNdaughters();i++){ - TGeoNode* daughter = vxd_vol->GetNode(i); - std::string nodeName = daughter->GetName(); - if(nodeName=="VXD_support_assembly_0"){ - TGeoNode* shell = FindNode(daughter, "SupportShell"); - if(shell){ - const TGeoShape* shape_shell = shell->GetVolume()->GetShape(); - //fucd: IsA() method does not always work for TGeoTube, sometimes, strange? - //if(shape_shell->IsA()==TGeoTube::Class()){ - if(shape_shell->TestShapeBit(TGeoTube::kGeoTube)){ - const TGeoTube* tube = (const TGeoTube*) shape_shell; - shellInnerRadius = tube->GetRmin()*CLHEP::cm; - shellOuterRadius = tube->GetRmax()*CLHEP::cm; - shellHalfLength = tube->GetDz()*CLHEP::cm; - } - else{ - error() << shell->GetName() << " is not a TGeoTube! Shape bits = " << shape_shell->TestShapeBits(0xFFFFFFFF) << endmsg; - } - TGeoMaterial* mat = shell->GetMedium()->GetMaterial(); - shellRadLength = mat->GetRadLen()*CLHEP::cm; - } - TGeoNode* block = FindNode(daughter, "BerylliumAnnulusBlock"); - if(block){ - const TGeoShape* shape_block = block->GetVolume()->GetShape(); - if(shape_block->IsA()==TGeoBBox::Class()){ - const TGeoBBox* box = (const TGeoBBox*) shape_block; - beryllium_ladder_block_length = box->GetDY()*CLHEP::cm; - } - else{ - error() << block->GetName() << " is not a TGeoTube! Shape bits = " << shape_block->TestShapeBits(0xFFFFFFFF) << endmsg; - } - } - TGeoNode* skin = FindNode(daughter, "CryostatAluSkinBarrel"); - if(skin){ - const TGeoShape* shape_skin = skin->GetVolume()->GetShape(); - if(shape_skin->TestShapeBit(TGeoTube::kGeoTube)){ - const TGeoTube* tube = (const TGeoTube*) shape_skin; - rAlu = tube->GetRmin()*CLHEP::cm; - drAlu = tube->GetRmax()*CLHEP::cm - rAlu; - aluHalfZ = tube->GetDz()*CLHEP::cm; - } - else{ - error() << skin->GetName() << " is not a TGeoTube! Shape bits = " << shape_skin->TestShapeBits(0xFFFFFFFF) << endmsg; - } - } - TGeoNode* foam = FindNode(daughter, "CryostatFoamBarrel"); - if(foam){ - const TGeoShape* shape_foam = foam->GetVolume()->GetShape(); - if(shape_foam->TestShapeBit(TGeoTube::kGeoTube)){ - const TGeoTube* tube = (const TGeoTube*) shape_foam; - rSty = tube->GetRmin()*CLHEP::cm; - drSty = tube->GetRmax()*CLHEP::cm - rSty; - dzSty = tube->GetDz()*CLHEP::cm; - } - else{ - error() << foam->GetName() << " is not a TGeoTube! Shape bits = " << shape_foam->TestShapeBits(0xFFFFFFFF) << endmsg; - } - TGeoMaterial* mat = foam->GetMedium()->GetMaterial(); - double Zeff = 0, ZAeff = 0; - for(int iEle = 0; iEle<mat->GetNelements(); iEle++){ - double A, Z, w; - mat->GetElementProp(A,Z,w,iEle); - Zeff += Z*w; - ZAeff += Z/A*w; - //std::cout << std::setprecision(16) << Z << " " << A << " " << w << std::endl; - } - styZeff = Zeff; - styAeff = Zeff/ZAeff; - styRadLen = mat->GetRadLen()*CLHEP::cm; - styIntLen = mat->GetIntLen()*CLHEP::cm; - styDensity = mat->GetDensity(); - } - TGeoNode* skinEnd = FindNode(daughter, "CryostatAluSkinEndPlateInner"); - if(skinEnd){ - const TGeoShape* shape_skinEnd = skinEnd->GetVolume()->GetShape(); - if(shape_skinEnd->TestShapeBit(TGeoTube::kGeoTube)){ - const TGeoTube* tube = (const TGeoTube*) shape_skinEnd; - rInner = tube->GetRmin()*CLHEP::cm; - double rmax = tube->GetRmax()*CLHEP::cm; - drAlu = tube->GetDz()*CLHEP::cm*2; - } - else{ - error() << skinEnd->GetName() << " is not a TGeoTube! Shape bits = " << shape_skinEnd->TestShapeBits(0xFFFFFFFF) << endmsg; - } - } - TGeoNode* shellEnd = FindNode(daughter, "EndPlateShell_outer"); - if(shellEnd){ - const TGeoShape* shape_shellEnd = shellEnd->GetVolume()->GetShape(); - if(shape_shellEnd->TestShapeBit(TGeoTube::kGeoTube)){ - const TGeoTube* tube = (const TGeoTube*) shape_shellEnd; - double rmin = tube->GetRmin()*CLHEP::cm; - double rmax = tube->GetRmax()*CLHEP::cm; - double zhalf = tube->GetDz()*CLHEP::cm; - } - else{ - error() << shellEnd->GetName() << " is not a TGeoTube! Shape bits = " << shape_shellEnd->TestShapeBits(0xFFFFFFFF) << endmsg; - } - } - +StatusCode GearSvc::convertStitching(dd4hep::DetElement& vtx){ + dd4hep::rec::CylindricalData* vtxData = nullptr; + try{ + vtxData = vtx.extension<dd4hep::rec::CylindricalData>(); + } + catch(std::runtime_error& e){ + warning() << e.what() << " " << vtxData << endmsg; + return StatusCode::RECOVERABLE; + } + if(vtxData){ + int vtxType = gear::ZPlanarParameters::CMOS; + gear::ZPlanarParametersImpl* gearVTX = new gear::ZPlanarParametersImpl(vtxType, vtxData->rInnerShell/dd4hep::mm, vtxData->rOuterShell/dd4hep::mm, + vtxData->zHalfShell/dd4hep::mm, vtxData->gapShell/dd4hep::mm, 0.); + std::vector<int> ids; + std::vector<double> zhalfs, rsens, tsens, rsups, tsups, phi0s, rgaps, dphis; + + for (unsigned i=0,n=vtxData->layers.size(); i<n; ++i) { + const dd4hep::rec::CylindricalData::LayerLayout& l = vtxData->layers[i]; + + ids.push_back(l.id); + zhalfs.push_back(l.zHalf/dd4hep::mm); + rsens.push_back(l.radiusSensitive/dd4hep::mm); + tsens.push_back(l.thicknessSensitive/dd4hep::mm); + rsups.push_back(l.radiusSupport/dd4hep::mm); + tsups.push_back(l.thicknessSupport/dd4hep::mm); + phi0s.push_back(l.phi0); + rgaps.push_back(l.rgap/dd4hep::mm); + dphis.push_back(l.dphi); } - else if(nodeName=="layer_assembly_0_1"){ - if(TGeoNode* side_band = FindNode(daughter, "ElectronicsBand")){ - const TGeoShape* shape_band = side_band->GetVolume()->GetShape(); - if(shape_band->IsA()==TGeoBBox::Class()){ - const TGeoBBox* box = (const TGeoBBox*) shape_band; - side_band_electronics_width = box->GetDX()*CLHEP::cm*2; - //info() << "fucd: "<< box->GetDX() << " " << box->GetDY() << " " << box->GetDZ() <<endmsg; - } - else{ - error() << "ElectronicsBand is not a TGeoBBox!!!"<< endmsg; - } - } - if(TGeoNode* end = FindNode(daughter, "ElectronicsEnd")){ - const TGeoShape* shape_end = end->GetVolume()->GetShape(); - if(shape_end->IsA()==TGeoBBox::Class()){ - const TGeoBBox* box = (const TGeoBBox*) shape_end; - end_electronics_half_z = box->GetDY()*CLHEP::cm; - //info() << "fucd: " << box->GetDX() << " " << box->GetDY() << " " << box->GetDZ() << endmsg; - } - else{ - error() << "ElectronicsEnd is not a TGeoBBox!!!"<< endmsg; - } - } + gearVTX->setIntVals("VTXLayerIds", ids); + gearVTX->setDoubleVals("VTXLayerHalfLengths", zhalfs); + gearVTX->setDoubleVals("VTXLayerSensitiveRadius", rsens); + gearVTX->setDoubleVals("VTXLayerSensitiveThickness", tsens); + gearVTX->setDoubleVals("VTXLayerSupperRadius", rsups); + gearVTX->setDoubleVals("VTXLayerSupperThickness", tsups); + gearVTX->setDoubleVals("VTXLayerPhi0", phi0s); + gearVTX->setDoubleVals("VTXLayerRadialGap", rgaps); + gearVTX->setDoubleVals("VTXLayerDeltaPhi", dphis); + + m_gearMgr->setVXDParameters(gearVTX); + + const dd4hep::rec::CylindricalData::LayerLayout& l = vtxData->layers[0] ; + dd4hep::rec::Vector3D a( l.radiusSupport, l.phi0 , 0. , dd4hep::rec::Vector3D::cylindrical ) ; + dd4hep::rec::Vector3D b( l.radiusSupport + l.thicknessSupport, l.phi0 , 0. , dd4hep::rec::Vector3D::cylindrical ) ; + gear::SimpleMaterialImpl* VXDSupportMaterial = CreateGearMaterial(a, b, "VXDSupportMaterial"); + m_gearMgr->registerSimpleMaterial(VXDSupportMaterial); + + if (vtxData->rOuterShell>vtxData->rInnerShell) { + dd4hep::rec::Vector3D a1( vtxData->rInnerShell, 0, 2.*dd4hep::mm); + dd4hep::rec::Vector3D b1( vtxData->rOuterShell, 0, 2.*dd4hep::mm); + gear::SimpleMaterialImpl* VXDShellMaterial = CreateGearMaterial(a1, b1, "VXDShellMaterial"); + m_gearMgr->registerSimpleMaterial(VXDShellMaterial); } } + return StatusCode::SUCCESS; +} - const std::map<std::string, dd4hep::DetElement>& components = vxd.children(); - for(std::map<std::string, dd4hep::DetElement>::const_iterator it=components.begin();it!=components.end();it++){ - dd4hep::DetElement component = it->second; - dd4hep::Volume vol = component.volume(); - dd4hep::PlacedVolume phys = component.placement(); - TGeoMaterial* mat = vol->GetMaterial(); - const TGeoShape* shape = vol->GetShape(); - const dd4hep::PlacedVolumeExtension::VolIDs& ids = phys.volIDs(); - if(vol.isSensitive()&&shape->IsA()==TGeoBBox::Class()){ - int iLayer = ids.find("layer")->second; - int iModule = ids.find("module")->second; - int iSide = ids.find("side")->second; - if(iModule==0&&iLayer==helpCount+1){ - helpCount++; - helpSensitives.push_back(thisLadder); - helpLadders.push_back(thisLadder); - helpNumberLadders.push_back(nLadders); - helpPhi0.push_back(phi0); - nLadders = 0; - thisLadder.length = 0; - } - if(iLayer == helpCount){ - if(iModule == 0){ - const TGeoBBox* box = (const TGeoBBox*) shape; - double width = box->GetDX()*CLHEP::cm; - double length = box->GetDY()*CLHEP::cm; - double thickness = box->GetDZ()*CLHEP::cm; - TGeoMatrix* matrix = phys->GetMatrix(); - const double* pos = matrix->GetTranslation(); - const double* rot_data = matrix->GetRotationMatrix(); - TGeoRotation rot; - rot.SetMatrix(rot_data); - double theta,phi,psi; - rot.GetAngles(phi,theta,psi); - phi *= deg_to_rad; - theta *= deg_to_rad; - psi *= deg_to_rad; - phi0 = -dd4hep::halfpi+phi; - double distance = fabs(cos(phi0)*sin(theta)*pos[0]+sin(phi0)*sin(theta)*pos[1]+cos(theta)*pos[2]); - double offset = sqrt(pos[0]*pos[0]+pos[1]*pos[1]-distance*distance)*pos[0]/fabs(pos[0])*CLHEP::cm; - distance *= CLHEP::cm; - distance -= thickness; - double radL = mat->GetRadLen()*CLHEP::cm; - //info() << " -> " << helpCount << ": " << distance << " " << cos(atan2(pos[1],pos[0])-phi)*sqrt(pos[0]*pos[0]+pos[1]*pos[1]) << endmsg; - thisLadder.distance = distance; - thisLadder.offset = offset; - thisLadder.thickness = thickness; - thisLadder.length += length; - thisLadder.width = width; - thisLadder.radLength = radL; - thisLadder.z = pos[2]*CLHEP::cm; - } - if(iModule==nLadders) nLadders++; - } +StatusCode GearSvc::convertComposite(dd4hep::DetElement& vtx){ + dd4hep::rec::CompositeData* vtxData = nullptr; + try{ + vtxData = vtx.extension<dd4hep::rec::CompositeData>(); + } + catch(std::runtime_error& e){ + warning() << e.what() << " " << vtxData << endmsg; + return StatusCode::RECOVERABLE; + } + if(vtxData){ + int vtxType = gear::ZPlanarParameters::CMOS; + gear::ZPlanarParametersImpl* gearVTX = new gear::ZPlanarParametersImpl(vtxType, vtxData->rInnerShell/dd4hep::mm, vtxData->rOuterShell/dd4hep::mm, + vtxData->zHalfShell/dd4hep::mm, vtxData->gapShell/dd4hep::mm, 0.); + for (unsigned i=0,n=vtxData->layersPlanar.size(); i<n; ++i) { + const dd4hep::rec::ZPlanarData::LayerLayout& l = vtxData->layersPlanar[i]; + // FIXME set rad lengths to 0 -> need to get from dd4hep .... + gearVTX->addLayer(l.ladderNumber, l.phi0, + l.distanceSupport/dd4hep::mm, l.offsetSupport/dd4hep::mm, l.thicknessSupport/dd4hep::mm, + l.zHalfSupport/dd4hep::mm, l.widthSupport/dd4hep::mm, 0., + l.distanceSensitive/dd4hep::mm, l.offsetSensitive/dd4hep::mm, l.thicknessSensitive/dd4hep::mm, + l.zHalfSensitive/dd4hep::mm, l.widthSensitive/dd4hep::mm, 0.); } - else if(it->first=="VXD_support"){ - helpCount++; - helpSensitives.push_back(thisLadder); - helpLadders.push_back(thisLadder); - helpNumberLadders.push_back(nLadders); - helpPhi0.push_back(phi0); - nLadders = 0; - if(vol->GetNdaughters()==0) error() << "!!!!!!!!!" << endmsg; - - int nFlexCable = 0, nFoamSpacer=0, nMetalTraces=0; - int currentLayer = -1; - double tFlexCable=0, tFoamSpacer=0, tMetalTraces=0; - double radLFlexCable=0, radLFoamSpacer=0, radLMetalTraces=0; - double intLFlexCable=0, intLFoamSpacer=0, intLMetalTraces=0; - double dFlexCable=0, dFoamSpacer=0, dMetalTraces=0; - double metalZeff=0, metalZAeff=0, foamZeff=0, foamZAeff=0, flexZeff=0, flexZAeff=0; - for(int i=0;i<vol->GetNdaughters();i++){ - TGeoNode* daughter = vol->GetNode(i); - TGeoMaterial* matDaughter = daughter->GetMedium()->GetMaterial(); - const TGeoShape* shape_sup = daughter->GetVolume()->GetShape(); - TGeoMatrix* matrix = daughter->GetMatrix(); - const double* pos = matrix->GetTranslation(); - const double* rot_data = matrix->GetRotationMatrix(); - TGeoRotation rot; - rot.SetMatrix(rot_data); - double theta,phi,psi; - rot.GetAngles(phi,theta,psi); - phi *= deg_to_rad; - theta *= deg_to_rad; - psi *= deg_to_rad; - phi0 = -CLHEP::halfpi+phi; - std::string phy_name = daughter->GetName(); - if(phy_name.find("FoamSpacer")==-1&&phy_name.find("FlexCable")==-1&&phy_name.find("MetalTraces")==-1){ - //info() << phy_name <<endmsg; - continue; - } - int iLayer = atoi(phy_name.substr(phy_name.find("_")+1,2).c_str()); - if(iLayer!=currentLayer){ - //info() << tFoamSpacer << "," << tFlexCable << "," << tMetalTraces << endmsg; - helpLadders[currentLayer].thickness = tFoamSpacer+tFlexCable+tMetalTraces; - helpLadders[currentLayer].radLength = helpLadders[currentLayer].thickness / (tFoamSpacer/radLFoamSpacer+tFlexCable/radLFlexCable+tMetalTraces/radLMetalTraces); - nFlexCable = 0; - nFoamSpacer=0; - nMetalTraces=0; - currentLayer=iLayer; - } - if(shape_sup->IsA()==TGeoBBox::Class()&&(nFoamSpacer==0||nFlexCable==0||nMetalTraces==0)){ - const TGeoBBox* box = (const TGeoBBox*) shape_sup; - double distance = fabs(cos(phi0)*sin(theta)*pos[0]+sin(phi0)*sin(theta)*pos[1]+cos(theta)*pos[2]); - double offset = sqrt(pos[0]*pos[0]+pos[1]*pos[1]-distance*distance)*pos[0]/fabs(pos[0])*CLHEP::cm; - distance -= box->GetDZ(); - distance *= CLHEP::cm; - if(helpLadders[iLayer].distance == helpSensitives[iLayer].distance) helpLadders[iLayer].distance = distance; - else helpLadders[iLayer].distance = TMath::Min(helpLadders[iLayer].distance, distance); - if(phy_name.find("FoamSpacer")!=-1&&nFoamSpacer==0){ - helpLadders[iLayer].offset = offset; - tFoamSpacer = box->GetDZ()*CLHEP::cm; - radLFoamSpacer = matDaughter->GetRadLen()*CLHEP::cm; - intLFoamSpacer = matDaughter->GetIntLen()*CLHEP::cm; - dFoamSpacer = matDaughter->GetDensity(); - double totalA = 0, Zeff = 0, ZAeff = 0; - for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){ - totalA += matDaughter->GetElement(iEle)->A(); - } - for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){ - double A, Z, w; - // by fucd: w!=A/totalA, strange! to fix - matDaughter->GetElementProp(A,Z,w,iEle); - Zeff += Z*w; - ZAeff += Z/A*w; - //info() << std::setprecision(16) << Z << " " << A << " " << A/totalA << " " << w << endmsg; - } - foamZeff = Zeff; - foamZAeff = ZAeff; - nFoamSpacer++; - } - if(phy_name.find("FlexCable")!=-1&&nFlexCable==0){ - helpLadders[iLayer].width = box->GetDX()*CLHEP::cm; - helpLadders[iLayer].length = box->GetDY()*CLHEP::cm-beryllium_ladder_block_length*2-end_electronics_half_z*2; - tFlexCable = box->GetDZ()*CLHEP::cm; - radLFlexCable = matDaughter->GetRadLen()*CLHEP::cm; - intLFlexCable = matDaughter->GetIntLen()*CLHEP::cm; - dFlexCable = matDaughter->GetDensity(); - double Zeff = 0, ZAeff = 0; - for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){ - double A, Z, w; - matDaughter->GetElementProp(A,Z,w,iEle); - Zeff += Z*w; - ZAeff += Z/A*w; - //std::cout << std::setprecision(16) << Z << " " << A << " " << w << std::endl; - } - flexZeff = Zeff; - flexZAeff = ZAeff; - nFlexCable++; - } - if(phy_name.find("MetalTraces")!=-1&&nMetalTraces==0){ - tMetalTraces = box->GetDZ()*CLHEP::cm; - radLMetalTraces = matDaughter->GetRadLen()*CLHEP::cm; - intLMetalTraces = matDaughter->GetIntLen()*CLHEP::cm; - dMetalTraces = matDaughter->GetDensity(); - double totalA = 0, Zeff = 0, ZAeff = 0; - for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){ - totalA += matDaughter->GetElement(iEle)->A(); - } - for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){ - double A, Z, w; - matDaughter->GetElementProp(A,Z,w,iEle); - Zeff += Z*w; - ZAeff += Z/A*w; - //info() << Z << " " << A << " " << w << endmsg; - } - metalZeff = Zeff; - metalZAeff = ZAeff; - nMetalTraces++; - } - } - } - { - //info() << tFoamSpacer << "," << tFlexCable << "," << tMetalTraces << endmsg; - double tSupport = tMetalTraces + tFoamSpacer + tFlexCable; - helpLadders[currentLayer].thickness = tSupport; - helpLadders[currentLayer].radLength = helpLadders[currentLayer].thickness / (tFoamSpacer/radLFoamSpacer+tFlexCable/radLFlexCable+tMetalTraces/radLMetalTraces); - nFlexCable = 0; - nFoamSpacer=0; - nMetalTraces=0; - - //calculations of thickness fractions of each layer of the support - double metalTF = tMetalTraces / tSupport; - double foamTF = tFoamSpacer / tSupport; - double flexTF = tFlexCable / tSupport; - //info() << foamTF << "," << flexTF << "," << metalTF << endmsg; - //info() << dFoamSpacer/(CLHEP::kg/CLHEP::cm3) << "," << dFlexCable/(CLHEP::kg/CLHEP::cm3) << "," << dMetalTraces/(CLHEP::kg/CLHEP::cm3) << endmsg; - //info() << foamZeff << " " << flexZeff << " " << metalZeff << endmsg; - //info() << foamZAeff << " " << flexZAeff << " " << metalZAeff << endmsg; - double elemVol = 1*CLHEP::cm3; - double VXDSupportMass = (foamTF*dFoamSpacer + flexTF*dFlexCable + metalTF*dMetalTraces)*elemVol; - VXDSupportDensity = VXDSupportMass/elemVol; - double foamFM = 100. * ((foamTF*(elemVol)*dFoamSpacer) / VXDSupportMass) ; - double kaptonFM = 100. * ((flexTF*(elemVol)*dFlexCable) / VXDSupportMass) ; - double metalFM = 100. * ((metalTF*(elemVol)*dMetalTraces) / VXDSupportMass) ; - - VXDSupportRadLen = helpLadders[currentLayer].radLength; - - VXDSupportZeff = (metalFM/100.)*metalZeff + (kaptonFM/100.)*flexZeff + (foamFM/100.)*foamZeff; - double VXDSupportZAeff = (metalFM/100.)*metalZAeff + (kaptonFM/100.)*flexZAeff + (foamFM/100.)*foamZAeff; - VXDSupportAeff = VXDSupportZeff / VXDSupportZAeff; - double VXDSupportIntLength = 1. / ((metalTF/intLMetalTraces) + (flexTF/intLFlexCable) + (foamTF/intLFoamSpacer)); - VXDSupportDensity = VXDSupportDensity*(CLHEP::g/CLHEP::cm3)/(CLHEP::kg/CLHEP::m3); - //info() << "fucd: " << VXDSupportZeff << " " << VXDSupportAeff << " " << VXDSupportRadLen << " " << VXDSupportIntLength << " " << VXDSupportDensity << endmsg; - //info() << intLMetalTraces << " " << intLFlexCable << " " << intLFoamSpacer <<endmsg; - } + + std::vector<int> ids; + std::vector<double> zhalfs, rsens, tsens, rsups, tsups, phi0s, rgaps, dphis; + + for (unsigned i=0,n=vtxData->layersBent.size(); i<n; ++i) { + const dd4hep::rec::CylindricalData::LayerLayout& l = vtxData->layersBent[i]; + + ids.push_back(l.id); + zhalfs.push_back(l.zHalf/dd4hep::mm); + rsens.push_back(l.radiusSensitive/dd4hep::mm); + tsens.push_back(l.thicknessSensitive/dd4hep::mm); + rsups.push_back(l.radiusSupport/dd4hep::mm); + tsups.push_back(l.thicknessSupport/dd4hep::mm); + phi0s.push_back(l.phi0); + rgaps.push_back(l.rgap/dd4hep::mm); + dphis.push_back(l.dphi); } - //info() << it->first << endmsg; - } - if(end_electronics_half_z>0 && side_band_electronics_width==0) type = gear::ZPlanarParametersImpl::CCD ; - if(side_band_electronics_width>0 && end_electronics_half_z==0 ) type = gear::ZPlanarParametersImpl::CMOS ; - if(side_band_electronics_width>0 && end_electronics_half_z>0) type = gear::ZPlanarParametersImpl::HYBRID ; - gear::ZPlanarParametersImpl* vxdParameters = new gear::ZPlanarParametersImpl(type, shellInnerRadius, shellOuterRadius, shellHalfLength, gap, shellRadLength); - // by fucd: debug info, if validated enough, merge them in future - info() << "=====================from convertor==============================" << endmsg; - info() << type << " " << shellInnerRadius << " " << shellOuterRadius << " " << shellHalfLength << " " << gap << " " << shellRadLength << endmsg; - for(int i=0;i<helpCount;i++){ - vxdParameters->addLayer(helpNumberLadders[i] , helpPhi0[i] , - helpLadders[i].distance , helpLadders[i].offset, helpLadders[i].thickness*2 , - helpLadders[i].length , helpLadders[i].width*2 , helpLadders[i].radLength , - helpSensitives[i].distance, helpSensitives[i].offset , helpSensitives[i].thickness*2 , - helpSensitives[i].length , helpSensitives[i].width*2 , helpSensitives[i].radLength ) ; - info() << "fucd " << i << ": " << helpNumberLadders[i] << ", " << helpPhi0[i] << ", " - << helpLadders[i].distance << ", " << helpLadders[i].offset << ", " << helpLadders[i].thickness*2 << ", " << helpLadders[i].length << ", " - << helpLadders[i].width*2 << ", " << helpLadders[i].radLength << ", " << helpSensitives[i].distance << ", " << helpSensitives[i].offset << ", " - << helpSensitives[i].thickness*2 << ", " << helpSensitives[i].length << ", " << helpSensitives[i].width*2 << ", " << helpSensitives[i].radLength << endmsg; - } - m_gearMgr->setVXDParameters(vxdParameters) ; - if(rAlu!=0){ - // rAlu=0, denote no cryostat - gear::GearParametersImpl* gearParameters = new gear::GearParametersImpl; - //CryostatAlRadius, CryostatAlThickness, CryostatAlInnerR, CryostatAlZEndCap, CryostatAlHalfZ - gearParameters->setDoubleVal("CryostatAlRadius", rAlu); - gearParameters->setDoubleVal("CryostatAlThickness", drAlu); - gearParameters->setDoubleVal("CryostatAlInnerR", rInner); - gearParameters->setDoubleVal("CryostatAlZEndCap", aluEndcapZ = dzSty+drSty+drAlu/2); - gearParameters->setDoubleVal("CryostatAlHalfZ", dzSty+drSty); - m_gearMgr->setGearParameters("VXDInfra", gearParameters); - } - dd4hep::rec::MaterialManager matMgr( dd4hep::Detector::getInstance().world().volume() ) ; - const gear::VXDLayerLayout& l = m_gearMgr->getVXDParameters().getVXDLayerLayout(); - dd4hep::rec::Vector3D a( l.getSensitiveDistance(0)+l.getSensitiveThickness(0), l.getPhi0(0), 0. , dd4hep::rec::Vector3D::cylindrical ) ; - dd4hep::rec::Vector3D b( l.getLadderDistance(0)+l.getLadderThickness(0), l.getPhi0(0), 0. , dd4hep::rec::Vector3D::cylindrical ) ; - const dd4hep::rec::MaterialVec& materials = matMgr.materialsBetween( a , b ) ; - dd4hep::rec::MaterialData mat = ( materials.size() > 1 ? matMgr.createAveragedMaterial( materials ) : materials[0].first ) ; + gearVTX->setIntVals("VTXLayerIds", ids); + gearVTX->setDoubleVals("VTXLayerHalfLengths", zhalfs); + gearVTX->setDoubleVals("VTXLayerSensitiveRadius", rsens); + gearVTX->setDoubleVals("VTXLayerSensitiveThickness", tsens); + gearVTX->setDoubleVals("VTXLayerSupperRadius", rsups); + gearVTX->setDoubleVals("VTXLayerSupperThickness", tsups); + gearVTX->setDoubleVals("VTXLayerPhi0", phi0s); + gearVTX->setDoubleVals("VTXLayerRadialGap", rgaps); + gearVTX->setDoubleVals("VTXLayerDeltaPhi", dphis); + + m_gearMgr->setVXDParameters(gearVTX); + + const dd4hep::rec::CylindricalData::LayerLayout& l = vtxData->layersBent[0] ; + dd4hep::rec::Vector3D a(l.radiusSupport, l.phi0, 0., dd4hep::rec::Vector3D::cylindrical); + dd4hep::rec::Vector3D b(l.radiusSupport + l.thicknessSupport, l.phi0, 0., dd4hep::rec::Vector3D::cylindrical); + gear::SimpleMaterialImpl* VXDSupportMaterial = CreateGearMaterial(a, b, "VXDSupportMaterial"); + m_gearMgr->registerSimpleMaterial(VXDSupportMaterial); - std::cout << " ####### found materials between points : " << a << " and " << b << " : " ; - for( unsigned i=0,n=materials.size();i<n;++i){ - std::cout << materials[i].first.name() << "[" << materials[i].second << "], " ; + if (vtxData->rOuterShell>vtxData->rInnerShell) { + dd4hep::rec::Vector3D a1( vtxData->rInnerShell, 0, 2.*dd4hep::mm); + dd4hep::rec::Vector3D b1( vtxData->rOuterShell, 0, 2.*dd4hep::mm); + gear::SimpleMaterialImpl* VXDShellMaterial = CreateGearMaterial(a1, b1, "VXDShellMaterial"); + m_gearMgr->registerSimpleMaterial(VXDShellMaterial); + } } - std::cout << std::endl ; - std::cout << " averaged material : " << mat << std::endl ; - gear::SimpleMaterialImpl* VXDSupportMaterial = new gear::SimpleMaterialImpl("VXDSupportMaterial", mat.A(), mat.Z(), - mat.density()/(dd4hep::kg/(dd4hep::g*dd4hep::m3)), - mat.radiationLength()/dd4hep::mm, - mat.interactionLength()/dd4hep::mm); - - //effective A different with what in Mokka, fix them as Mokka's - gear::SimpleMaterialImpl* VXDFoamShellMaterial_old = new gear::SimpleMaterialImpl("VXDFoamShellMaterial_old", 1.043890843e+01, 5.612886646e+00, 2.500000000e+01, 1.751650267e+04, 0); - m_gearMgr->registerSimpleMaterial(VXDFoamShellMaterial_old); - gear::SimpleMaterialImpl* VXDFoamShellMaterial = new gear::SimpleMaterialImpl("VXDFoamShellMaterial", styAeff, styZeff, styDensity*(CLHEP::g/CLHEP::cm3)/(CLHEP::kg/CLHEP::m3), - styRadLen, styIntLen); - m_gearMgr->registerSimpleMaterial(VXDFoamShellMaterial); - m_gearMgr->registerSimpleMaterial(VXDSupportMaterial); - info() << "cryostat parameters: " << rAlu << " " << drAlu << " " << rInner << " " << aluEndcapZ << " " << aluHalfZ << endmsg; - - return sc; + return StatusCode::SUCCESS; } StatusCode GearSvc::convertFTD(dd4hep::DetElement& ftd){ @@ -799,8 +483,9 @@ StatusCode GearSvc::convertTPC(dd4hep::DetElement& tpc){ gear::PadRowLayout2D *padLayout = new gear::FixedPadSizeDiskLayout(tpcData->rMinReadout*CLHEP::cm, tpcData->rMaxReadout*CLHEP::cm, tpcData->padHeight*CLHEP::cm, tpcData->padWidth*CLHEP::cm, tpcData->maxRow, 0.0); tpcParameters->setPadLayout(padLayout); + always() << tpcParameters->getPlaneExtent()[0] << " " << tpcParameters->getPlaneExtent()[1] << " " << tpcData->rMinReadout*CLHEP::cm << endmsg; tpcParameters->setMaxDriftLength(tpcData->driftLength*CLHEP::cm); - tpcParameters->setDriftVelocity( 0.0); // SJA: not set in Mokka so set to 0.0 + tpcParameters->setDriftVelocity( 0.0); // SJA: not set in Mokka so set to 0.0 tpcParameters->setReadoutFrequency( 0.0); tpcParameters->setDoubleVal( "tpcOuterRadius" , tpcData->rMax*CLHEP::cm ) ; tpcParameters->setDoubleVal( "tpcInnerRadius", tpcData->rMin*CLHEP::cm ) ; @@ -824,44 +509,19 @@ StatusCode GearSvc::convertTPC(dd4hep::DetElement& tpc){ tpcParameters->setDoubleVal( "tpcEndplateZmin", supportData->sections[1].zPos*CLHEP::cm + 1*CLHEP::um ); tpcParameters->setDoubleVal( "tpcEndplateZmax", supportData->sections[2].zPos*CLHEP::cm ); - dd4hep::rec::MaterialManager matMgr( dd4hep::Detector::getInstance().world().volume() ); double x = (supportData->sections[0].rInner + supportData->sections[0].rOuter)/2.; // Readout { dd4hep::rec::Vector3D a(x, 0, supportData->sections[0].zPos); dd4hep::rec::Vector3D b(x, 0, supportData->sections[1].zPos); - const dd4hep::rec::MaterialVec& materials = matMgr.materialsBetween(a, b); - dd4hep::rec::MaterialData mat = (materials.size()>1) ? matMgr.createAveragedMaterial(materials) : materials[0].first; - - std::cout << " ####### found materials between points : " << a << " and " << b << " : " ; - for( unsigned i=0,n=materials.size();i<n;++i){ - std::cout << materials[i].first.name() << "[" << materials[i].second << "], " ; - } - std::cout << std::endl ; - std::cout << " averaged material : " << mat << std::endl ; - gear::SimpleMaterialImpl* TPCReadoutMaterial = new gear::SimpleMaterialImpl("TPCReadoutMaterial", mat.A(), mat.Z(), - mat.density()/(dd4hep::kg/(dd4hep::g*dd4hep::m3)), - mat.radiationLength()/dd4hep::mm, - mat.interactionLength()/dd4hep::mm); + gear::SimpleMaterialImpl* TPCReadoutMaterial = CreateGearMaterial(a, b, "TPCReadoutMaterial"); m_gearMgr->registerSimpleMaterial(TPCReadoutMaterial); } // Endplate { dd4hep::rec::Vector3D a(x, 0, supportData->sections[1].zPos); dd4hep::rec::Vector3D b(x, 0, supportData->sections[2].zPos); - const dd4hep::rec::MaterialVec& materials = matMgr.materialsBetween(a, b); - dd4hep::rec::MaterialData mat = (materials.size()>1) ? matMgr.createAveragedMaterial(materials) : materials[0].first; - - std::cout << " ####### found materials between points : " << a << " and " << b << " : " ; - for( unsigned i=0,n=materials.size();i<n;++i){ - std::cout << materials[i].first.name() << "[" << materials[i].second << "], " ; - } - std::cout << std::endl ; - std::cout << " averaged material : " << mat << std::endl ; - gear::SimpleMaterialImpl* TPCEndplateMaterial = new gear::SimpleMaterialImpl("TPCEndplateMaterial", mat.A(), mat.Z(), - mat.density()/(dd4hep::kg/(dd4hep::g*dd4hep::m3)), - mat.radiationLength()/dd4hep::mm, - mat.interactionLength()/dd4hep::mm); + gear::SimpleMaterialImpl* TPCEndplateMaterial = CreateGearMaterial(a, b, "TPCEndplateMaterial"); m_gearMgr->registerSimpleMaterial(TPCEndplateMaterial); } @@ -1062,12 +722,10 @@ StatusCode GearSvc::convertCal(dd4hep::DetElement& cal) { std::string name = cal.name(); dd4hep::rec::LayeredCalorimeterData* calData = nullptr; - bool extensionDataValid = true; try{ calData = cal.extension<dd4hep::rec::LayeredCalorimeterData>(); } catch(std::runtime_error& e){ - extensionDataValid = false; info() << e.what() << " " << calData << endmsg; } if(calData){ @@ -1147,3 +805,23 @@ TGeoNode* GearSvc::FindNode(TGeoNode* mother, char* name) { } return next; } + +gear::SimpleMaterialImpl* GearSvc::CreateGearMaterial(const dd4hep::rec::Vector3D& a, const dd4hep::rec::Vector3D& b, + const std::string name) { + // TODO: move to GeomSvc + dd4hep::rec::MaterialManager matMgr( dd4hep::Detector::getInstance().world().volume() ); + + const dd4hep::rec::MaterialVec& materials = matMgr.materialsBetween(a, b); + dd4hep::rec::MaterialData mat = (materials.size() > 1) ? matMgr.createAveragedMaterial(materials) : materials[0].first; + + debug() << " ####### found materials between points : " << a << " and " << b << " ######" << endmsg; + for (unsigned i=0,n=materials.size(); i<n; ++i) { + debug() << materials[i].first.name() << " [" << materials[i].second << "]" << endmsg; + } + debug() << " averaged material : " << mat << endmsg; + gear::SimpleMaterialImpl* gearMaterial = new gear::SimpleMaterialImpl(name.c_str(), mat.A(), mat.Z(), + mat.density()/(dd4hep::kg/(dd4hep::g*dd4hep::m3)), + mat.radiationLength()/dd4hep::mm, + mat.interactionLength()/dd4hep::mm); + return gearMaterial; +} diff --git a/Service/GearSvc/src/GearSvc.h b/Service/GearSvc/src/GearSvc.h index 3e2e250b9daf25feba27d6972a95f61f56b45ea0..38e831de25a31fd15fde20990cfc6ff412bd6695 100644 --- a/Service/GearSvc/src/GearSvc.h +++ b/Service/GearSvc/src/GearSvc.h @@ -4,6 +4,8 @@ #include "GearSvc/IGearSvc.h" #include <GaudiKernel/Service.h> #include "DD4hep/Detector.h" +#include "DDRec/Vector3D.h" +#include "gearimpl/SimpleMaterialImpl.h" class TGeoNode; class GearSvc : public extends<Service, IGearSvc> @@ -20,6 +22,8 @@ class GearSvc : public extends<Service, IGearSvc> private: StatusCode convertBeamPipe(dd4hep::DetElement& pipe); StatusCode convertVXD(dd4hep::DetElement& vxd); + StatusCode convertStitching(dd4hep::DetElement& vtx); + StatusCode convertComposite(dd4hep::DetElement& vtx); StatusCode convertSIT(dd4hep::DetElement& sit); StatusCode convertTPC(dd4hep::DetElement& tpc); StatusCode convertDC (dd4hep::DetElement& dc); @@ -27,6 +31,7 @@ class GearSvc : public extends<Service, IGearSvc> StatusCode convertFTD(dd4hep::DetElement& ftd); StatusCode convertCal(dd4hep::DetElement& cal); TGeoNode* FindNode(TGeoNode* mother, char* name); + gear::SimpleMaterialImpl* CreateGearMaterial(const dd4hep::rec::Vector3D& a, const dd4hep::rec::Vector3D& b, const std::string name); Gaudi::Property<std::string> m_gearFile{this, "GearXMLFile", ""}; Gaudi::Property<float> m_field{this, "MagneticField", 0}; diff --git a/Service/TrackSystemSvc/src/MarlinKalTest.cc b/Service/TrackSystemSvc/src/MarlinKalTest.cc index 5931143ffd6bd50ba136097eb32f5bf95abf62e2..302675bf183d95f2cd05d7f3b0c793b56f6483b9 100644 --- a/Service/TrackSystemSvc/src/MarlinKalTest.cc +++ b/Service/TrackSystemSvc/src/MarlinKalTest.cc @@ -9,6 +9,7 @@ #include "kaldet/ILDSupportKalDetector.h" #include "kaldet/ILDVXDKalDetector.h" +#include "kaldet/CEPCVTXKalDetector.h" #include "kaldet/ILDSITKalDetector.h" #include "kaldet/ILDSITCylinderKalDetector.h" #include "kaldet/ILDSETKalDetector.h" @@ -114,7 +115,8 @@ namespace MarlinTrk{ } try{ - ILDVXDKalDetector* vxddet = new ILDVXDKalDetector( *_gearMgr, _geoSvc ) ; + //ILDVXDKalDetector* vxddet = new ILDVXDKalDetector( *_gearMgr, _geoSvc ) ; + CEPCVTXKalDetector* vxddet = new CEPCVTXKalDetector(*_gearMgr, _geoSvc); // store the measurement layer id's for the active layers this->storeActiveMeasurementModuleIDs(vxddet); _det->Install( *vxddet ) ; diff --git a/Utilities/KalDet/include/kaldet/CEPCVTXKalDetector.h b/Utilities/KalDet/include/kaldet/CEPCVTXKalDetector.h new file mode 100644 index 0000000000000000000000000000000000000000..7184f78c5c8b25ab1b2373fa3be15cf7f36a5b38 --- /dev/null +++ b/Utilities/KalDet/include/kaldet/CEPCVTXKalDetector.h @@ -0,0 +1,65 @@ +#ifndef __CEPCVTXKALDETECTOR__ +#define __CEPCVTXKALDETECTOR__ + +/** Ladder based VXD to be used for ILD DBD studies + * + * @author S.Aplin DESY + */ + +#include "kaltest/TVKalDetector.h" + +#include "TMath.h" + +class TNode; +class IGeomSvc; + +namespace gear{ + class GearMgr ; +} + + +class CEPCVTXKalDetector : public TVKalDetector { + +public: + + CEPCVTXKalDetector( const gear::GearMgr& gearMgr, IGeomSvc* geoSvc); + + +private: + + void setupGearGeom( const gear::GearMgr& gearMgr ); + void setupGearGeom( IGeomSvc* geoSvc) ; + + int _nLayers[2]; + double _bZ; + + double _relative_position_of_measurement_surface; + + struct VXD_Layer { + int nLadders; + double phi0; + double dphi; + double senRMin; + double supRMin; + double length; + double width; + double offset; + double senThickness; + double supThickness; + }; + std::vector<VXD_Layer> _VXDgeo; + + struct STT_Layer { + int id; + double length; + double senRMin; + double senThickness; + double supRMin; + double supThickness; + double phi0; + double rgap; + double dphi; + }; + std::vector<STT_Layer> _STTgeo; +}; +#endif diff --git a/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.cc b/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.cc new file mode 100644 index 0000000000000000000000000000000000000000..3b7528b5258cf58a02fb428c4839c9f59e057b5f --- /dev/null +++ b/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.cc @@ -0,0 +1,383 @@ + +#include "CEPCVTXKalDetector.h" + +#include "MaterialDataBase.h" + +#include "ILDParallelPlanarMeasLayer.h" +#include "ILDCylinderMeasLayer.h" +#include "ILDDiscMeasLayer.h" + +#include <UTIL/BitField64.h> +#include <UTIL/ILDConf.h> + +#include "DetInterface/IGeomSvc.h" +#include "DD4hep/Detector.h" +#include "DDRec/DetectorData.h" +#include "CLHEP/Units/SystemOfUnits.h" +#include "DD4hep/DD4hepUnits.h" + +#include <gear/GEAR.h> +#include "gear/BField.h" +#include <gearimpl/ZPlanarParametersImpl.h> +#include <gear/VXDParameters.h> +#include <gear/VXDLayerLayout.h> +#include "gearimpl/Util.h" +#include "DetInterface/IGeomSvc.h" + +#include "TMath.h" + +#include "math.h" +#include <sstream> + +// #include "streamlog/streamlog.h" + +CEPCVTXKalDetector::CEPCVTXKalDetector( const gear::GearMgr& gearMgr, IGeomSvc* geoSvc ) +: TVKalDetector(300) // SJA:FIXME initial size, 300 looks reasonable for ILD, though this would be better stored as a const somewhere +{ + + // streamlog_out(DEBUG1) << "CEPCVTXKalDetector building VXD detector using GEAR " << std::endl ; + + MaterialDataBase::Instance().registerForService(gearMgr, geoSvc); + + TMaterial & air = *MaterialDataBase::Instance().getMaterial("air"); + TMaterial & silicon = *MaterialDataBase::Instance().getMaterial("silicon"); + TMaterial & carbon = *MaterialDataBase::Instance().getMaterial("VXDSupportMaterial"); + TMaterial & beryllium = *MaterialDataBase::Instance().getMaterial("beryllium"); + + // needed for cryostat + TMaterial & aluminium = *MaterialDataBase::Instance().getMaterial("aluminium"); + + if(geoSvc){ + this->setupGearGeom(geoSvc) ; + } + else{ + this->setupGearGeom(gearMgr) ; + } + + //--The Ladder structure (realistic ladder)-- + int nLadders; + + Bool_t active = true; + Bool_t dummy = false; + + static const double eps = 1e-6; + + UTIL::BitField64 encoder( lcio::ILDCellID0::encoder_string ) ; + + for (int layer=0; layer<_nLayers[0]; ++layer) { + + nLadders = _VXDgeo[layer].nLadders ; + + double phi0 = _VXDgeo[layer].phi0 ; + + double ladder_distance = _VXDgeo[layer].supRMin ; + double ladder_thickness = _VXDgeo[layer].supThickness ; + + double sensitive_distance = _VXDgeo[layer].senRMin ; + double sensitive_thickness = _VXDgeo[layer].senThickness ; + + double width = _VXDgeo[layer].width ; + double length = _VXDgeo[layer].length; + double offset = _VXDgeo[layer].offset; + + double pos_xi_nonoverlap_width = (2.0 * (( width / 2.0 ) - fabs(offset))); + + double currPhi; + double dphi = _VXDgeo[layer].dphi ; + + static const double z_offset = 0.0; // all VXD planes are centred at zero + + for (int ladder=0; ladder<nLadders; ++ladder) { + + currPhi = phi0 + (dphi * ladder); + + encoder.reset() ; // reset to 0 + + encoder[lcio::ILDCellID0::subdet] = lcio::ILDDetID::VXD ; + encoder[lcio::ILDCellID0::side] = 0 ; + encoder[lcio::ILDCellID0::layer] = layer + _nLayers[1]; + encoder[lcio::ILDCellID0::module] = ladder ; + encoder[lcio::ILDCellID0::sensor] = 0 ; + + int CellID = encoder.lowWord() ; + + // even layers have the senstive side facing the IP + if(layer%2 == 0 ){ // overlap section of ladder0 is defined after the last ladder, + + + double sen_front_sorting_policy = sensitive_distance + (4 * ladder+0) * eps ; + double measurement_plane_sorting_policy = sensitive_distance + (4 * ladder+1) * eps ; + double sen_back_sorting_policy = sensitive_distance + (4 * ladder+2) * eps ; + double sup_back_sorting_policy = sensitive_distance + (4 * ladder+3) * eps ; + + + if(ladder==0){ // bacause overlap section of ladder0 is further outer than the last ladder. + + // streamlog_out(DEBUG0) << "CEPCVTXKalDetector add surface with CellID = " + // << CellID + // << std::endl ; + + // non overlapping region + // air - sensitive boundary + Add(new ILDParallelPlanarMeasLayer(air, silicon, sensitive_distance, currPhi, _bZ, sen_front_sorting_policy, pos_xi_nonoverlap_width, length, 0.0, z_offset, offset, dummy,-1,"VXDSenFront_non_overlap_even" )) ; + + // measurement plane defined as the middle of the sensitive volume - unless "relative_position_of_measurement_surface" parameter given in GEAR + Add(new ILDParallelPlanarMeasLayer(silicon, silicon, sensitive_distance+sensitive_thickness*_relative_position_of_measurement_surface, currPhi, _bZ, measurement_plane_sorting_policy, pos_xi_nonoverlap_width, length, 0.0, z_offset, offset, active, CellID, "VXDMeasLayer_non_overlap_even" )) ; + + // sensitive - support boundary + Add(new ILDParallelPlanarMeasLayer(silicon, carbon, sensitive_distance+sensitive_thickness, currPhi, _bZ, sen_back_sorting_policy, pos_xi_nonoverlap_width, length, 0.0, z_offset, offset, dummy,-1,"VXDSenSuppportIntf_non_overlap_even" )) ; + + // support - air boundary + Add(new ILDParallelPlanarMeasLayer(carbon, air, ladder_distance+ladder_thickness, currPhi, _bZ, sup_back_sorting_policy, pos_xi_nonoverlap_width, length, 0.0, z_offset, offset, dummy,-1,"VXDSupRear_non_overlap_even" )) ; + + + // overlapping region + double overlap_region_width = width - pos_xi_nonoverlap_width ; + double overlap_region_offset = -(overlap_region_width/2.0) - (pos_xi_nonoverlap_width)/2.0 ; + + + // overlap sorting policy uses nLadders as the overlapping "ladder" is the order i.e. there will now be nLadders+1 + double overlap_front_sorting_policy = sensitive_distance + (4* nLadders+0) * eps; + double overlap_measurement_plane_sorting_policy = sensitive_distance + (4* nLadders+1) * eps; + double overlap_back_sorting_policy = sensitive_distance + (4* nLadders+2) * eps; + double overlap_sup_back_sorting_policy = sensitive_distance + (4* nLadders+3) * eps; + + // streamlog_out(DEBUG0) << "CEPCVTXKalDetector add surface with CellID = " + // << CellID + // << std::endl ; + + // air - sensitive boundary + Add(new ILDParallelPlanarMeasLayer(air, silicon, sensitive_distance, currPhi, _bZ, overlap_front_sorting_policy, overlap_region_width, length, overlap_region_offset, z_offset, offset, dummy,-1,"VXDSenFront_overlap_even")) ; + + // measurement plane defined as the middle of the sensitive volume - unless "relative_position_of_measurement_surface" parameter given in GEAR + Add(new ILDParallelPlanarMeasLayer(silicon, silicon, sensitive_distance+sensitive_thickness*_relative_position_of_measurement_surface, currPhi, _bZ, overlap_measurement_plane_sorting_policy, overlap_region_width, length, overlap_region_offset, z_offset, offset, active, CellID, "VXDMeasLayer_overlap_even" )) ; + + + // sensitive - support boundary + Add(new ILDParallelPlanarMeasLayer(silicon, carbon, sensitive_distance+sensitive_thickness, currPhi, _bZ, overlap_back_sorting_policy, overlap_region_width, length, overlap_region_offset, z_offset, offset, dummy,-1,"VXDSenSuppportIntf_overlap_even")) ; + + // support - air boundary + Add(new ILDParallelPlanarMeasLayer(carbon, air, ladder_distance+ladder_thickness, currPhi, _bZ, overlap_sup_back_sorting_policy, overlap_region_width, length, overlap_region_offset, z_offset, offset, dummy,-1,"VXDSupRear_overlap_even")) ; + + } + else{ + + // streamlog_out(DEBUG0) << "CEPCVTXKalDetector (ILDParallelPlanarMeasLayer) add surface with CellID = " + // << CellID + // << std::endl ; + + + // air - sensitive boundary + Add(new ILDParallelPlanarMeasLayer(air, silicon, sensitive_distance, currPhi, _bZ, sen_front_sorting_policy, width, length, offset, z_offset, offset, dummy,-1, "VXDSenFront_even")) ; + + + // measurement plane defined as the middle of the sensitive volume - unless "relative_position_of_measurement_surface" parameter given in GEAR - even layers face outwards ! + Add(new ILDParallelPlanarMeasLayer(silicon, silicon, sensitive_distance+sensitive_thickness*( 1.-_relative_position_of_measurement_surface ), currPhi, _bZ, measurement_plane_sorting_policy, width, length, offset, z_offset, offset, active, CellID, "VXDMeaslayer_even" )) ; + + // sensitive - support boundary + Add(new ILDParallelPlanarMeasLayer(silicon, carbon, sensitive_distance+sensitive_thickness, currPhi, _bZ, sen_back_sorting_policy, width, length, offset, z_offset, offset, dummy,-1,"VXDSenSuppportIntf_even" )) ; + + // support - air boundary + Add(new ILDParallelPlanarMeasLayer(carbon, air, ladder_distance+ladder_thickness, currPhi, _bZ, sup_back_sorting_policy, width, length, offset, z_offset, offset, dummy,-1,"VXDSupRear_even" )) ; + + } + } + else{ // counting from 0, odd numbered layers are placed with the support closer to the IP than the sensitive + + + + double sup_forward_sorting_policy = ladder_distance + (4 * ladder+0) * eps ; + double sup_back_sorting_policy = ladder_distance + (4 * ladder+1) * eps ; + double measurement_plane_sorting_policy = ladder_distance + (4 * ladder+2) * eps ; + double sen_back_sorting_policy = ladder_distance + (4 * ladder+3) * eps ; + + // streamlog_out(DEBUG0) << "CEPCVTXKalDetector (ILDPlanarMeasLayer) add surface with CellID = " + // << CellID + // << std::endl ; + + // air - support boundary + Add(new ILDParallelPlanarMeasLayer(air, carbon, ladder_distance, currPhi, _bZ, sup_forward_sorting_policy, width, length, offset, z_offset, offset, dummy,-1,"VXDSupFront_odd" )) ; + + // support - sensitive boundary + Add(new ILDParallelPlanarMeasLayer(carbon, silicon, (ladder_distance+ladder_thickness), currPhi, _bZ, sup_back_sorting_policy, width, length, offset, z_offset, offset, dummy,-1,"VXDSenSuppportIntf_odd")) ; + + // measurement plane defined as the middle of the sensitive volume + Add(new ILDParallelPlanarMeasLayer(silicon, silicon, (sensitive_distance+sensitive_thickness*0.5), currPhi, _bZ, measurement_plane_sorting_policy, width, length, offset, z_offset, offset, active, CellID, "VXDMeaslayer_odd")) ; + + // sensitive air - sensitive boundary + Add(new ILDParallelPlanarMeasLayer(silicon, air, (sensitive_distance+sensitive_thickness), currPhi, _bZ, sen_back_sorting_policy, width, length, offset, z_offset, offset, dummy,-1,"VXDSenRear_odd")) ; + + + } + } + } + + for (int layer=0; layer<_nLayers[1]; ++layer) { + //std::cout << "add stitching ... " << layer << std::endl; + double phi0 = _STTgeo[layer].phi0; + + double ladder_distance = _STTgeo[layer].supRMin; + double ladder_thickness = _STTgeo[layer].supThickness; + + double sensitive_distance = _STTgeo[layer].senRMin; + double sensitive_thickness = _STTgeo[layer].senThickness; + + //std::cout << "sensor radius: " << sensitive_distance << " sensitive_thickness: " << sensitive_thickness << std::endl; + //double width = _STTgeo[layer].width; + double halfz = _STTgeo[layer].length/2.0; + + double currPhi; + double dphi = _STTgeo[layer].dphi ; + + static const double z_offset = 0.0; // all VXD planes are centred at zero + double x0(0), y0(0), z0(0); + + for (int im=0; im<2; im++) { + currPhi = phi0 + (dphi * im); + + encoder.reset() ; // reset to 0 + encoder[lcio::ILDCellID0::subdet] = lcio::ILDDetID::VXD; + encoder[lcio::ILDCellID0::side] = 0; + encoder[lcio::ILDCellID0::layer] = layer; + encoder[lcio::ILDCellID0::module] = im; + encoder[lcio::ILDCellID0::sensor] = 0; + + int CellID = encoder.lowWord() ; + + if (im==1) { + sensitive_distance += _STTgeo[layer].rgap; + ladder_distance += _STTgeo[layer].rgap; + } + + // air - sensitive boundary + Add(new ILDCylinderMeasLayer(air, silicon, sensitive_distance, halfz, x0, y0, z0, _bZ, dummy, -1, "STTMeasL_0")); + // measurement plane defined as the middle of the sensitive volume + // - unless "relative_position_of_measurement_surface" parameter given in GEAR - even layers face outwards ! + Add(new ILDCylinderMeasLayer(silicon, silicon, sensitive_distance+sensitive_thickness*(1.-_relative_position_of_measurement_surface), + halfz, x0, y0, z0, _bZ, active, CellID, "STTMeaslayer_0")); + // sensitive - support boundary + Add(new ILDCylinderMeasLayer(silicon, carbon, sensitive_distance+sensitive_thickness, halfz, x0, y0, z0, _bZ, dummy, -1, "STTSenSuppportIntf_0")); + // support - air boundary + Add(new ILDCylinderMeasLayer(carbon, air, ladder_distance+ladder_thickness, halfz, x0, y0, z0, _bZ, dummy, -1, "STTSupRear_0" )) ; + } + } + + SetOwner(); +} + +void CEPCVTXKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){ + const gear::VXDParameters& pVXDDetMain = gearMgr.getVXDParameters(); + const gear::VXDLayerLayout& pVXDLayerLayout = pVXDDetMain.getVXDLayerLayout(); + + _bZ = gearMgr.getBField().at( gear::Vector3D( 0.,0.,0.) ).z() ; + + _nLayers[0] = pVXDLayerLayout.getNLayers(); + _VXDgeo.resize(_nLayers[0]); + + for( int layer=0; layer < _nLayers[0]; ++layer){ + _VXDgeo[layer].nLadders = pVXDLayerLayout.getNLadders(layer); + _VXDgeo[layer].phi0 = pVXDLayerLayout.getPhi0(layer); + _VXDgeo[layer].dphi = 2*M_PI / _VXDgeo[layer].nLadders; + _VXDgeo[layer].senRMin = pVXDLayerLayout.getSensitiveDistance(layer); + _VXDgeo[layer].supRMin = pVXDLayerLayout.getLadderDistance(layer); + _VXDgeo[layer].length = pVXDLayerLayout.getSensitiveLength(layer) * 2.0 ; // note: gear for historical reasons uses the halflength + _VXDgeo[layer].width = pVXDLayerLayout.getSensitiveWidth(layer); + _VXDgeo[layer].offset = pVXDLayerLayout.getSensitiveOffset(layer); + _VXDgeo[layer].senThickness = pVXDLayerLayout.getSensitiveThickness(layer); + _VXDgeo[layer].supThickness = pVXDLayerLayout.getLadderThickness(layer); + //std::cout << layer << ": " << _VXDgeo[layer].nLadders << " " << _VXDgeo[layer].phi0 << " " << _VXDgeo[layer].dphi << " " << _VXDgeo[layer].senRMin + // << " " << _VXDgeo[layer].supRMin << " " << _VXDgeo[layer].length << " " << _VXDgeo[layer].width << " " << _VXDgeo[layer].offset + // << " " << _VXDgeo[layer].senThickness << " " << _VXDgeo[layer].supThickness << std::endl; + } + + _relative_position_of_measurement_surface = 0.5 ; + + try { + _relative_position_of_measurement_surface = pVXDDetMain.getDoubleVal( "relative_position_of_measurement_surface" ); + } + catch (gear::UnknownParameterException& e) {} + + _nLayers[1] = 0; + try { + const std::vector<int> ids = pVXDDetMain.getIntVals("VTXLayerIds"); + const std::vector<double> zhalfs = pVXDDetMain.getDoubleVals("VTXLayerHalfLengths"); + const std::vector<double> rsens = pVXDDetMain.getDoubleVals("VTXLayerSensitiveRadius"); + const std::vector<double> tsens = pVXDDetMain.getDoubleVals("VTXLayerSensitiveThickness"); + const std::vector<double> rsups = pVXDDetMain.getDoubleVals("VTXLayerSupperRadius"); + const std::vector<double> tsups = pVXDDetMain.getDoubleVals("VTXLayerSupperThickness"); + const std::vector<double> phi0s = pVXDDetMain.getDoubleVals("VTXLayerPhi0"); + const std::vector<double> rgaps = pVXDDetMain.getDoubleVals("VTXLayerRadialGap"); + const std::vector<double> dphis = pVXDDetMain.getDoubleVals("VTXLayerDeltaPhi"); + + _nLayers[1] = ids.size(); + _STTgeo.resize(_nLayers[1]); + for (int ilayer=0; ilayer<_nLayers[1]; ilayer++) { + _STTgeo[ilayer].id = ids[ilayer]; + _STTgeo[ilayer].length = zhalfs[ilayer]*2; + _STTgeo[ilayer].senRMin = rsens[ilayer]; + _STTgeo[ilayer].senThickness = tsens[ilayer]; + _STTgeo[ilayer].supRMin = rsups[ilayer]; + _STTgeo[ilayer].supThickness = tsups[ilayer]; + _STTgeo[ilayer].phi0 = phi0s[ilayer]; + _STTgeo[ilayer].rgap = rgaps[ilayer]; + _STTgeo[ilayer].dphi = dphis[ilayer]; + } + } + catch (gear::UnknownParameterException& e) {} + + //std::cout << "planar: " << _nLayers[0] << " bent: " << _nLayers[1] << std::endl; +} + + +void CEPCVTXKalDetector::setupGearGeom( IGeomSvc* geoSvc){ + dd4hep::DetElement world = geoSvc->getDD4HepGeo(); + dd4hep::DetElement vxd; + const std::map<std::string, dd4hep::DetElement>& subs = world.children(); + for(std::map<std::string, dd4hep::DetElement>::const_iterator it=subs.begin();it!=subs.end();it++){ + if(it->first!="VXD") continue; + vxd = it->second; + } + dd4hep::rec::ZPlanarData* vxdData = nullptr; + try{ + vxdData = vxd.extension<dd4hep::rec::ZPlanarData>(); + } + catch(std::runtime_error& e){ + std::cout << e.what() << " " << vxdData << std::endl; + throw GaudiException(e.what(), "FATAL", StatusCode::FAILURE); + } + + const dd4hep::Direction& field = geoSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0)); + _bZ = field.z()/dd4hep::tesla; + + const std::vector<dd4hep::rec::ZPlanarData::LayerLayout>& layers = vxdData->layers; + + _nLayers[0] = layers.size(); + _VXDgeo.resize(_nLayers[0]); + + //SJA:FIXME: for now the support is taken as the same size the sensitive + // if this is not done then the exposed areas of the support would leave a carbon - air boundary, + // which if traversed in the reverse direction to the next boundary then the track will be propagated through carbon + // for a significant distance + + for( int layer=0; layer < _nLayers[0]; ++layer){ + const dd4hep::rec::ZPlanarData::LayerLayout& thisLayer = layers[layer]; + _VXDgeo[layer].nLadders = thisLayer.ladderNumber; + _VXDgeo[layer].phi0 = thisLayer.phi0; + _VXDgeo[layer].dphi = 2*M_PI / _VXDgeo[layer].nLadders; + _VXDgeo[layer].senRMin = thisLayer.distanceSensitive; + _VXDgeo[layer].supRMin = thisLayer.distanceSupport; + _VXDgeo[layer].length = thisLayer.zHalfSensitive * 2.0 ; // note: gear for historical reasons uses the halflength + _VXDgeo[layer].width = thisLayer.widthSensitive; + _VXDgeo[layer].offset = thisLayer.offsetSensitive; + _VXDgeo[layer].senThickness = thisLayer.thicknessSensitive; + _VXDgeo[layer].supThickness = thisLayer.thicknessSupport; + //std::cout << layer << ": " << _VXDgeo[layer].nLadders << " " << _VXDgeo[layer].phi0 << " " << _VXDgeo[layer].dphi << " " << _VXDgeo[layer].senRMin + // << " " << _VXDgeo[layer].supRMin << " " << _VXDgeo[layer].length << " " << _VXDgeo[layer].width << " " << _VXDgeo[layer].offset + // << " " << _VXDgeo[layer].senThickness << " " << _VXDgeo[layer].supThickness << std::endl; + } + // by default, we put the measurement surface in the middle of the sensitive + // layer, this can optionally be changed, e.g. in the case of the FPCCD where the + // epitaxial layer is 15 mu thick (in a 50 mu wafer) + _relative_position_of_measurement_surface = 0.5 ; + +} diff --git a/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.h b/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.h new file mode 100644 index 0000000000000000000000000000000000000000..7184f78c5c8b25ab1b2373fa3be15cf7f36a5b38 --- /dev/null +++ b/Utilities/KalDet/src/ild/vxd/CEPCVTXKalDetector.h @@ -0,0 +1,65 @@ +#ifndef __CEPCVTXKALDETECTOR__ +#define __CEPCVTXKALDETECTOR__ + +/** Ladder based VXD to be used for ILD DBD studies + * + * @author S.Aplin DESY + */ + +#include "kaltest/TVKalDetector.h" + +#include "TMath.h" + +class TNode; +class IGeomSvc; + +namespace gear{ + class GearMgr ; +} + + +class CEPCVTXKalDetector : public TVKalDetector { + +public: + + CEPCVTXKalDetector( const gear::GearMgr& gearMgr, IGeomSvc* geoSvc); + + +private: + + void setupGearGeom( const gear::GearMgr& gearMgr ); + void setupGearGeom( IGeomSvc* geoSvc) ; + + int _nLayers[2]; + double _bZ; + + double _relative_position_of_measurement_surface; + + struct VXD_Layer { + int nLadders; + double phi0; + double dphi; + double senRMin; + double supRMin; + double length; + double width; + double offset; + double senThickness; + double supThickness; + }; + std::vector<VXD_Layer> _VXDgeo; + + struct STT_Layer { + int id; + double length; + double senRMin; + double senThickness; + double supRMin; + double supThickness; + double phi0; + double rgap; + double dphi; + }; + std::vector<STT_Layer> _STTgeo; +}; +#endif